20 #include "libmesh/libmesh_config.h" 22 #ifdef LIBMESH_TRILINOS_HAVE_EPETRA 25 #include "libmesh/trilinos_epetra_matrix.h" 26 #include "libmesh/trilinos_epetra_vector.h" 27 #include "libmesh/dof_map.h" 28 #include "libmesh/dense_matrix.h" 29 #include "libmesh/parallel.h" 30 #include "libmesh/sparsity_pattern.h" 31 #include "libmesh/int_range.h" 51 (sparsity_pattern.size());
55 const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(this->processor_id());
61 libmesh_assert_equal_to (n, m);
62 libmesh_assert_equal_to (n_l, m_l);
68 this->comm().sum (summed_m_l);
69 this->comm().sum (summed_n_l);
71 libmesh_assert_equal_to (m, summed_m_l);
72 libmesh_assert_equal_to (n, summed_n_l);
77 _map =
new Epetra_Map (static_cast<int>(m),
80 Epetra_MpiComm (this->comm().
get()));
82 libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->NumGlobalPoints()), m);
83 libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->MaxAllGID()+1), m);
85 const std::vector<numeric_index_type> & n_nz = this->_sp->get_n_nz();
86 const std::vector<numeric_index_type> & n_oz = this->_sp->get_n_oz();
89 libmesh_assert_equal_to (n_nz.size(), n_l);
90 libmesh_assert_equal_to (n_oz.size(), n_l);
93 std::vector<int> n_nz_tot; n_nz_tot.reserve(n_nz.size());
96 n_nz_tot.push_back(std::min(n_nz[i] + n_oz[i], n));
101 _graph =
new Epetra_CrsGraph(Copy, *_map, n_nz_tot.data());
106 _graph->InsertGlobalIndices(_graph->GRID(i),
107 cast_int<numeric_index_type>(sparsity_pattern[i].size()),
108 const_cast<int *>(reinterpret_cast<const int *>(sparsity_pattern[i].data())));
110 _graph->FillComplete();
120 template <
typename T>
129 if ((m==0) || (n==0))
146 libmesh_assert_equal_to (n, m);
147 libmesh_assert_equal_to (n_l, m_l);
153 this->comm().sum (summed_m_l);
154 this->comm().sum (summed_n_l);
156 libmesh_assert_equal_to (m, summed_m_l);
157 libmesh_assert_equal_to (n, summed_n_l);
162 _map =
new Epetra_Map (static_cast<int>(m),
165 Epetra_MpiComm (this->comm().
get()));
167 libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->NumGlobalPoints()), m);
168 libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->MaxAllGID()+1), m);
170 _mat =
new Epetra_FECrsMatrix (Copy, *_map, nnz + noz);
176 template <
typename T>
190 _mat =
new Epetra_FECrsMatrix (Copy, *_graph);
195 template <
typename T>
205 template <
typename T>
210 libmesh_not_implemented();
213 auto mat_copy = std::make_unique<EpetraMatrix<T>>(this->comm());
219 return std::unique_ptr<SparseMatrix<T>>(mat_copy.release());
224 template <
typename T>
229 auto mat_copy = this->zero_clone();
230 mat_copy->add(1., *
this);
234 return std::unique_ptr<SparseMatrix<T>>(mat_copy.release());
239 template <
typename T>
252 template <
typename T>
259 return static_cast<Real>(_mat->NormOne());
264 template <
typename T>
272 return static_cast<Real>(_mat->NormInf());
277 template <
typename T>
279 const std::vector<numeric_index_type> & rows,
280 const std::vector<numeric_index_type> & cols)
287 libmesh_assert_equal_to (rows.size(), m);
288 libmesh_assert_equal_to (cols.size(), n);
300 template <
typename T>
307 _mat->ExtractDiagonalCopy(*(epetra_dest->
vec()));
312 template <
typename T>
322 if (&epetra_dest !=
this)
323 libmesh_not_implemented();
331 template <
typename T>
333 std::vector<numeric_index_type> & indices,
334 std::vector<T> & values)
const 339 libmesh_assert_greater_equal (i, this->row_start());
340 libmesh_assert_less (i, this->row_stop());
346 _mat->ExtractMyRowView (i-this->row_start(),
351 indices.resize(row_length);
352 values.resize(row_length);
356 indices[i] = row_indices[i];
357 values[i] = row_values[i];
363 template <
typename T>
366 _destroy_mat_on_exit(true),
367 _use_transpose(false)
373 template <
typename T>
377 _destroy_mat_on_exit(false),
378 _use_transpose(false)
387 template <
typename T>
395 template <
typename T>
400 _mat->GlobalAssemble();
405 template <
typename T>
415 template <
typename T>
425 template <
typename T>
436 template <
typename T>
447 template <
typename T>
458 template <
typename T>
469 template <
typename T>
477 epetra_i =
static_cast<int>(i),
478 epetra_j = static_cast<int>(j);
480 T epetra_value =
value;
483 _mat->ReplaceGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
485 _mat->InsertGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
490 template <
typename T>
498 epetra_i =
static_cast<int>(i),
499 epetra_j = static_cast<int>(j);
501 T epetra_value =
value;
503 _mat->SumIntoGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
508 template <
typename T>
510 const std::vector<numeric_index_type> & dof_indices)
512 this->add_matrix (dm, dof_indices, dof_indices);
517 template <
typename T>
520 #ifdef LIBMESH_TRILINOS_HAVE_EPETRAEXT 525 libmesh_assert_equal_to (this->m(), X_in.
m());
526 libmesh_assert_equal_to (this->n(), X_in.
n());
529 cast_ptr<const EpetraMatrix<T> *> (&X_in);
531 EpetraExt::MatrixMatrix::Add (*X->_mat,
false, a_in, *_mat, 1.);
533 libmesh_error_msg(
"ERROR: EpetraExt is required for EpetraMatrix::add()!");
540 template <
typename T>
547 libmesh_assert_greater_equal (i, this->row_start());
548 libmesh_assert_less (i, this->row_stop());
555 _mat->ExtractMyRowView (i-this->row_start(),
562 int * index = std::lower_bound (row_indices, row_indices+row_length, j);
564 libmesh_assert_less (*index, row_length);
565 libmesh_assert_equal_to (static_cast<numeric_index_type>(row_indices[*index]), j);
569 return values[*index];
575 template <
typename T>
581 return this->_mat->Filled();
585 template <
typename T>
588 std::swap(_mat, m.
_mat);
596 template <
typename T>
614 #endif // LIBMESH_TRILINOS_HAVE_EPETRA This class provides a nice interface to the Trilinos Epetra_Vector object.
virtual std::unique_ptr< SparseMatrix< T > > clone() const override
virtual numeric_index_type row_stop() const override
virtual void get_transpose(SparseMatrix< T > &dest) const override
Copies the transpose of the matrix into dest, which may be *this.
virtual void clear() noexcept override
clear() is called from the destructor, so it should not throw.
bool _destroy_mat_on_exit
This boolean value should only be set to false for the constructor which takes an Epetra_FECrsMatrix ...
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
bool _use_transpose
Epetra has no GetUseTranspose so we need to keep track of whether we're transposed manually...
virtual std::unique_ptr< SparseMatrix< T > > zero_clone() const override
The libMesh namespace provides an interface to certain functionality in the 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 numeric_index_type col_stop() const override
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 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.
virtual numeric_index_type n() const override
virtual void set(const numeric_index_type i, const numeric_index_type j, const T value) override
Set the element (i,j) to value.
virtual numeric_index_type col_start() const override
dof_id_type numeric_index_type
bool _is_initialized
Flag that tells if init() has been called.
bool _is_initialized
Flag indicating whether or not the matrix has been initialized.
virtual numeric_index_type m() const override
virtual numeric_index_type m() const =0
void init(triangulateio &t)
Initializes the fields of t to nullptr/0 as necessary.
std::vector< T > & get_values()
virtual numeric_index_type row_start() const override
virtual void close() override
Calls the SparseMatrix's internal assembly routines, ensuring that the values are consistent across p...
virtual void print_personal(std::ostream &os=libMesh::out) const override
Print the contents of the matrix to the screen in a package-personalized style, if available...
EpetraMatrix(const Parallel::Communicator &comm)
Constructor; initializes the matrix to be empty, without any structure, i.e.
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual Real l1_norm() const override
virtual Real linfty_norm() const override
virtual void get_diagonal(NumericVector< T > &dest) const override
Copies the diagonal part of the matrix into dest.
virtual bool closed() const override
virtual void zero() override
Set all entries to 0.
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.
int * numeric_trilinos_cast(const numeric_index_type *p)
Epetra_FECrsMatrix * _mat
Actual Epetra datatype to hold matrix entries.
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value) override
Add value to the element (i,j).
Defines a dense matrix for use in Finite Element-type computations.
void swap(EpetraMatrix< T > &)
Swaps the internal data pointers, no actual values are swapped.
This class provides a nice interface to the Epetra data structures for parallel, sparse matrices...
virtual T operator()(const numeric_index_type i, const numeric_index_type j) const override
virtual void update_sparsity_pattern(const SparsityPattern::Graph &) override
Updates the matrix sparsity pattern.
auto index_range(const T &sizable)
Helper function that returns an IntRange<std::size_t> representing all the indices of the passed-in v...
virtual numeric_index_type n() const =0
ParallelType
Defines an enum for parallel data structure types.