libMesh
trilinos_epetra_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2024 Benjamin S. Kirk, John W. Peterson, Roy H. Stogner
3 
4 // This library is free software; you can redistribute it and/or
5 // modify it under the terms of the GNU Lesser General Public
6 // License as published by the Free Software Foundation; either
7 // version 2.1 of the License, or (at your option) any later version.
8 
9 // This library is distributed in the hope that it will be useful,
10 // but WITHOUT ANY WARRANTY; without even the implied warranty of
11 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
12 // Lesser General Public License for more details.
13 
14 // You should have received a copy of the GNU Lesser General Public
15 // License along with this library; if not, write to the Free Software
16 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307 USA
17 
18 
19 // C++ includes
20 #include "libmesh/libmesh_config.h"
21 
22 #ifdef LIBMESH_TRILINOS_HAVE_EPETRA
23 
24 // Local includes
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"
32 
33 namespace libMesh
34 {
35 
36 
37 
38 //-----------------------------------------------------------------------
39 //EpetraMatrix members
40 
41 template <typename T>
43 {
44  // clear data, start over
45  this->clear ();
46 
47  // big trouble if this fails!
48  libmesh_assert(this->_dof_map);
49 
50  const numeric_index_type n_rows = cast_int<numeric_index_type>
51  (sparsity_pattern.size());
52 
53  const numeric_index_type m = this->_dof_map->n_dofs();
54  const numeric_index_type n = m;
55  const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(this->processor_id());
56  const numeric_index_type m_l = n_l;
57 
58  // error checking
59 #ifndef NDEBUG
60  {
61  libmesh_assert_equal_to (n, m);
62  libmesh_assert_equal_to (n_l, m_l);
63 
65  summed_m_l = m_l,
66  summed_n_l = n_l;
67 
68  this->comm().sum (summed_m_l);
69  this->comm().sum (summed_n_l);
70 
71  libmesh_assert_equal_to (m, summed_m_l);
72  libmesh_assert_equal_to (n, summed_n_l);
73  }
74 #endif
75 
76  // build a map defining the data distribution
77  _map = new Epetra_Map (static_cast<int>(m),
78  m_l,
79  0,
80  Epetra_MpiComm (this->comm().get()));
81 
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);
84 
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();
87 
88  // Make sure the sparsity pattern isn't empty
89  libmesh_assert_equal_to (n_nz.size(), n_l);
90  libmesh_assert_equal_to (n_oz.size(), n_l);
91 
92  // Epetra wants the total number of nonzeros, both local and remote.
93  std::vector<int> n_nz_tot; n_nz_tot.reserve(n_nz.size());
94 
95  for (auto i : index_range(n_nz))
96  n_nz_tot.push_back(std::min(n_nz[i] + n_oz[i], n));
97 
98  if (m==0)
99  return;
100 
101  _graph = new Epetra_CrsGraph(Copy, *_map, n_nz_tot.data());
102 
103  // Tell the matrix about its structure. Initialize it
104  // to zero.
105  for (numeric_index_type i=0; i<n_rows; i++)
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())));
109 
110  _graph->FillComplete();
111 
112  //Initialize the matrix
113  libmesh_assert (!this->initialized());
114  this->init ();
115  libmesh_assert (this->initialized());
116 }
117 
118 
119 
120 template <typename T>
122  const numeric_index_type n,
123  const numeric_index_type m_l,
124  const numeric_index_type libmesh_dbg_var(n_l),
125  const numeric_index_type nnz,
126  const numeric_index_type noz,
127  const numeric_index_type /* blocksize */)
128 {
129  if ((m==0) || (n==0))
130  return;
131 
132  {
133  // Clear initialized matrices
134  if (this->initialized())
135  this->clear();
136 
137  libmesh_assert (!this->_mat);
138  libmesh_assert (!this->_map);
139 
140  this->_is_initialized = true;
141  }
142 
143  // error checking
144 #ifndef NDEBUG
145  {
146  libmesh_assert_equal_to (n, m);
147  libmesh_assert_equal_to (n_l, m_l);
148 
150  summed_m_l = m_l,
151  summed_n_l = n_l;
152 
153  this->comm().sum (summed_m_l);
154  this->comm().sum (summed_n_l);
155 
156  libmesh_assert_equal_to (m, summed_m_l);
157  libmesh_assert_equal_to (n, summed_n_l);
158  }
159 #endif
160 
161  // build a map defining the data distribution
162  _map = new Epetra_Map (static_cast<int>(m),
163  m_l,
164  0,
165  Epetra_MpiComm (this->comm().get()));
166 
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);
169 
170  _mat = new Epetra_FECrsMatrix (Copy, *_map, nnz + noz);
171 }
172 
173 
174 
175 
176 template <typename T>
178 {
179  libmesh_assert(this->_dof_map);
180 
181  {
182  // Clear initialized matrices
183  if (this->initialized())
184  this->clear();
185 
186  this->_is_initialized = true;
187  }
188 
189 
190  _mat = new Epetra_FECrsMatrix (Copy, *_graph);
191 }
192 
193 
194 
195 template <typename T>
197 {
198  libmesh_assert (this->initialized());
199 
200  _mat->Scale(0.0);
201 }
202 
203 
204 
205 template <typename T>
206 std::unique_ptr<SparseMatrix<T>> EpetraMatrix<T>::zero_clone () const
207 {
208  // This function is marked as "not implemented" since it hasn't been
209  // tested, the code below might serve as a possible implementation.
210  libmesh_not_implemented();
211 
212  // Make empty copy with matching comm, initialize, and return.
213  auto mat_copy = std::make_unique<EpetraMatrix<T>>(this->comm());
214  mat_copy->init();
215  mat_copy->zero();
216 
217  // Work around an issue on older compilers. We are able to simply
218  // "return mat_copy;" on newer compilers
219  return std::unique_ptr<SparseMatrix<T>>(mat_copy.release());
220 }
221 
222 
223 
224 template <typename T>
225 std::unique_ptr<SparseMatrix<T>> EpetraMatrix<T>::clone () const
226 {
227  // We don't currently have a faster implementation than making a
228  // zero clone and then filling in the values.
229  auto mat_copy = this->zero_clone();
230  mat_copy->add(1., *this);
231 
232  // Work around an issue on older compilers. We are able to simply
233  // "return mat_copy;" on newer compilers
234  return std::unique_ptr<SparseMatrix<T>>(mat_copy.release());
235 }
236 
237 
238 
239 template <typename T>
240 void EpetraMatrix<T>::clear () noexcept
241 {
242  // FIXME: clear() doesn't actually free the memory managed by this
243  // class, so it probably leaks memory.
244  // delete _mat;
245  // delete _map;
246 
247  this->_is_initialized = false;
248 }
249 
250 
251 
252 template <typename T>
254 {
255  libmesh_assert (this->initialized());
256 
257  libmesh_assert(_mat);
258 
259  return static_cast<Real>(_mat->NormOne());
260 }
261 
262 
263 
264 template <typename T>
266 {
267  libmesh_assert (this->initialized());
268 
269 
270  libmesh_assert(_mat);
271 
272  return static_cast<Real>(_mat->NormInf());
273 }
274 
275 
276 
277 template <typename T>
279  const std::vector<numeric_index_type> & rows,
280  const std::vector<numeric_index_type> & cols)
281 {
282  libmesh_assert (this->initialized());
283 
284  const numeric_index_type m = dm.m();
285  const numeric_index_type n = dm.n();
286 
287  libmesh_assert_equal_to (rows.size(), m);
288  libmesh_assert_equal_to (cols.size(), n);
289 
290  _mat->SumIntoGlobalValues(m, numeric_trilinos_cast(rows.data()),
291  n, numeric_trilinos_cast(cols.data()),
292  dm.get_values().data());
293 }
294 
295 
296 
297 
298 
299 
300 template <typename T>
302 {
303  // Convert vector to EpetraVector.
304  EpetraVector<T> * epetra_dest = cast_ptr<EpetraVector<T> *>(&dest);
305 
306  // Call Epetra function.
307  _mat->ExtractDiagonalCopy(*(epetra_dest->vec()));
308 }
309 
310 
311 
312 template <typename T>
314 {
315  // Make sure the SparseMatrix passed in is really a EpetraMatrix
316  EpetraMatrix<T> & epetra_dest = cast_ref<EpetraMatrix<T> &>(dest);
317 
318  // We currently only support calling get_transpose() with ourself
319  // as the destination. Previously, this called the default copy
320  // constructor which was not safe because this class manually
321  // manages memory.
322  if (&epetra_dest != this)
323  libmesh_not_implemented();
324 
325  epetra_dest._use_transpose = !epetra_dest._use_transpose;
326  epetra_dest._mat->SetUseTranspose(epetra_dest._use_transpose);
327 }
328 
329 
330 
331 template <typename T>
333  std::vector<numeric_index_type> & indices,
334  std::vector<T> & values) const
335 {
336  libmesh_assert (this->initialized());
337  libmesh_assert(this->_mat);
338  libmesh_assert (this->_mat->MyGlobalRow(static_cast<int>(i)));
339  libmesh_assert_greater_equal (i, this->row_start());
340  libmesh_assert_less (i, this->row_stop());
341 
342  int row_length;
343  int * row_indices;
344  double * row_values;
345 
346  _mat->ExtractMyRowView (i-this->row_start(),
347  row_length,
348  row_values,
349  row_indices);
350 
351  indices.resize(row_length);
352  values.resize(row_length);
353 
354  for (auto i : make_range(row_length))
355  {
356  indices[i] = row_indices[i];
357  values[i] = row_values[i];
358  }
359 }
360 
361 
362 
363 template <typename T>
365  SparseMatrix<T>(comm),
366  _destroy_mat_on_exit(true),
367  _use_transpose(false)
368 {}
369 
370 
371 
372 
373 template <typename T>
374 EpetraMatrix<T>::EpetraMatrix(Epetra_FECrsMatrix * m,
375  const Parallel::Communicator & comm) :
376  SparseMatrix<T>(comm),
377  _destroy_mat_on_exit(false),
378  _use_transpose(false) // dumb guess is the best we can do...
379 {
380  this->_mat = m;
381  this->_is_initialized = true;
382 }
383 
384 
385 
386 
387 template <typename T>
389 {
390  this->clear();
391 }
392 
393 
394 
395 template <typename T>
397 {
398  libmesh_assert(_mat);
399 
400  _mat->GlobalAssemble();
401 }
402 
403 
404 
405 template <typename T>
407 {
408  libmesh_assert (this->initialized());
409 
410  return static_cast<numeric_index_type>(_mat->NumGlobalRows());
411 }
412 
413 
414 
415 template <typename T>
417 {
418  libmesh_assert (this->initialized());
419 
420  return static_cast<numeric_index_type>(_mat->NumGlobalCols());
421 }
422 
423 
424 
425 template <typename T>
427 {
428  libmesh_assert (this->initialized());
429  libmesh_assert(_map);
430 
431  return static_cast<numeric_index_type>(_map->MinMyGID());
432 }
433 
434 
435 
436 template <typename T>
438 {
439  libmesh_assert (this->initialized());
440  libmesh_assert(_map);
441 
442  return static_cast<numeric_index_type>(_map->MaxMyGID())+1;
443 }
444 
445 
446 
447 template <typename T>
449 {
450  libmesh_assert (this->initialized());
451  libmesh_assert(_map);
452 
453  return static_cast<numeric_index_type>(_map->MinMyGID());
454 }
455 
456 
457 
458 template <typename T>
460 {
461  libmesh_assert (this->initialized());
462  libmesh_assert(_map);
463 
464  return static_cast<numeric_index_type>(_map->MaxMyGID())+1;
465 }
466 
467 
468 
469 template <typename T>
471  const numeric_index_type j,
472  const T value)
473 {
474  libmesh_assert (this->initialized());
475 
476  int
477  epetra_i = static_cast<int>(i),
478  epetra_j = static_cast<int>(j);
479 
480  T epetra_value = value;
481 
482  if (_mat->Filled())
483  _mat->ReplaceGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
484  else
485  _mat->InsertGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
486 }
487 
488 
489 
490 template <typename T>
492  const numeric_index_type j,
493  const T value)
494 {
495  libmesh_assert (this->initialized());
496 
497  int
498  epetra_i = static_cast<int>(i),
499  epetra_j = static_cast<int>(j);
500 
501  T epetra_value = value;
502 
503  _mat->SumIntoGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
504 }
505 
506 
507 
508 template <typename T>
510  const std::vector<numeric_index_type> & dof_indices)
511 {
512  this->add_matrix (dm, dof_indices, dof_indices);
513 }
514 
515 
516 
517 template <typename T>
518 void EpetraMatrix<T>::add (const T a_in, const SparseMatrix<T> & X_in)
519 {
520 #ifdef LIBMESH_TRILINOS_HAVE_EPETRAEXT
521  libmesh_assert (this->initialized());
522 
523  // sanity check. but this cannot avoid
524  // crash due to incompatible sparsity structure...
525  libmesh_assert_equal_to (this->m(), X_in.m());
526  libmesh_assert_equal_to (this->n(), X_in.n());
527 
528  const EpetraMatrix<T> * X =
529  cast_ptr<const EpetraMatrix<T> *> (&X_in);
530 
531  EpetraExt::MatrixMatrix::Add (*X->_mat, false, a_in, *_mat, 1.);
532 #else
533  libmesh_error_msg("ERROR: EpetraExt is required for EpetraMatrix::add()!");
534 #endif
535 }
536 
537 
538 
539 
540 template <typename T>
542  const numeric_index_type j) const
543 {
544  libmesh_assert (this->initialized());
545  libmesh_assert(this->_mat);
546  libmesh_assert (this->_mat->MyGlobalRow(static_cast<int>(i)));
547  libmesh_assert_greater_equal (i, this->row_start());
548  libmesh_assert_less (i, this->row_stop());
549 
550 
551  int row_length;
552  int * row_indices;
553  double * values;
554 
555  _mat->ExtractMyRowView (i-this->row_start(),
556  row_length,
557  values,
558  row_indices);
559 
560  //libMesh::out << "row_length=" << row_length << std::endl;
561 
562  int * index = std::lower_bound (row_indices, row_indices+row_length, j);
563 
564  libmesh_assert_less (*index, row_length);
565  libmesh_assert_equal_to (static_cast<numeric_index_type>(row_indices[*index]), j);
566 
567  //libMesh::out << "val=" << values[*index] << std::endl;
568 
569  return values[*index];
570 }
571 
572 
573 
574 
575 template <typename T>
577 {
578  libmesh_assert (this->initialized());
579  libmesh_assert(this->_mat);
580 
581  return this->_mat->Filled();
582 }
583 
584 
585 template <typename T>
587 {
588  std::swap(_mat, m._mat);
589  std::swap(_destroy_mat_on_exit, m._destroy_mat_on_exit);
590 }
591 
592 
593 
594 
595 
596 template <typename T>
597 void EpetraMatrix<T>::print_personal(std::ostream & os) const
598 {
599  libmesh_assert (this->initialized());
600  libmesh_assert(_mat);
601 
602  os << *_mat;
603 }
604 
605 
606 
607 //------------------------------------------------------------------
608 // Explicit instantiations
609 template class LIBMESH_EXPORT EpetraMatrix<Number>;
610 
611 } // namespace libMesh
612 
613 
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...
Definition: vector_fe_ex5.C:43
bool _use_transpose
Epetra has no GetUseTranspose so we need to keep track of whether we&#39;re transposed manually...
unsigned int m() const
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.
Generic sparse matrix.
Definition: vector_fe_ex5.C:45
virtual numeric_index_type col_start() const override
dof_id_type numeric_index_type
Definition: id_types.h:99
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:247
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.
libmesh_assert(ctx)
std::vector< T > & get_values()
Definition: dense_matrix.h:382
virtual numeric_index_type row_start() const override
virtual void close() override
Calls the SparseMatrix&#39;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.
static const bool value
Definition: xdr_io.C:54
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...
Definition: int_range.h:134
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:266
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).
unsigned int n() const
Defines a dense matrix for use in Finite Element-type computations.
Definition: dof_map.h:65
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...
Definition: int_range.h:111
virtual numeric_index_type n() const =0
ParallelType
Defines an enum for parallel data structure types.