21 #include "libmesh/libmesh_config.h" 23 #ifdef LIBMESH_HAVE_EIGEN 25 #include "libmesh/eigen_sparse_vector.h" 26 #include "libmesh/eigen_sparse_matrix.h" 27 #include "libmesh/dense_matrix.h" 28 #include "libmesh/dof_map.h" 29 #include "libmesh/sparsity_pattern.h" 51 libmesh_assert_equal_to (m_in, m_l);
52 libmesh_assert_equal_to (n_in, n_l);
53 libmesh_assert_equal_to (m_in, n_in);
54 libmesh_assert_greater (nnz, 0);
56 _mat.resize(m_in, n_in);
57 _mat.reserve(Eigen::Matrix<numeric_index_type, Eigen::Dynamic, 1>::Constant(m_in,nnz));
89 libmesh_assert_equal_to (n_rows, n_cols);
90 libmesh_assert_equal_to (m_l, n_rows);
91 libmesh_assert_equal_to (n_l, n_cols);
93 const std::vector<numeric_index_type> & n_nz = this->_sp->get_n_nz();
98 const std::vector<numeric_index_type> & n_oz = this->_sp->get_n_oz();
102 libmesh_assert_equal_to (n_nz.size(), n_l);
103 libmesh_assert_equal_to (n_oz.size(), n_l);
111 _mat.resize(n_rows,n_cols);
116 libmesh_assert_equal_to (n_rows, this->m());
117 libmesh_assert_equal_to (n_cols, this->n());
122 template <
typename T>
124 const std::vector<numeric_index_type> & rows,
125 const std::vector<numeric_index_type> & cols)
129 unsigned int n_rows = cast_int<unsigned int>(rows.size());
130 unsigned int n_cols = cast_int<unsigned int>(cols.size());
131 libmesh_assert_equal_to (dm.
m(), n_rows);
132 libmesh_assert_equal_to (dm.
n(), n_cols);
135 for (
unsigned int i=0; i<n_rows; i++)
136 for (
unsigned int j=0; j<n_cols; j++)
137 this->add(rows[i],cols[j],dm(i,j));
142 template <
typename T>
147 dest.
_vec = _mat.diagonal();
152 template <
typename T>
157 dest.
_mat = _mat.transpose();
162 template <
typename T>
171 template <
typename T>
182 template <
typename T>
191 const std::vector<numeric_index_type> & n_nz = this->_sp->get_n_nz();
198 template <
typename T>
203 auto ret = std::make_unique<EigenSparseMatrix<T>>(*this);
208 return std::unique_ptr<SparseMatrix<T>>(ret.release());
213 template <
typename T>
216 return std::make_unique<EigenSparseMatrix<T>>(*this);
221 template <
typename T>
226 return cast_int<numeric_index_type>(_mat.rows());
231 template <
typename T>
236 return cast_int<numeric_index_type>(_mat.cols());
241 template <
typename T>
249 template <
typename T>
257 template <
typename T>
265 template <
typename T>
273 template <
typename T>
279 libmesh_assert_less (i, this->m());
280 libmesh_assert_less (j, this->n());
282 _mat.coeffRef(i,j) =
value;
287 template <
typename T>
293 libmesh_assert_less (i, this->m());
294 libmesh_assert_less (j, this->n());
296 _mat.coeffRef(i,j) +=
value;
301 template <
typename T>
303 const std::vector<numeric_index_type> & dof_indices)
305 this->add_matrix (dm, dof_indices, dof_indices);
310 template <
typename T>
314 libmesh_assert_equal_to (this->m(), X_in.
m());
315 libmesh_assert_equal_to (this->n(), X_in.
n());
318 cast_ref<const EigenSparseMatrix<T> &> (X_in);
326 template <
typename T>
331 libmesh_assert_less (i, this->m());
332 libmesh_assert_less (j, this->n());
334 return _mat.coeff(i,j);
339 template <
typename T>
346 std::vector<Real> abs_col_sums(this->n());
352 EigenSM::InnerIterator it(_mat, row);
354 abs_col_sums[it.col()] +=
std::abs(it.value());
357 return *(std::max_element(abs_col_sums.begin(), abs_col_sums.end()));
362 template <
typename T>
365 Real max_abs_row_sum = 0.;
371 Real current_abs_row_sum = 0.;
372 EigenSM::InnerIterator it(_mat, row);
374 current_abs_row_sum +=
std::abs(it.value());
376 max_abs_row_sum = std::max(max_abs_row_sum, current_abs_row_sum);
379 return max_abs_row_sum;
384 template <
typename T>
386 std::vector<numeric_index_type> & indices,
387 std::vector<T> & values)
const 393 static_assert(EigenSM::IsRowMajor);
395 for (EigenSM::InnerIterator it(_mat, i); it; ++it)
397 indices.push_back(it.col());
398 values.push_back(it.value());
411 #endif // #ifdef LIBMESH_HAVE_EIGEN virtual void set(const numeric_index_type i, const numeric_index_type j, const T value) override
Set the element (i,j) to value.
The EigenSparseMatrix class wraps a sparse matrix object from the Eigen library.
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 void add(const numeric_index_type i, const numeric_index_type j, const T value) override
Add value to the element (i,j).
virtual numeric_index_type m() const override
virtual numeric_index_type row_stop() const override
DataType _vec
Actual Eigen::SparseVector<> we are wrapping.
virtual std::unique_ptr< SparseMatrix< T > > clone() const override
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
The libMesh namespace provides an interface to certain functionality in the library.
DataType _mat
Actual Eigen::SparseMatrix<> we are wrapping.
virtual void get_transpose(SparseMatrix< T > &dest) const override
Copies the transpose of the matrix into dest, which may be *this.
ADRealEigenVector< T, D, asd > abs(const ADRealEigenVector< T, D, asd > &)
virtual void zero() override
Set all entries to 0.
virtual numeric_index_type n() const override
virtual void get_diagonal(NumericVector< T > &dest) const override
Copies the diagonal part of the matrix into dest.
dof_id_type numeric_index_type
bool _is_initialized
Flag that tells if init() has been called.
EigenSparseMatrix(const Parallel::Communicator &comm)
Constructor; initializes the matrix to be empty, without any structure, i.e.
virtual numeric_index_type m() const =0
virtual numeric_index_type col_stop() const override
virtual T operator()(const numeric_index_type i, const numeric_index_type j) const override
virtual numeric_index_type row_start() const override
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual Real linfty_norm() const override
virtual numeric_index_type col_start() const override
virtual void clear() override
Restores the SparseMatrix<T> to a pristine state.
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.
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 std::unique_ptr< SparseMatrix< T > > zero_clone() const override
virtual Real l1_norm() const override
Defines a dense matrix for use in Finite Element-type computations.
This class provides a nice interface to the Eigen C++-based data structures for serial vectors...
virtual numeric_index_type n() const =0
ParallelType
Defines an enum for parallel data structure types.
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.