libMesh
trilinos_epetra_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2017 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 
32 namespace libMesh
33 {
34 
35 
36 
37 //-----------------------------------------------------------------------
38 //EpetraMatrix members
39 
40 template <typename T>
41 void EpetraMatrix<T>::update_sparsity_pattern (const SparsityPattern::Graph & sparsity_pattern)
42 {
43  // clear data, start over
44  this->clear ();
45 
46  // big trouble if this fails!
47  libmesh_assert(this->_dof_map);
48 
49  const numeric_index_type n_rows = sparsity_pattern.size();
50 
51  const numeric_index_type m = this->_dof_map->n_dofs();
52  const numeric_index_type n = m;
53  const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(this->processor_id());
54  const numeric_index_type m_l = n_l;
55 
56  // error checking
57 #ifndef NDEBUG
58  {
59  libmesh_assert_equal_to (n, m);
60  libmesh_assert_equal_to (n_l, m_l);
61 
63  summed_m_l = m_l,
64  summed_n_l = n_l;
65 
66  this->comm().sum (summed_m_l);
67  this->comm().sum (summed_n_l);
68 
69  libmesh_assert_equal_to (m, summed_m_l);
70  libmesh_assert_equal_to (n, summed_n_l);
71  }
72 #endif
73 
74  // build a map defining the data distribution
75  _map = new Epetra_Map (static_cast<int>(m),
76  m_l,
77  0,
78  Epetra_MpiComm (this->comm().get()));
79 
80  libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->NumGlobalPoints()), m);
81  libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->MaxAllGID()+1), m);
82 
83  const std::vector<numeric_index_type> & n_nz = this->_dof_map->get_n_nz();
84  const std::vector<numeric_index_type> & n_oz = this->_dof_map->get_n_oz();
85 
86  // Make sure the sparsity pattern isn't empty
87  libmesh_assert_equal_to (n_nz.size(), n_l);
88  libmesh_assert_equal_to (n_oz.size(), n_l);
89 
90  // Epetra wants the total number of nonzeros, both local and remote.
91  std::vector<int> n_nz_tot; n_nz_tot.reserve(n_nz.size());
92 
93  for (numeric_index_type i=0; i<n_nz.size(); i++)
94  n_nz_tot.push_back(std::min(n_nz[i] + n_oz[i], n));
95 
96  if (m==0)
97  return;
98 
99  _graph = new Epetra_CrsGraph(Copy, *_map, &n_nz_tot[0]);
100 
101  // Tell the matrix about its structure. Initialize it
102  // to zero.
103  for (numeric_index_type i=0; i<n_rows; i++)
104  _graph->InsertGlobalIndices(_graph->GRID(i),
105  sparsity_pattern[i].size(),
106  const_cast<int *>((const int *)&sparsity_pattern[i][0]));
107 
108  _graph->FillComplete();
109 
110  //Initialize the matrix
111  libmesh_assert (!this->initialized());
112  this->init ();
113  libmesh_assert (this->initialized());
114 }
115 
116 
117 
118 template <typename T>
120  const numeric_index_type n,
121  const numeric_index_type m_l,
123  const numeric_index_type nnz,
124  const numeric_index_type noz,
125  const numeric_index_type /* blocksize */)
126 {
127  if ((m==0) || (n==0))
128  return;
129 
130  {
131  // Clear initialized matrices
132  if (this->initialized())
133  this->clear();
134 
135  libmesh_assert (!this->_mat);
136  libmesh_assert (!this->_map);
137 
138  this->_is_initialized = true;
139  }
140 
141  // error checking
142 #ifndef NDEBUG
143  {
144  libmesh_assert_equal_to (n, m);
145  libmesh_assert_equal_to (n_l, m_l);
146 
148  summed_m_l = m_l,
149  summed_n_l = n_l;
150 
151  this->comm().sum (summed_m_l);
152  this->comm().sum (summed_n_l);
153 
154  libmesh_assert_equal_to (m, summed_m_l);
155  libmesh_assert_equal_to (n, summed_n_l);
156  }
157 #endif
158 
159  // build a map defining the data distribution
160  _map = new Epetra_Map (static_cast<int>(m),
161  m_l,
162  0,
163  Epetra_MpiComm (this->comm().get()));
164 
165  libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->NumGlobalPoints()), m);
166  libmesh_assert_equal_to (static_cast<numeric_index_type>(_map->MaxAllGID()+1), m);
167 
168  _mat = new Epetra_FECrsMatrix (Copy, *_map, nnz + noz);
169 }
170 
171 
172 
173 
174 template <typename T>
175 void EpetraMatrix<T>::init ()
176 {
177  libmesh_assert(this->_dof_map);
178 
179  {
180  // Clear initialized matrices
181  if (this->initialized())
182  this->clear();
183 
184  this->_is_initialized = true;
185  }
186 
187 
188  _mat = new Epetra_FECrsMatrix (Copy, *_graph);
189 }
190 
191 
192 
193 template <typename T>
194 void EpetraMatrix<T>::zero ()
195 {
196  libmesh_assert (this->initialized());
197 
198  _mat->Scale(0.0);
199 }
200 
201 
202 
203 template <typename T>
204 void EpetraMatrix<T>::clear ()
205 {
206  // delete _mat;
207  // delete _map;
208 
209  this->_is_initialized = false;
210 
211  libmesh_assert (!this->initialized());
212 }
213 
214 
215 
216 template <typename T>
217 Real EpetraMatrix<T>::l1_norm () const
218 {
219  libmesh_assert (this->initialized());
220 
221  libmesh_assert(_mat);
222 
223  return static_cast<Real>(_mat->NormOne());
224 }
225 
226 
227 
228 template <typename T>
229 Real EpetraMatrix<T>::linfty_norm () const
230 {
231  libmesh_assert (this->initialized());
232 
233 
234  libmesh_assert(_mat);
235 
236  return static_cast<Real>(_mat->NormInf());
237 }
238 
239 
240 
241 template <typename T>
242 void EpetraMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
243  const std::vector<numeric_index_type> & rows,
244  const std::vector<numeric_index_type> & cols)
245 {
246  libmesh_assert (this->initialized());
247 
248  const numeric_index_type m = dm.m();
249  const numeric_index_type n = dm.n();
250 
251  libmesh_assert_equal_to (rows.size(), m);
252  libmesh_assert_equal_to (cols.size(), n);
253 
254  _mat->SumIntoGlobalValues(m, (int *)&rows[0], n, (int *)&cols[0], &dm.get_values()[0]);
255 }
256 
257 
258 
259 
260 
261 
262 template <typename T>
263 void EpetraMatrix<T>::get_diagonal (NumericVector<T> & dest) const
264 {
265  // Convert vector to EpetraVector.
266  EpetraVector<T> * epetra_dest = cast_ptr<EpetraVector<T> *>(&dest);
267 
268  // Call Epetra function.
269  _mat->ExtractDiagonalCopy(*(epetra_dest->vec()));
270 }
271 
272 
273 
274 template <typename T>
275 void EpetraMatrix<T>::get_transpose (SparseMatrix<T> & dest) const
276 {
277  // Make sure the SparseMatrix passed in is really a EpetraMatrix
278  EpetraMatrix<T> & epetra_dest = cast_ref<EpetraMatrix<T> &>(dest);
279 
280  if (&epetra_dest != this)
281  epetra_dest = *this;
282 
283  epetra_dest._use_transpose = !epetra_dest._use_transpose;
284  epetra_dest._mat->SetUseTranspose(epetra_dest._use_transpose);
285 }
286 
287 
288 
289 template <typename T>
290 EpetraMatrix<T>::EpetraMatrix(const Parallel::Communicator & comm) :
291  SparseMatrix<T>(comm),
292  _destroy_mat_on_exit(true),
293  _use_transpose(false)
294 {}
295 
296 
297 
298 
299 template <typename T>
300 EpetraMatrix<T>::EpetraMatrix(Epetra_FECrsMatrix * m,
301  const Parallel::Communicator & comm) :
302  SparseMatrix<T>(comm),
303  _destroy_mat_on_exit(false),
304  _use_transpose(false) // dumb guess is the best we can do...
305 {
306  this->_mat = m;
307  this->_is_initialized = true;
308 }
309 
310 
311 
312 
313 template <typename T>
314 EpetraMatrix<T>::~EpetraMatrix()
315 {
316  this->clear();
317 }
318 
319 
320 
321 template <typename T>
322 void EpetraMatrix<T>::close ()
323 {
324  libmesh_assert(_mat);
325 
326  _mat->GlobalAssemble();
327 }
328 
329 
330 
331 template <typename T>
332 numeric_index_type EpetraMatrix<T>::m () const
333 {
334  libmesh_assert (this->initialized());
335 
336  return static_cast<numeric_index_type>(_mat->NumGlobalRows());
337 }
338 
339 
340 
341 template <typename T>
342 numeric_index_type EpetraMatrix<T>::n () const
343 {
344  libmesh_assert (this->initialized());
345 
346  return static_cast<numeric_index_type>(_mat->NumGlobalCols());
347 }
348 
349 
350 
351 template <typename T>
352 numeric_index_type EpetraMatrix<T>::row_start () const
353 {
354  libmesh_assert (this->initialized());
355  libmesh_assert(_map);
356 
357  return static_cast<numeric_index_type>(_map->MinMyGID());
358 }
359 
360 
361 
362 template <typename T>
363 numeric_index_type EpetraMatrix<T>::row_stop () const
364 {
365  libmesh_assert (this->initialized());
366  libmesh_assert(_map);
367 
368  return static_cast<numeric_index_type>(_map->MaxMyGID())+1;
369 }
370 
371 
372 
373 template <typename T>
374 void EpetraMatrix<T>::set (const numeric_index_type i,
375  const numeric_index_type j,
376  const T value)
377 {
378  libmesh_assert (this->initialized());
379 
380  int
381  epetra_i = static_cast<int>(i),
382  epetra_j = static_cast<int>(j);
383 
384  T epetra_value = value;
385 
386  if (_mat->Filled())
387  _mat->ReplaceGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
388  else
389  _mat->InsertGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
390 }
391 
392 
393 
394 template <typename T>
395 void EpetraMatrix<T>::add (const numeric_index_type i,
396  const numeric_index_type j,
397  const T value)
398 {
399  libmesh_assert (this->initialized());
400 
401  int
402  epetra_i = static_cast<int>(i),
403  epetra_j = static_cast<int>(j);
404 
405  T epetra_value = value;
406 
407  _mat->SumIntoGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
408 }
409 
410 
411 
412 template <typename T>
413 void EpetraMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
414  const std::vector<numeric_index_type> & dof_indices)
415 {
416  this->add_matrix (dm, dof_indices, dof_indices);
417 }
418 
419 
420 
421 template <typename T>
422 void EpetraMatrix<T>::add (const T a_in, SparseMatrix<T> & X_in)
423 {
424 #ifdef LIBMESH_TRILINOS_HAVE_EPETRAEXT
425  libmesh_assert (this->initialized());
426 
427  // sanity check. but this cannot avoid
428  // crash due to incompatible sparsity structure...
429  libmesh_assert_equal_to (this->m(), X_in.m());
430  libmesh_assert_equal_to (this->n(), X_in.n());
431 
432  EpetraMatrix<T> * X = cast_ptr<EpetraMatrix<T> *> (&X_in);
433 
434  EpetraExt::MatrixMatrix::Add (*X->_mat, false, a_in, *_mat, 1.);
435 #else
436  libmesh_error_msg("ERROR: EpetraExt is required for EpetraMatrix::add()!");
437 #endif
438 }
439 
440 
441 
442 
443 template <typename T>
444 T EpetraMatrix<T>::operator () (const numeric_index_type i,
445  const numeric_index_type j) const
446 {
447  libmesh_assert (this->initialized());
448  libmesh_assert(this->_mat);
449  libmesh_assert (this->_mat->MyGlobalRow(static_cast<int>(i)));
450  libmesh_assert_greater_equal (i, this->row_start());
451  libmesh_assert_less (i, this->row_stop());
452 
453 
454  int row_length;
455  int * row_indices;
456  double * values;
457 
458  _mat->ExtractMyRowView (i-this->row_start(),
459  row_length,
460  values,
461  row_indices);
462 
463  //libMesh::out << "row_length=" << row_length << std::endl;
464 
465  int * index = std::lower_bound (row_indices, row_indices+row_length, j);
466 
467  libmesh_assert_less (*index, row_length);
468  libmesh_assert_equal_to (static_cast<numeric_index_type>(row_indices[*index]), j);
469 
470  //libMesh::out << "val=" << values[*index] << std::endl;
471 
472  return values[*index];
473 }
474 
475 
476 
477 
478 template <typename T>
479 bool EpetraMatrix<T>::closed() const
480 {
481  libmesh_assert (this->initialized());
482  libmesh_assert(this->_mat);
483 
484  return this->_mat->Filled();
485 }
486 
487 
488 template <typename T>
489 void EpetraMatrix<T>::swap(EpetraMatrix<T> & m)
490 {
491  std::swap(_mat, m._mat);
492  std::swap(_destroy_mat_on_exit, m._destroy_mat_on_exit);
493 }
494 
495 
496 
497 
498 
499 template <typename T>
500 void EpetraMatrix<T>::print_personal(std::ostream & os) const
501 {
502  libmesh_assert (this->initialized());
503  libmesh_assert(_mat);
504 
505  os << *_mat;
506 }
507 
508 
509 
510 //------------------------------------------------------------------
511 // Explicit instantiations
512 template class EpetraMatrix<Number>;
513 
514 } // namespace libMesh
515 
516 
517 #endif // LIBMESH_TRILINOS_HAVE_EPETRA
bool closed()
Checks that the library has been closed.
Definition: libmesh.C:281
The libMesh namespace provides an interface to certain functionality in the library.
const Number zero
.
Definition: libmesh.h:178
libmesh_assert(j)
dof_id_type numeric_index_type
Definition: id_types.h:92
bool _is_initialized
Flag that tells if init() has been called.
Definition: libmesh.C:255
void init(triangulateio &t)
Initializes the fields of t to NULL/0 as necessary.
PetscErrorCode Vec Mat libmesh_dbg_var(j)
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
void swap(Iterator &lhs, Iterator &rhs)
swap, used to implement op=
static const bool value
Definition: xdr_io.C:108
bool initialized()
Checks that library initialization has been done.
Definition: libmesh.C:274
processor_id_type processor_id()
Definition: libmesh_base.h:96
long double min(long double a, double b)