trilinos_epetra_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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_HAVE_TRILINOS
23 
24 // Local includes
27 #include "libmesh/dof_map.h"
28 #include "libmesh/dense_matrix.h"
29 #include "libmesh/parallel.h"
31 
32 namespace libMesh
33 {
34 
35 
36 
37 //-----------------------------------------------------------------------
38 //EpetraMatrix members
39 
40 template <typename T>
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,
122  const numeric_index_type n_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>
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>
195 {
196  libmesh_assert (this->initialized());
197 
198  _mat->Scale(0.0);
199 }
200 
201 
202 
203 template <typename T>
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>
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>
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>::print_matlab (const std::string) const
243 {
244  libmesh_assert (this->initialized());
245 
246  // libmesh_assert (this->closed());
247  this->close();
248 
249  libmesh_not_implemented();
250 }
251 
252 
253 
254 
255 template <typename T>
257  const std::vector<numeric_index_type>& rows,
258  const std::vector<numeric_index_type>& cols)
259 {
260  libmesh_assert (this->initialized());
261 
262  const numeric_index_type m = dm.m();
263  const numeric_index_type n = dm.n();
264 
265  libmesh_assert_equal_to (rows.size(), m);
266  libmesh_assert_equal_to (cols.size(), n);
267 
268  _mat->SumIntoGlobalValues(m, (int *)&rows[0], n, (int *)&cols[0], &dm.get_values()[0]);
269 }
270 
271 
272 
273 
274 
275 
276 template <typename T>
278 {
279  // Convert vector to EpetraVector.
280  EpetraVector<T>* epetra_dest = libmesh_cast_ptr<EpetraVector<T>*>(&dest);
281 
282  // Call Epetra function.
283  _mat->ExtractDiagonalCopy(*(epetra_dest->vec()));
284 }
285 
286 
287 
288 template <typename T>
290 {
291  // Make sure the SparseMatrix passed in is really a EpetraMatrix
292  EpetraMatrix<T>& epetra_dest = libmesh_cast_ref<EpetraMatrix<T>&>(dest);
293 
294  if(&epetra_dest != this)
295  epetra_dest = *this;
296 
297  epetra_dest._use_transpose = !epetra_dest._use_transpose;
298  epetra_dest._mat->SetUseTranspose(epetra_dest._use_transpose);
299 }
300 
301 
302 
303 template <typename T>
305  : SparseMatrix<T>(comm),
306  _destroy_mat_on_exit(true),
307  _use_transpose(false)
308 {}
309 
310 
311 
312 
313 template <typename T>
314 EpetraMatrix<T>::EpetraMatrix(Epetra_FECrsMatrix * m,
315  const Parallel::Communicator &comm)
316  : SparseMatrix<T>(comm),
317  _destroy_mat_on_exit(false),
318  _use_transpose(false) // dumb guess is the best we can do...
319 {
320  this->_mat = m;
321  this->_is_initialized = true;
322 }
323 
324 
325 
326 
327 template <typename T>
329 {
330  this->clear();
331 }
332 
333 
334 
335 template <typename T>
337 {
338  libmesh_assert(_mat);
339 
340  _mat->GlobalAssemble();
341 }
342 
343 
344 
345 template <typename T>
347 {
348  libmesh_assert (this->initialized());
349 
350  return static_cast<numeric_index_type>(_mat->NumGlobalRows());
351 }
352 
353 
354 
355 template <typename T>
357 {
358  libmesh_assert (this->initialized());
359 
360  return static_cast<numeric_index_type>(_mat->NumGlobalCols());
361 }
362 
363 
364 
365 template <typename T>
367 {
368  libmesh_assert (this->initialized());
369  libmesh_assert(_map);
370 
371  return static_cast<numeric_index_type>(_map->MinMyGID());
372 }
373 
374 
375 
376 template <typename T>
378 {
379  libmesh_assert (this->initialized());
380  libmesh_assert(_map);
381 
382  return static_cast<numeric_index_type>(_map->MaxMyGID())+1;
383 }
384 
385 
386 
387 template <typename T>
389  const numeric_index_type j,
390  const T value)
391 {
392  libmesh_assert (this->initialized());
393 
394  int
395  epetra_i = static_cast<int>(i),
396  epetra_j = static_cast<int>(j);
397 
398  T epetra_value = value;
399 
400  if (_mat->Filled())
401  _mat->ReplaceGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
402  else
403  _mat->InsertGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
404 }
405 
406 
407 
408 template <typename T>
410  const numeric_index_type j,
411  const T value)
412 {
413  libmesh_assert (this->initialized());
414 
415  int
416  epetra_i = static_cast<int>(i),
417  epetra_j = static_cast<int>(j);
418 
419  T epetra_value = value;
420 
421  _mat->SumIntoGlobalValues (epetra_i, 1, &epetra_value, &epetra_j);
422 }
423 
424 
425 
426 template <typename T>
428  const std::vector<numeric_index_type>& dof_indices)
429 {
430  this->add_matrix (dm, dof_indices, dof_indices);
431 }
432 
433 
434 
435 template <typename T>
436 void EpetraMatrix<T>::add (const T a_in, SparseMatrix<T> &X_in)
437 {
438  libmesh_assert (this->initialized());
439 
440  // sanity check. but this cannot avoid
441  // crash due to incompatible sparsity structure...
442  libmesh_assert_equal_to (this->m(), X_in.m());
443  libmesh_assert_equal_to (this->n(), X_in.n());
444 
445  EpetraMatrix<T>* X = libmesh_cast_ptr<EpetraMatrix<T>*> (&X_in);
446 
447  EpetraExt::MatrixMatrix::Add (*X->_mat, false, a_in, *_mat, 1.);
448 }
449 
450 
451 
452 
453 template <typename T>
455  const numeric_index_type j) const
456 {
457  libmesh_assert (this->initialized());
458  libmesh_assert(this->_mat);
459  libmesh_assert (this->_mat->MyGlobalRow(static_cast<int>(i)));
460  libmesh_assert_greater_equal (i, this->row_start());
461  libmesh_assert_less (i, this->row_stop());
462 
463 
464  int row_length, *row_indices;
465  double *values;
466 
467  _mat->ExtractMyRowView (i-this->row_start(),
468  row_length,
469  values,
470  row_indices);
471 
472  //libMesh::out << "row_length=" << row_length << std::endl;
473 
474  int *index = std::lower_bound (row_indices, row_indices+row_length, j);
475 
476  libmesh_assert_less (*index, row_length);
477  libmesh_assert_equal_to (static_cast<numeric_index_type>(row_indices[*index]), j);
478 
479  //libMesh::out << "val=" << values[*index] << std::endl;
480 
481  return values[*index];
482 }
483 
484 
485 
486 
487 template <typename T>
489 {
490  libmesh_assert (this->initialized());
491  libmesh_assert(this->_mat);
492 
493  return this->_mat->Filled();
494 }
495 
496 
497 template <typename T>
499 {
500  std::swap(_mat, m._mat);
501  std::swap(_destroy_mat_on_exit, m._destroy_mat_on_exit);
502 }
503 
504 
505 
506 
507 
508 template <typename T>
509 void EpetraMatrix<T>::print_personal(std::ostream& os) const
510 {
511  libmesh_assert (this->initialized());
512  libmesh_assert(_mat);
513 
514  os << *_mat;
515 }
516 
517 
518 
519 //------------------------------------------------------------------
520 // Explicit instantiations
521 template class EpetraMatrix<Number>;
522 
523 } // namespace libMesh
524 
525 
526 #endif // #ifdef LIBMESH_HAVE_TRILINOS

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:07 UTC

Hosted By:
SourceForge.net Logo