21 #include "libmesh/libmesh_config.h" 23 #ifdef LIBMESH_HAVE_LASPACK 25 #include "libmesh/laspack_matrix.h" 26 #include "libmesh/dense_matrix.h" 27 #include "libmesh/dof_map.h" 28 #include "libmesh/numeric_vector.h" 29 #include "libmesh/sparsity_pattern.h" 59 size += sparsity_pattern[row].size();
62 _row_start.reserve(n_rows + 1);
68 std::vector<numeric_index_type>::iterator position = _csr.begin();
70 _row_start.push_back (position);
75 for (
const auto & col : sparsity_pattern[row])
82 _row_start.push_back (position);
93 libmesh_assert_equal_to (n_rows, this->m());
99 auto rs = _row_start[i];
103 Q_SetLen (&_QMat, i+1, length);
117 libmesh_assert_equal_to (this->pos(i,j), l);
118 Q_SetEntry (&_QMat, i+1, l, j+1, 0.);
128 template <
typename T>
138 libmesh_assert_equal_to (m_in, m_l);
139 libmesh_assert_equal_to (n_in, n_l);
140 libmesh_assert_equal_to (m_in, n_in);
141 libmesh_assert_greater (nnz, 0);
145 libmesh_not_implemented_msg(
"Laspack does not support distributed matrices");
148 libmesh_not_implemented_msg(
"Laspack does not support rectangular matrices");
151 libmesh_warning(
"Using inefficient LaspackMatrix allocation via this init() method");
158 _csr.resize (m_in * m_in);
159 _row_start.reserve(m_in + 1);
164 std::vector<numeric_index_type>::iterator position = _csr.begin();
166 _row_start.push_back (position);
178 _row_start.push_back (position);
182 Q_Constr(&_QMat, const_cast<char *>(
"Mat"), m_in, _LPFalse, Rowws, Normal, _LPTrue);
190 auto rs = _row_start[i];
194 Q_SetLen (&_QMat, i+1, length);
200 libmesh_assert_equal_to (this->pos(i,j), l);
201 Q_SetEntry (&_QMat, i+1, l, j+1, 0.);
208 template <
typename T>
232 libmesh_assert_equal_to (n_rows, n_cols);
233 libmesh_assert_equal_to (m_l, n_rows);
234 libmesh_assert_equal_to (n_l, n_cols);
239 const std::vector<numeric_index_type> & n_nz = this->_sp->get_n_nz();
240 const std::vector<numeric_index_type> & n_oz = this->_sp->get_n_oz();
244 libmesh_assert_equal_to (n_nz.size(), n_l);
245 libmesh_assert_equal_to (n_oz.size(), n_l);
250 Q_Constr(&_QMat, const_cast<char *>(
"Mat"), n_rows, _LPFalse, Rowws, Normal, _LPTrue);
254 libmesh_assert_equal_to (n_rows, this->m());
259 template <
typename T>
261 const std::vector<numeric_index_type> & rows,
262 const std::vector<numeric_index_type> & cols)
266 unsigned int n_rows = cast_int<unsigned int>(rows.size());
267 unsigned int n_cols = cast_int<unsigned int>(cols.size());
268 libmesh_assert_equal_to (dm.
m(), n_rows);
269 libmesh_assert_equal_to (dm.
n(), n_cols);
272 for (
unsigned int i=0; i<n_rows; i++)
273 for (
unsigned int j=0; j<n_cols; j++)
274 this->add(rows[i],cols[j],dm(i,j));
279 template <
typename T>
285 std::vector<Real> abs_col_sums(this->n());
291 auto r_start = _row_start[row];
296 libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
303 libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
305 abs_col_sums[j] += std::abs(Q_GetVal (&_QMat, row+1, l));
309 return *(std::max_element(abs_col_sums.begin(), abs_col_sums.end()));
314 template <
typename T>
319 Real max_row_sum = 0;
328 libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
333 libmesh_assert_equal_to ((*(_row_start[row] + l)+1), Q_GetPos
336 row_sum += std::abs(Q_GetVal (&_QMat, row+1, l));
339 max_row_sum = std::max(max_row_sum, row_sum);
347 template <
typename T>
352 dest.
init(n_rows, n_rows);
355 dest.
set(i, Q_GetEl (const_cast<QMatrix*>(&_QMat), i+1, i+1));
360 template <
typename T>
367 libmesh_assert_equal_to(N, this->m());
369 std::vector<numeric_index_type> col_sizes(N);
373 auto r_start = _row_start[row];
381 target.
_csr.resize(_csr.size());
383 std::vector<std::vector<numeric_index_type>::iterator> writeable_row_start(N+1);
384 writeable_row_start[0] = target.
_csr.begin();
388 writeable_row_start[col+1] = writeable_row_start[col] + col_sizes[col];
392 std::vector<numeric_index_type> col_offsets;
394 col_offsets.swap(col_sizes);
395 std::fill(col_offsets.begin(), col_offsets.end(), 0);
397 Q_Constr(&target.
_QMat, const_cast<char *>(
"Mat"), N, _LPFalse, Rowws, Normal, _LPTrue);
403 auto r_start = _row_start[row];
407 *(writeable_row_start[col]+col_offsets[col]) = row;
414 auto len = col_offsets[col];
417 Q_SetLen(&target.
_QMat, col+1, len);
423 const auto value = Q_GetEl(const_cast<QMatrix*>(&_QMat), j+1, col+1);
424 Q_SetEntry (&target.
_QMat, col+1, l, j+1,
value);
431 template <
typename T>
433 std::vector<numeric_index_type> & indices,
434 std::vector<T> & values)
const 442 libmesh_assert_equal_to (len, Q_GetLen(&_QMat, i+1));
444 auto r_start = _row_start[i];
451 libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, l));
453 indices.push_back(j);
455 values.push_back(Q_GetVal (&_QMat, i+1, l));
461 template <
typename T>
470 template <
typename T>
478 template <
typename T>
494 template <
typename T>
501 auto r_start = _row_start[row];
506 libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
513 libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
515 Q_SetEntry (&_QMat, row+1, l, j+1, 0.);
524 template <
typename T>
528 auto mat_copy = std::make_unique<LaspackMatrix<T>>(this->comm());
530 mat_copy->attach_dof_map(*this->_dof_map);
532 mat_copy->init(this->m(), this->n(), this->local_m(),
533 this->local_n(), this->local_n());
542 template <
typename T>
547 auto mat_copy = this->zero_clone();
548 mat_copy->add(1., *
this);
553 template <
typename T>
563 template <
typename T>
573 template <
typename T>
581 template <
typename T>
589 template <
typename T>
597 template <
typename T>
605 template <
typename T>
611 libmesh_assert_less (i, this->m());
612 libmesh_assert_less (j, this->n());
617 libmesh_assert_equal_to (*(_row_start[i]+position), j);
618 libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, i+1, position));
620 Q_SetEntry (&_QMat, i+1, position, j+1,
value);
625 template <
typename T>
631 libmesh_assert_less (i, this->m());
632 libmesh_assert_less (j, this->n());
637 libmesh_assert_equal_to (*(_row_start[i]+position), j);
639 Q_AddVal (&_QMat, i+1, position,
value);
644 template <
typename T>
646 const std::vector<numeric_index_type> & dof_indices)
648 this->add_matrix (dm, dof_indices, dof_indices);
653 template <
typename T>
657 libmesh_assert_equal_to (this->m(), X_in.
m());
658 libmesh_assert_equal_to (this->n(), X_in.
n());
661 cast_ptr<const LaspackMatrix<T> *> (&X_in);
663 _LPNumber a =
static_cast<_LPNumber
> (a_in);
673 auto r_start = _row_start[row];
678 libmesh_assert_equal_to (len, Q_GetLen(&_QMat, row+1));
680 libmesh_assert_equal_to (len, Q_GetLen(&(X->_QMat), row+1));
688 libmesh_assert_equal_to ((j+1), Q_GetPos (&_QMat, row+1, l));
690 const _LPNumber
value = a * Q_GetEl(const_cast<QMatrix*>(&(X->_QMat)), row+1, j+1);
691 Q_AddVal (&_QMat, row+1, l,
value);
699 template <
typename T>
704 libmesh_assert_less (i, this->m());
705 libmesh_assert_less (j, this->n());
707 return Q_GetEl (const_cast<QMatrix*>(&_QMat), i+1, j+1);
712 template <
typename T>
716 libmesh_assert_less (i, this->m());
717 libmesh_assert_less (j, this->n());
718 libmesh_assert_less (i+1, _row_start.size());
722 auto p = std::equal_range (_row_start[i], _row_start[i+1], j);
736 template <
typename T>
741 this->_closed =
true;
745 *_QMat.DiagElAlloc = _LPFalse;
746 *_QMat.ElSorted = _LPFalse;
747 if (*_QMat.ILUExists)
749 *_QMat.ILUExists = _LPFalse;
763 #endif // #ifdef LIBMESH_HAVE_LASPACK virtual void init(const numeric_index_type m, const numeric_index_type n, const numeric_index_type m_l, const numeric_index_type n_l, const numeric_index_type nnz=30, const numeric_index_type noz=10, const numeric_index_type blocksize=1) override
Initialize SparseMatrix with the specified sizes.
virtual std::unique_ptr< SparseMatrix< T > > zero_clone() const override
virtual void get_diagonal(NumericVector< T > &dest) const override
Copies the diagonal part of the matrix into dest.
virtual numeric_index_type col_start() const override
std::vector< numeric_index_type > _csr
The compressed row indices.
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
QMatrix _QMat
The Laspack sparse matrix pointer.
The libMesh namespace provides an interface to certain functionality in the library.
Real distance(const Point &p)
virtual void init(const numeric_index_type n, const numeric_index_type n_local, const bool fast=false, const ParallelType ptype=AUTOMATIC)=0
Change the dimension of the vector to n.
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value) override
Add value to the element (i,j).
virtual std::unique_ptr< SparseMatrix< T > > clone() const override
numeric_index_type pos(const numeric_index_type i, const numeric_index_type j) const
virtual void get_transpose(SparseMatrix< T > &dest) const override
Copies the transpose of the matrix into dest, which may be *this.
virtual void close() override
Calls the SparseMatrix's internal assembly routines, ensuring that the values are consistent across p...
dof_id_type numeric_index_type
bool _is_initialized
Flag that tells if init() has been called.
virtual Real l1_norm() const override
bool _is_initialized
Flag indicating whether or not the matrix has been initialized.
virtual numeric_index_type m() const =0
virtual T operator()(const numeric_index_type i, const numeric_index_type j) const override
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
virtual Real linfty_norm() const override
virtual numeric_index_type row_start() const override
virtual void zero() override
Set all entries to 0.
virtual void clear() override
Restores the SparseMatrix<T> to a pristine state.
std::vector< std::vector< numeric_index_type >::const_iterator > _row_start
The start of each row in the compressed row index data structure.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual numeric_index_type m() const override
virtual numeric_index_type n() const override
IntRange< T > make_range(T beg, T end)
The 2-parameter make_range() helper function returns an IntRange<T> when both input parameters are of...
bool initialized()
Checks that library initialization has been done.
virtual void update_sparsity_pattern(const SparsityPattern::Graph &) override
Updates the matrix sparsity pattern.
virtual numeric_index_type row_stop() const override
virtual void set(const numeric_index_type i, const T value)=0
Sets v(i) = value.
The LaspackMatrix class wraps a QMatrix object from the Laspack library.
virtual void set(const numeric_index_type i, const numeric_index_type j, const T value) override
Set the element (i,j) to value.
Defines a dense matrix for use in Finite Element-type computations.
virtual void get_row(numeric_index_type i, std::vector< numeric_index_type > &indices, std::vector< T > &values) const override
Get a row from the matrix.
virtual numeric_index_type n() const =0
virtual numeric_index_type col_stop() const override
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols) override
Add the full matrix dm to the SparseMatrix.
LaspackMatrix(const Parallel::Communicator &comm)
Constructor; initializes the matrix to be empty, without any structure, i.e.
ParallelType
Defines an enum for parallel data structure types.