petsc_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2018 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 <unistd.h> // mkstemp
21 #include <fstream>
22 
23 #include "libmesh/libmesh_config.h"
24 
25 #ifdef LIBMESH_HAVE_PETSC
26 
27 // Local includes
28 #include "libmesh/petsc_matrix.h"
29 #include "libmesh/dof_map.h"
30 #include "libmesh/dense_matrix.h"
31 #include "libmesh/petsc_vector.h"
32 
33 // For some reason, the blocked matrix API calls below seem to break with PETSC 3.2 & presumably earlier.
34 // For example:
35 // [0]PETSC ERROR: --------------------- Error Message ------------------------------------
36 // [0]PETSC ERROR: Nonconforming object sizes!
37 // [0]PETSC ERROR: Attempt to set block size 3 with BAIJ 1!
38 // [0]PETSC ERROR: ------------------------------------------------------------------------
39 // so as a cowardly workaround disable the functionality in this translation unit for older PETSc's
40 #if PETSC_VERSION_LESS_THAN(3,3,0)
41 # undef LIBMESH_ENABLE_BLOCKED_STORAGE
42 #endif
43 
44 
45 
46 #ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
47 
48 namespace
49 {
50 using namespace libMesh;
51 
52 // historic libMesh n_nz & n_oz arrays are set up for PETSc's AIJ format.
53 // however, when the blocksize is >1, we need to transform these into
54 // their BAIJ counterparts.
55 inline
56 void transform_preallocation_arrays (const PetscInt blocksize,
57  const std::vector<numeric_index_type> & n_nz,
58  const std::vector<numeric_index_type> & n_oz,
59  std::vector<numeric_index_type> & b_n_nz,
60  std::vector<numeric_index_type> & b_n_oz)
61 {
62  libmesh_assert_equal_to (n_nz.size(), n_oz.size());
63  libmesh_assert_equal_to (n_nz.size()%blocksize, 0);
64 
65  b_n_nz.clear(); b_n_nz.reserve(n_nz.size()/blocksize);
66  b_n_oz.clear(); b_n_oz.reserve(n_oz.size()/blocksize);
67 
68  for (std::size_t nn=0; nn<n_nz.size(); nn += blocksize)
69  {
70  b_n_nz.push_back (n_nz[nn]/blocksize);
71  b_n_oz.push_back (n_oz[nn]/blocksize);
72  }
73 }
74 }
75 
76 #endif
77 
78 
79 
80 namespace libMesh
81 {
82 
83 
84 //-----------------------------------------------------------------------
85 // PetscMatrix members
86 
87 
88 // Constructor
89 template <typename T>
91  SparseMatrix<T>(comm_in),
92  _destroy_mat_on_exit(true),
93  _mat_type(AIJ)
94 {}
95 
96 
97 
98 // Constructor taking an existing Mat but not the responsibility
99 // for destroying it
100 template <typename T>
102  const Parallel::Communicator & comm_in) :
103  SparseMatrix<T>(comm_in),
104  _destroy_mat_on_exit(false),
105  _mat_type(AIJ)
106 {
107  this->_mat = mat_in;
108  this->_is_initialized = true;
109 }
110 
111 
112 
113 // Destructor
114 template <typename T>
116 {
117  this->clear();
118 }
119 
120 template <typename T>
122 {
123  _mat_type = mat_type;
124 }
125 
126 template <typename T>
128  const numeric_index_type n_in,
129  const numeric_index_type m_l,
130  const numeric_index_type n_l,
131  const numeric_index_type nnz,
132  const numeric_index_type noz,
133  const numeric_index_type blocksize_in)
134 {
135  // So compilers don't warn when !LIBMESH_ENABLE_BLOCKED_STORAGE
136  libmesh_ignore(blocksize_in);
137 
138  // Clear initialized matrices
139  if (this->initialized())
140  this->clear();
141 
142  this->_is_initialized = true;
143 
144 
145  PetscErrorCode ierr = 0;
146  PetscInt m_global = static_cast<PetscInt>(m_in);
147  PetscInt n_global = static_cast<PetscInt>(n_in);
148  PetscInt m_local = static_cast<PetscInt>(m_l);
149  PetscInt n_local = static_cast<PetscInt>(n_l);
150  PetscInt n_nz = static_cast<PetscInt>(nnz);
151  PetscInt n_oz = static_cast<PetscInt>(noz);
152 
153  ierr = MatCreate(this->comm().get(), &_mat);
154  LIBMESH_CHKERR(ierr);
155  ierr = MatSetSizes(_mat, m_local, n_local, m_global, n_global);
156  LIBMESH_CHKERR(ierr);
157 
158 #ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
159  PetscInt blocksize = static_cast<PetscInt>(blocksize_in);
160  if (blocksize > 1)
161  {
162  // specified blocksize, bs>1.
163  // double check sizes.
164  libmesh_assert_equal_to (m_local % blocksize, 0);
165  libmesh_assert_equal_to (n_local % blocksize, 0);
166  libmesh_assert_equal_to (m_global % blocksize, 0);
167  libmesh_assert_equal_to (n_global % blocksize, 0);
168  libmesh_assert_equal_to (n_nz % blocksize, 0);
169  libmesh_assert_equal_to (n_oz % blocksize, 0);
170 
171  ierr = MatSetType(_mat, MATBAIJ); // Automatically chooses seqbaij or mpibaij
172  LIBMESH_CHKERR(ierr);
173  ierr = MatSetBlockSize(_mat, blocksize);
174  LIBMESH_CHKERR(ierr);
175  ierr = MatSeqBAIJSetPreallocation(_mat, blocksize, n_nz/blocksize, PETSC_NULL);
176  LIBMESH_CHKERR(ierr);
177  ierr = MatMPIBAIJSetPreallocation(_mat, blocksize,
178  n_nz/blocksize, PETSC_NULL,
179  n_oz/blocksize, PETSC_NULL);
180  LIBMESH_CHKERR(ierr);
181  }
182  else
183 #endif
184  {
185  switch (_mat_type) {
186  case AIJ:
187  ierr = MatSetType(_mat, MATAIJ); // Automatically chooses seqaij or mpiaij
188  LIBMESH_CHKERR(ierr);
189  ierr = MatSeqAIJSetPreallocation(_mat, n_nz, PETSC_NULL);
190  LIBMESH_CHKERR(ierr);
191  ierr = MatMPIAIJSetPreallocation(_mat, n_nz, PETSC_NULL, n_oz, PETSC_NULL);
192  LIBMESH_CHKERR(ierr);
193  break;
194 
195  case HYPRE:
196 #if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
197  ierr = MatSetType(_mat, MATHYPRE);
198  LIBMESH_CHKERR(ierr);
199  ierr = MatHYPRESetPreallocation(_mat, n_nz, PETSC_NULL, n_oz, PETSC_NULL);
200  LIBMESH_CHKERR(ierr);
201 #else
202  libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
203 #endif
204  break;
205 
206  default: libmesh_error_msg("Unsupported petsc matrix type");
207  }
208  }
209 
210  // Make it an error for PETSc to allocate new nonzero entries during assembly
211 #if PETSC_VERSION_LESS_THAN(3,0,0)
212  ierr = MatSetOption(_mat, MAT_NEW_NONZERO_ALLOCATION_ERR);
213 #else
214  ierr = MatSetOption(_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
215 #endif
216  LIBMESH_CHKERR(ierr);
217 
218  // Is prefix information available somewhere? Perhaps pass in the system name?
219  ierr = MatSetOptionsPrefix(_mat, "");
220  LIBMESH_CHKERR(ierr);
221  ierr = MatSetFromOptions(_mat);
222  LIBMESH_CHKERR(ierr);
223 
224  this->zero ();
225 }
226 
227 
228 template <typename T>
230  const numeric_index_type n_in,
231  const numeric_index_type m_l,
232  const numeric_index_type n_l,
233  const std::vector<numeric_index_type> & n_nz,
234  const std::vector<numeric_index_type> & n_oz,
235  const numeric_index_type blocksize_in)
236 {
237  // So compilers don't warn when !LIBMESH_ENABLE_BLOCKED_STORAGE
238  libmesh_ignore(blocksize_in);
239 
240  // Clear initialized matrices
241  if (this->initialized())
242  this->clear();
243 
244  this->_is_initialized = true;
245 
246  // Make sure the sparsity pattern isn't empty unless the matrix is 0x0
247  libmesh_assert_equal_to (n_nz.size(), m_l);
248  libmesh_assert_equal_to (n_oz.size(), m_l);
249 
250  PetscErrorCode ierr = 0;
251  PetscInt m_global = static_cast<PetscInt>(m_in);
252  PetscInt n_global = static_cast<PetscInt>(n_in);
253  PetscInt m_local = static_cast<PetscInt>(m_l);
254  PetscInt n_local = static_cast<PetscInt>(n_l);
255 
256  ierr = MatCreate(this->comm().get(), &_mat);
257  LIBMESH_CHKERR(ierr);
258  ierr = MatSetSizes(_mat, m_local, n_local, m_global, n_global);
259  LIBMESH_CHKERR(ierr);
260 
261 #ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
262  PetscInt blocksize = static_cast<PetscInt>(blocksize_in);
263  if (blocksize > 1)
264  {
265  // specified blocksize, bs>1.
266  // double check sizes.
267  libmesh_assert_equal_to (m_local % blocksize, 0);
268  libmesh_assert_equal_to (n_local % blocksize, 0);
269  libmesh_assert_equal_to (m_global % blocksize, 0);
270  libmesh_assert_equal_to (n_global % blocksize, 0);
271 
272  ierr = MatSetType(_mat, MATBAIJ); // Automatically chooses seqbaij or mpibaij
273  LIBMESH_CHKERR(ierr);
274  ierr = MatSetBlockSize(_mat, blocksize);
275  LIBMESH_CHKERR(ierr);
276 
277  // transform the per-entry n_nz and n_oz arrays into their block counterparts.
278  std::vector<numeric_index_type> b_n_nz, b_n_oz;
279 
280  transform_preallocation_arrays (blocksize,
281  n_nz, n_oz,
282  b_n_nz, b_n_oz);
283 
284  ierr = MatSeqBAIJSetPreallocation (_mat,
285  blocksize,
286  0,
287  numeric_petsc_cast(b_n_nz.empty() ? nullptr : b_n_nz.data()));
288  LIBMESH_CHKERR(ierr);
289 
290  ierr = MatMPIBAIJSetPreallocation (_mat,
291  blocksize,
292  0,
293  numeric_petsc_cast(b_n_nz.empty() ? nullptr : b_n_nz.data()),
294  0,
295  numeric_petsc_cast(b_n_oz.empty() ? nullptr : b_n_oz.data()));
296  LIBMESH_CHKERR(ierr);
297  }
298  else
299 #endif
300  {
301  switch (_mat_type) {
302  case AIJ:
303  ierr = MatSetType(_mat, MATAIJ); // Automatically chooses seqaij or mpiaij
304  LIBMESH_CHKERR(ierr);
305  ierr = MatSeqAIJSetPreallocation (_mat,
306  0,
307  numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()));
308  LIBMESH_CHKERR(ierr);
309  ierr = MatMPIAIJSetPreallocation (_mat,
310  0,
311  numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
312  0,
313  numeric_petsc_cast(n_oz.empty() ? nullptr : n_oz.data()));
314  break;
315 
316  case HYPRE:
317 #if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
318  ierr = MatSetType(_mat, MATHYPRE);
319  LIBMESH_CHKERR(ierr);
320  ierr = MatHYPRESetPreallocation (_mat,
321  0,
322  numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
323  0,
324  numeric_petsc_cast(n_oz.empty() ? nullptr : n_oz.data()));
325  LIBMESH_CHKERR(ierr);
326 #else
327  libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
328 #endif
329  break;
330 
331  default: libmesh_error_msg("Unsupported petsc matrix type");
332  }
333 
334  }
335 
336  // Is prefix information available somewhere? Perhaps pass in the system name?
337  ierr = MatSetOptionsPrefix(_mat, "");
338  LIBMESH_CHKERR(ierr);
339  ierr = MatSetFromOptions(_mat);
340  LIBMESH_CHKERR(ierr);
341 
342 
343  this->zero();
344 }
345 
346 
347 template <typename T>
349 {
350  libmesh_assert(this->_dof_map);
351 
352  // Clear initialized matrices
353  if (this->initialized())
354  this->clear();
355 
356  this->_is_initialized = true;
357 
358 
359  const numeric_index_type my_m = this->_dof_map->n_dofs();
360  const numeric_index_type my_n = my_m;
361  const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(this->processor_id());
362  const numeric_index_type m_l = n_l;
363 
364 
365  const std::vector<numeric_index_type> & n_nz = this->_dof_map->get_n_nz();
366  const std::vector<numeric_index_type> & n_oz = this->_dof_map->get_n_oz();
367 
368  // Make sure the sparsity pattern isn't empty unless the matrix is 0x0
369  libmesh_assert_equal_to (n_nz.size(), m_l);
370  libmesh_assert_equal_to (n_oz.size(), m_l);
371 
372  PetscErrorCode ierr = 0;
373  PetscInt m_global = static_cast<PetscInt>(my_m);
374  PetscInt n_global = static_cast<PetscInt>(my_n);
375  PetscInt m_local = static_cast<PetscInt>(m_l);
376  PetscInt n_local = static_cast<PetscInt>(n_l);
377 
378  ierr = MatCreate(this->comm().get(), &_mat);
379  LIBMESH_CHKERR(ierr);
380  ierr = MatSetSizes(_mat, m_local, n_local, m_global, n_global);
381  LIBMESH_CHKERR(ierr);
382 
383 #ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
384  PetscInt blocksize = static_cast<PetscInt>(this->_dof_map->block_size());
385  if (blocksize > 1)
386  {
387  // specified blocksize, bs>1.
388  // double check sizes.
389  libmesh_assert_equal_to (m_local % blocksize, 0);
390  libmesh_assert_equal_to (n_local % blocksize, 0);
391  libmesh_assert_equal_to (m_global % blocksize, 0);
392  libmesh_assert_equal_to (n_global % blocksize, 0);
393 
394  ierr = MatSetType(_mat, MATBAIJ); // Automatically chooses seqbaij or mpibaij
395  LIBMESH_CHKERR(ierr);
396  ierr = MatSetBlockSize(_mat, blocksize);
397  LIBMESH_CHKERR(ierr);
398 
399  // transform the per-entry n_nz and n_oz arrays into their block counterparts.
400  std::vector<numeric_index_type> b_n_nz, b_n_oz;
401 
402  transform_preallocation_arrays (blocksize,
403  n_nz, n_oz,
404  b_n_nz, b_n_oz);
405 
406  ierr = MatSeqBAIJSetPreallocation (_mat,
407  blocksize,
408  0,
409  numeric_petsc_cast(b_n_nz.empty() ? nullptr : b_n_nz.data()));
410  LIBMESH_CHKERR(ierr);
411 
412  ierr = MatMPIBAIJSetPreallocation (_mat,
413  blocksize,
414  0,
415  numeric_petsc_cast(b_n_nz.empty() ? nullptr : b_n_nz.data()),
416  0,
417  numeric_petsc_cast(b_n_oz.empty() ? nullptr : b_n_oz.data()));
418  LIBMESH_CHKERR(ierr);
419  }
420  else
421 #endif
422  {
423  switch (_mat_type) {
424  case AIJ:
425  ierr = MatSetType(_mat, MATAIJ); // Automatically chooses seqaij or mpiaij
426  LIBMESH_CHKERR(ierr);
427  ierr = MatSeqAIJSetPreallocation (_mat,
428  0,
429  numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()));
430  LIBMESH_CHKERR(ierr);
431  ierr = MatMPIAIJSetPreallocation (_mat,
432  0,
433  numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
434  0,
435  numeric_petsc_cast(n_oz.empty() ? nullptr : n_oz.data()));
436  break;
437 
438  case HYPRE:
439 #if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
440  ierr = MatSetType(_mat, MATHYPRE);
441  LIBMESH_CHKERR(ierr);
442  ierr = MatHYPRESetPreallocation (_mat,
443  0,
444  numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
445  0,
446  numeric_petsc_cast(n_oz.empty() ? nullptr : n_oz.data()));
447  LIBMESH_CHKERR(ierr);
448 #else
449  libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
450 #endif
451  break;
452 
453  default: libmesh_error_msg("Unsupported petsc matrix type");
454  }
455  }
456 
457  // Is prefix information available somewhere? Perhaps pass in the system name?
458  ierr = MatSetOptionsPrefix(_mat, "");
459  LIBMESH_CHKERR(ierr);
460  ierr = MatSetFromOptions(_mat);
461  LIBMESH_CHKERR(ierr);
462 
463  this->zero();
464 }
465 
466 
467 template <typename T>
469 {
470  libmesh_not_implemented();
471 }
472 
473 
474 template <typename T>
476 {
477  libmesh_assert (this->initialized());
478 
479  semiparallel_only();
480 
481  PetscErrorCode ierr=0;
482 
483  PetscInt m_l, n_l;
484 
485  ierr = MatGetLocalSize(_mat,&m_l,&n_l);
486  LIBMESH_CHKERR(ierr);
487 
488  if (n_l)
489  {
490  ierr = MatZeroEntries(_mat);
491  LIBMESH_CHKERR(ierr);
492  }
493 }
494 
495 template <typename T>
496 void PetscMatrix<T>::zero_rows (std::vector<numeric_index_type> & rows, T diag_value)
497 {
498  libmesh_assert (this->initialized());
499 
500  semiparallel_only();
501 
502  PetscErrorCode ierr=0;
503 
504 #if PETSC_RELEASE_LESS_THAN(3,1,1)
505  if (!rows.empty())
506  ierr = MatZeroRows(_mat, rows.size(),
507  numeric_petsc_cast(rows.data()), diag_value);
508  else
509  ierr = MatZeroRows(_mat, 0, PETSC_NULL, diag_value);
510 #else
511  // As of petsc-dev at the time of 3.1.0, MatZeroRows now takes two additional
512  // optional arguments. The optional arguments (x,b) can be used to specify the
513  // solutions for the zeroed rows (x) and right hand side (b) to update.
514  // Could be useful for setting boundary conditions...
515  if (!rows.empty())
516  ierr = MatZeroRows(_mat, cast_int<PetscInt>(rows.size()),
517  numeric_petsc_cast(rows.data()), diag_value,
518  PETSC_NULL, PETSC_NULL);
519  else
520  ierr = MatZeroRows(_mat, 0, PETSC_NULL, diag_value, PETSC_NULL,
521  PETSC_NULL);
522 #endif
523 
524  LIBMESH_CHKERR(ierr);
525 }
526 
527 template <typename T>
529 {
530  PetscErrorCode ierr=0;
531 
532  if ((this->initialized()) && (this->_destroy_mat_on_exit))
533  {
534  semiparallel_only();
535 
536  ierr = LibMeshMatDestroy (&_mat);
537  LIBMESH_CHKERR(ierr);
538 
539  this->_is_initialized = false;
540  }
541 }
542 
543 
544 
545 template <typename T>
547 {
548  libmesh_assert (this->initialized());
549 
550  semiparallel_only();
551 
552  PetscErrorCode ierr=0;
553  PetscReal petsc_value;
554  Real value;
555 
556  libmesh_assert (this->closed());
557 
558  ierr = MatNorm(_mat, NORM_1, &petsc_value);
559  LIBMESH_CHKERR(ierr);
560 
561  value = static_cast<Real>(petsc_value);
562 
563  return value;
564 }
565 
566 
567 
568 template <typename T>
570 {
571  libmesh_assert (this->initialized());
572 
573  semiparallel_only();
574 
575  PetscErrorCode ierr=0;
576  PetscReal petsc_value;
577  Real value;
578 
579  libmesh_assert (this->closed());
580 
581  ierr = MatNorm(_mat, NORM_INFINITY, &petsc_value);
582  LIBMESH_CHKERR(ierr);
583 
584  value = static_cast<Real>(petsc_value);
585 
586  return value;
587 }
588 
589 
590 
591 template <typename T>
592 void PetscMatrix<T>::print_matlab (const std::string & name) const
593 {
594  libmesh_assert (this->initialized());
595 
596  semiparallel_only();
597 
598  if (!this->closed())
599  {
600  libmesh_deprecated();
601  libmesh_warning("The matrix must be assembled before calling PetscMatrix::print_matlab().\n"
602  "Please update your code, as this warning will become an error in a future release.");
603  const_cast<PetscMatrix<T> *>(this)->close();
604  }
605 
606  PetscErrorCode ierr=0;
607  PetscViewer petsc_viewer;
608 
609 
610  ierr = PetscViewerCreate (this->comm().get(),
611  &petsc_viewer);
612  LIBMESH_CHKERR(ierr);
613 
618  if (name != "")
619  {
620  ierr = PetscViewerASCIIOpen( this->comm().get(),
621  name.c_str(),
622  &petsc_viewer);
623  LIBMESH_CHKERR(ierr);
624 
625 #if PETSC_VERSION_LESS_THAN(3,7,0)
626  ierr = PetscViewerSetFormat (petsc_viewer,
627  PETSC_VIEWER_ASCII_MATLAB);
628 #else
629  ierr = PetscViewerPushFormat (petsc_viewer,
630  PETSC_VIEWER_ASCII_MATLAB);
631 #endif
632 
633  LIBMESH_CHKERR(ierr);
634 
635  ierr = MatView (_mat, petsc_viewer);
636  LIBMESH_CHKERR(ierr);
637  }
638 
642  else
643  {
644 #if PETSC_VERSION_LESS_THAN(3,7,0)
645  ierr = PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
646  PETSC_VIEWER_ASCII_MATLAB);
647 #else
648  ierr = PetscViewerPushFormat (PETSC_VIEWER_STDOUT_WORLD,
649  PETSC_VIEWER_ASCII_MATLAB);
650 #endif
651 
652  LIBMESH_CHKERR(ierr);
653 
654  ierr = MatView (_mat, PETSC_VIEWER_STDOUT_WORLD);
655  LIBMESH_CHKERR(ierr);
656  }
657 
658 
662  ierr = LibMeshPetscViewerDestroy (&petsc_viewer);
663  LIBMESH_CHKERR(ierr);
664 }
665 
666 
667 
668 
669 
670 template <typename T>
671 void PetscMatrix<T>::print_personal(std::ostream & os) const
672 {
673  libmesh_assert (this->initialized());
674 
675  // Routine must be called in parallel on parallel matrices
676  // and serial on serial matrices.
677  semiparallel_only();
678 
679  // #ifndef NDEBUG
680  // if (os != std::cout)
681  // libMesh::err << "Warning! PETSc can only print to std::cout!" << std::endl;
682  // #endif
683 
684  // Matrix must be in an assembled state to be printed
685  if (!this->closed())
686  {
687  libmesh_deprecated();
688  libmesh_warning("The matrix must be assembled before calling PetscMatrix::print_personal().\n"
689  "Please update your code, as this warning will become an error in a future release.");
690  const_cast<PetscMatrix<T> *>(this)->close();
691  }
692 
693  PetscErrorCode ierr=0;
694 
695  // Print to screen if ostream is stdout
696  if (os.rdbuf() == std::cout.rdbuf())
697  {
698  ierr = MatView(_mat, PETSC_VIEWER_STDOUT_SELF);
699  LIBMESH_CHKERR(ierr);
700  }
701 
702  // Otherwise, print to the requested file, in a roundabout way...
703  else
704  {
705  // We will create a temporary filename, and file, for PETSc to
706  // write to.
707  std::string temp_filename;
708 
709  {
710  // Template for temporary filename
711  char c[] = "temp_petsc_matrix.XXXXXX";
712 
713  // Generate temporary, unique filename only on processor 0. We will
714  // use this filename for PetscViewerASCIIOpen, before copying it into
715  // the user's stream
716  if (this->processor_id() == 0)
717  {
718  int fd = mkstemp(c);
719 
720  // Check to see that mkstemp did not fail.
721  if (fd == -1)
722  libmesh_error_msg("mkstemp failed in PetscMatrix::print_personal()");
723 
724  // mkstemp returns a file descriptor for an open file,
725  // so let's close it before we hand it to PETSc!
726  ::close (fd);
727  }
728 
729  // Store temporary filename as string, makes it easier to broadcast
730  temp_filename = c;
731  }
732 
733  // Now broadcast the filename from processor 0 to all processors.
734  this->comm().broadcast(temp_filename);
735 
736  // PetscViewer object for passing to MatView
737  PetscViewer petsc_viewer;
738 
739  // This PETSc function only takes a string and handles the opening/closing
740  // of the file internally. Since print_personal() takes a reference to
741  // an ostream, we have to do an extra step... print_personal() should probably
742  // have a version that takes a string to get rid of this problem.
743  ierr = PetscViewerASCIIOpen( this->comm().get(),
744  temp_filename.c_str(),
745  &petsc_viewer);
746  LIBMESH_CHKERR(ierr);
747 
748  // Probably don't need to set the format if it's default...
749  // ierr = PetscViewerSetFormat (petsc_viewer,
750  // PETSC_VIEWER_DEFAULT);
751  // LIBMESH_CHKERR(ierr);
752 
753  // Finally print the matrix using the viewer
754  ierr = MatView (_mat, petsc_viewer);
755  LIBMESH_CHKERR(ierr);
756 
757  if (this->processor_id() == 0)
758  {
759  // Now the inefficient bit: open temp_filename as an ostream and copy the contents
760  // into the user's desired ostream. We can't just do a direct file copy, we don't even have the filename!
761  std::ifstream input_stream(temp_filename.c_str());
762  os << input_stream.rdbuf(); // The "most elegant" way to copy one stream into another.
763  // os.close(); // close not defined in ostream
764 
765  // Now remove the temporary file
766  input_stream.close();
767  std::remove(temp_filename.c_str());
768  }
769  }
770 }
771 
772 
773 
774 
775 
776 
777 template <typename T>
779  const std::vector<numeric_index_type> & rows,
780  const std::vector<numeric_index_type> & cols)
781 {
782  libmesh_assert (this->initialized());
783 
784  const numeric_index_type n_rows = dm.m();
785  const numeric_index_type n_cols = dm.n();
786 
787  libmesh_assert_equal_to (rows.size(), n_rows);
788  libmesh_assert_equal_to (cols.size(), n_cols);
789 
790  PetscErrorCode ierr=0;
791  ierr = MatSetValues(_mat,
792  n_rows, numeric_petsc_cast(rows.data()),
793  n_cols, numeric_petsc_cast(cols.data()),
794  const_cast<PetscScalar *>(dm.get_values().data()),
795  ADD_VALUES);
796  LIBMESH_CHKERR(ierr);
797 }
798 
799 
800 
801 
802 
803 
804 template <typename T>
806  const std::vector<numeric_index_type> & brows,
807  const std::vector<numeric_index_type> & bcols)
808 {
809  libmesh_assert (this->initialized());
810 
811  const numeric_index_type n_brows =
812  cast_int<numeric_index_type>(brows.size());
813  const numeric_index_type n_bcols =
814  cast_int<numeric_index_type>(bcols.size());
815 
816  PetscErrorCode ierr=0;
817 
818 #ifndef NDEBUG
819  const numeric_index_type n_rows =
820  cast_int<numeric_index_type>(dm.m());
821  const numeric_index_type n_cols =
822  cast_int<numeric_index_type>(dm.n());
823  const numeric_index_type blocksize = n_rows / n_brows;
824 
825  libmesh_assert_equal_to (n_cols / n_bcols, blocksize);
826  libmesh_assert_equal_to (blocksize*n_brows, n_rows);
827  libmesh_assert_equal_to (blocksize*n_bcols, n_cols);
828 
829  PetscInt petsc_blocksize;
830  ierr = MatGetBlockSize(_mat, &petsc_blocksize);
831  LIBMESH_CHKERR(ierr);
832  libmesh_assert_equal_to (blocksize, static_cast<numeric_index_type>(petsc_blocksize));
833 #endif
834 
835  // These casts are required for PETSc <= 2.1.5
836  ierr = MatSetValuesBlocked(_mat,
837  n_brows, numeric_petsc_cast(brows.data()),
838  n_bcols, numeric_petsc_cast(bcols.data()),
839  const_cast<PetscScalar *>(dm.get_values().data()),
840  ADD_VALUES);
841  LIBMESH_CHKERR(ierr);
842 }
843 
844 
845 
846 
847 
848 template <typename T>
850  const std::vector<numeric_index_type> & rows,
851  const std::vector<numeric_index_type> & cols,
852  const bool reuse_submatrix) const
853 {
854  if (!this->closed())
855  {
856  libmesh_deprecated();
857  libmesh_warning("The matrix must be assembled before calling PetscMatrix::create_submatrix().\n"
858  "Please update your code, as this warning will become an error in a future release.");
859  const_cast<PetscMatrix<T> *>(this)->close();
860  }
861 
862  // Make sure the SparseMatrix passed in is really a PetscMatrix
863  PetscMatrix<T> * petsc_submatrix = cast_ptr<PetscMatrix<T> *>(&submatrix);
864 
865  // If we're not reusing submatrix and submatrix is already initialized
866  // then we need to clear it, otherwise we get a memory leak.
867  if (!reuse_submatrix && submatrix.initialized())
868  submatrix.clear();
869 
870  // Construct row and column index sets.
871  PetscErrorCode ierr=0;
872  IS isrow, iscol;
873 
874  ierr = ISCreateLibMesh(this->comm().get(),
875  cast_int<PetscInt>(rows.size()),
876  numeric_petsc_cast(rows.data()),
878  &isrow); LIBMESH_CHKERR(ierr);
879 
880  ierr = ISCreateLibMesh(this->comm().get(),
881  cast_int<PetscInt>(cols.size()),
882  numeric_petsc_cast(cols.data()),
884  &iscol); LIBMESH_CHKERR(ierr);
885 
886  // Extract submatrix
887  ierr = LibMeshCreateSubMatrix(_mat,
888  isrow,
889  iscol,
890 #if PETSC_RELEASE_LESS_THAN(3,0,1)
891  PETSC_DECIDE,
892 #endif
893  (reuse_submatrix ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX),
894  &(petsc_submatrix->_mat)); LIBMESH_CHKERR(ierr);
895 
896  // Specify that the new submatrix is initialized and close it.
897  petsc_submatrix->_is_initialized = true;
898  petsc_submatrix->close();
899 
900  // Clean up PETSc data structures
901  ierr = LibMeshISDestroy(&isrow); LIBMESH_CHKERR(ierr);
902  ierr = LibMeshISDestroy(&iscol); LIBMESH_CHKERR(ierr);
903 }
904 
905 
906 
907 template <typename T>
909 {
910  // Make sure the NumericVector passed in is really a PetscVector
911  PetscVector<T> & petsc_dest = cast_ref<PetscVector<T> &>(dest);
912 
913  // Needs a const_cast since PETSc does not work with const.
914  PetscErrorCode ierr =
915  MatGetDiagonal(const_cast<PetscMatrix<T> *>(this)->mat(),petsc_dest.vec()); LIBMESH_CHKERR(ierr);
916 }
917 
918 
919 
920 template <typename T>
922 {
923  // Make sure the SparseMatrix passed in is really a PetscMatrix
924  PetscMatrix<T> & petsc_dest = cast_ref<PetscMatrix<T> &>(dest);
925 
926  // If we aren't reusing the matrix then need to clear dest,
927  // otherwise we get a memory leak
928  if (&petsc_dest != this)
929  dest.clear();
930 
931  PetscErrorCode ierr;
932 #if PETSC_VERSION_LESS_THAN(3,0,0)
933  if (&petsc_dest == this)
934  ierr = MatTranspose(_mat,PETSC_NULL);
935  else
936  ierr = MatTranspose(_mat,&petsc_dest._mat);
937  LIBMESH_CHKERR(ierr);
938 #else
939  // FIXME - we can probably use MAT_REUSE_MATRIX in more situations
940  if (&petsc_dest == this)
941  ierr = MatTranspose(_mat,MAT_REUSE_MATRIX,&petsc_dest._mat);
942  else
943  ierr = MatTranspose(_mat,MAT_INITIAL_MATRIX,&petsc_dest._mat);
944  LIBMESH_CHKERR(ierr);
945 #endif
946 
947  // Specify that the transposed matrix is initialized and close it.
948  petsc_dest._is_initialized = true;
949  petsc_dest.close();
950 }
951 
952 
953 
954 
955 
956 template <typename T>
958 {
959  semiparallel_only();
960 
961  // BSK - 1/19/2004
962  // strictly this check should be OK, but it seems to
963  // fail on matrix-free matrices. Do they falsely
964  // state they are assembled? Check with the developers...
965  // if (this->closed())
966  // return;
967 
968  PetscErrorCode ierr=0;
969 
970  ierr = MatAssemblyBegin (_mat, MAT_FINAL_ASSEMBLY);
971  LIBMESH_CHKERR(ierr);
972  ierr = MatAssemblyEnd (_mat, MAT_FINAL_ASSEMBLY);
973  LIBMESH_CHKERR(ierr);
974 }
975 
976 template <typename T>
978 {
979  semiparallel_only();
980 
981  PetscErrorCode ierr=0;
982 
983  ierr = MatAssemblyBegin (_mat, MAT_FLUSH_ASSEMBLY);
984  LIBMESH_CHKERR(ierr);
985  ierr = MatAssemblyEnd (_mat, MAT_FLUSH_ASSEMBLY);
986  LIBMESH_CHKERR(ierr);
987 }
988 
989 
990 
991 template <typename T>
993 {
994  libmesh_assert (this->initialized());
995 
996  PetscInt petsc_m=0, petsc_n=0;
997  PetscErrorCode ierr=0;
998 
999  ierr = MatGetSize (_mat, &petsc_m, &petsc_n);
1000  LIBMESH_CHKERR(ierr);
1001 
1002  return static_cast<numeric_index_type>(petsc_m);
1003 }
1004 
1005 
1006 
1007 template <typename T>
1009 {
1010  libmesh_assert (this->initialized());
1011 
1012  PetscInt petsc_m=0, petsc_n=0;
1013  PetscErrorCode ierr=0;
1014 
1015  ierr = MatGetSize (_mat, &petsc_m, &petsc_n);
1016  LIBMESH_CHKERR(ierr);
1017 
1018  return static_cast<numeric_index_type>(petsc_n);
1019 }
1020 
1021 
1022 
1023 template <typename T>
1025 {
1026  libmesh_assert (this->initialized());
1027 
1028  PetscInt start=0, stop=0;
1029  PetscErrorCode ierr=0;
1030 
1031  ierr = MatGetOwnershipRange(_mat, &start, &stop);
1032  LIBMESH_CHKERR(ierr);
1033 
1034  return static_cast<numeric_index_type>(start);
1035 }
1036 
1037 
1038 
1039 template <typename T>
1041 {
1042  libmesh_assert (this->initialized());
1043 
1044  PetscInt start=0, stop=0;
1045  PetscErrorCode ierr=0;
1046 
1047  ierr = MatGetOwnershipRange(_mat, &start, &stop);
1048  LIBMESH_CHKERR(ierr);
1049 
1050  return static_cast<numeric_index_type>(stop);
1051 }
1052 
1053 
1054 
1055 template <typename T>
1057  const numeric_index_type j,
1058  const T value)
1059 {
1060  libmesh_assert (this->initialized());
1061 
1062  PetscErrorCode ierr=0;
1063  PetscInt i_val=i, j_val=j;
1064 
1065  PetscScalar petsc_value = static_cast<PetscScalar>(value);
1066  ierr = MatSetValues(_mat, 1, &i_val, 1, &j_val,
1067  &petsc_value, INSERT_VALUES);
1068  LIBMESH_CHKERR(ierr);
1069 }
1070 
1071 
1072 
1073 template <typename T>
1075  const numeric_index_type j,
1076  const T value)
1077 {
1078  libmesh_assert (this->initialized());
1079 
1080  PetscErrorCode ierr=0;
1081  PetscInt i_val=i, j_val=j;
1082 
1083  PetscScalar petsc_value = static_cast<PetscScalar>(value);
1084  ierr = MatSetValues(_mat, 1, &i_val, 1, &j_val,
1085  &petsc_value, ADD_VALUES);
1086  LIBMESH_CHKERR(ierr);
1087 }
1088 
1089 
1090 
1091 template <typename T>
1093  const std::vector<numeric_index_type> & dof_indices)
1094 {
1095  this->add_matrix (dm, dof_indices, dof_indices);
1096 }
1097 
1098 
1099 
1100 
1101 
1102 
1103 
1104 template <typename T>
1105 void PetscMatrix<T>::add (const T a_in, const SparseMatrix<T> & X_in)
1106 {
1107  libmesh_assert (this->initialized());
1108 
1109  // sanity check. but this cannot avoid
1110  // crash due to incompatible sparsity structure...
1111  libmesh_assert_equal_to (this->m(), X_in.m());
1112  libmesh_assert_equal_to (this->n(), X_in.n());
1113 
1114  PetscScalar a = static_cast<PetscScalar> (a_in);
1115  const PetscMatrix<T> * X = cast_ptr<const PetscMatrix<T> *> (&X_in);
1116 
1117  libmesh_assert (X);
1118 
1119  PetscErrorCode ierr=0;
1120 
1121  // the matrix from which we copy the values has to be assembled/closed
1122  libmesh_assert(X->closed());
1123 
1124  semiparallel_only();
1125 
1126  ierr = MatAXPY(_mat, a, X->_mat, DIFFERENT_NONZERO_PATTERN);
1127  LIBMESH_CHKERR(ierr);
1128 }
1129 
1130 
1131 
1132 
1133 template <typename T>
1135  const numeric_index_type j_in) const
1136 {
1137  libmesh_assert (this->initialized());
1138 
1139  // PETSc 2.2.1 & newer
1140  const PetscScalar * petsc_row;
1141  const PetscInt * petsc_cols;
1142 
1143  // If the entry is not in the sparse matrix, it is 0.
1144  T value=0.;
1145 
1146  PetscErrorCode
1147  ierr=0;
1148  PetscInt
1149  ncols=0,
1150  i_val=static_cast<PetscInt>(i_in),
1151  j_val=static_cast<PetscInt>(j_in);
1152 
1153 
1154  // the matrix needs to be closed for this to work
1155  // this->close();
1156  // but closing it is a semiparallel operation; we want operator()
1157  // to run on one processor.
1158  libmesh_assert(this->closed());
1159 
1160  ierr = MatGetRow(_mat, i_val, &ncols, &petsc_cols, &petsc_row);
1161  LIBMESH_CHKERR(ierr);
1162 
1163  // Perform a binary search to find the contiguous index in
1164  // petsc_cols (resp. petsc_row) corresponding to global index j_val
1165  std::pair<const PetscInt *, const PetscInt *> p =
1166  std::equal_range (petsc_cols, petsc_cols + ncols, j_val);
1167 
1168  // Found an entry for j_val
1169  if (p.first != p.second)
1170  {
1171  // The entry in the contiguous row corresponding
1172  // to the j_val column of interest
1173  const std::size_t j =
1174  std::distance (const_cast<PetscInt *>(petsc_cols),
1175  const_cast<PetscInt *>(p.first));
1176 
1177  libmesh_assert_less (static_cast<PetscInt>(j), ncols);
1178  libmesh_assert_equal_to (petsc_cols[j], j_val);
1179 
1180  value = static_cast<T> (petsc_row[j]);
1181  }
1182 
1183  ierr = MatRestoreRow(_mat, i_val,
1184  &ncols, &petsc_cols, &petsc_row);
1185  LIBMESH_CHKERR(ierr);
1186 
1187  return value;
1188 }
1189 
1190 
1191 
1192 
1193 template <typename T>
1195 {
1196  libmesh_assert (this->initialized());
1197 
1198  PetscErrorCode ierr=0;
1199  PetscBool assembled;
1200 
1201  ierr = MatAssembled(_mat, &assembled);
1202  LIBMESH_CHKERR(ierr);
1203 
1204  return (assembled == PETSC_TRUE);
1205 }
1206 
1207 
1208 
1209 template <typename T>
1211 {
1212  std::swap(_mat, m_in._mat);
1213  std::swap(_destroy_mat_on_exit, m_in._destroy_mat_on_exit);
1214 }
1215 
1216 
1217 
1218 //------------------------------------------------------------------
1219 // Explicit instantiations
1220 template class PetscMatrix<Number>;
1221 
1222 } // namespace libMesh
1223 
1224 
1225 #endif // #ifdef LIBMESH_HAVE_PETSC
std::string name(const ElemQuality q)
Definition: elem_quality.C:42
virtual bool initialized() const
bool closed()
Definition: libmesh.C:265
virtual T operator()(const numeric_index_type i, const numeric_index_type j) const override
NumericVector interface to PETSc Vec.
Definition: petsc_vector.h:64
void set_matrix_type(PetscMatrixType mat_type)
Definition: petsc_matrix.C:121
virtual void print_matlab(const std::string &name="") const override
Definition: petsc_matrix.C:592
virtual void zero() override
Definition: petsc_matrix.C:475
PetscMatrix(const Parallel::Communicator &comm_in)
Definition: petsc_matrix.C:90
virtual numeric_index_type n() const override
Provides a uniform interface to vector storage schemes for different linear algebra libraries...
Definition: diff_context.h:40
void update_preallocation_and_zero()
Definition: petsc_matrix.C:468
virtual void init() override
Definition: petsc_matrix.C:348
virtual void add_block_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &brows, const std::vector< numeric_index_type > &bcols) override
Definition: petsc_matrix.C:805
unsigned int m() const
virtual void print_personal(std::ostream &os=libMesh::out) const override
Definition: petsc_matrix.C:671
const Number zero
Definition: libmesh.h:239
void swap(PetscMatrix< T > &)
virtual void clear()=0
PetscInt * numeric_petsc_cast(const numeric_index_type *p)
void libmesh_ignore(const Args &...)
dof_id_type numeric_index_type
Definition: id_types.h:92
virtual numeric_index_type m() const =0
void stop(const char *file, int line, const char *date, const char *time)
virtual numeric_index_type row_stop() const override
std::vector< T > & get_values()
Definition: dense_matrix.h:341
virtual numeric_index_type m() const override
Definition: petsc_matrix.C:992
virtual void add(const numeric_index_type i, const numeric_index_type j, const T value) override
virtual void get_diagonal(NumericVector< T > &dest) const override
Definition: petsc_matrix.C:908
PetscTruth PetscBool
Definition: petsc_macro.h:67
virtual void clear() override
Definition: petsc_matrix.C:528
DIE A HORRIBLE DEATH HERE typedef LIBMESH_DEFAULT_SCALAR_TYPE Real
virtual bool closed() const override
virtual numeric_index_type row_start() const override
virtual void get_transpose(SparseMatrix< T > &dest) const override
Definition: petsc_matrix.C:921
PetscErrorCode ierr
virtual Real l1_norm() const override
Definition: petsc_matrix.C:546
void swap(Iterator &lhs, Iterator &rhs)
static const bool value
Definition: xdr_io.C:109
virtual void add_matrix(const DenseMatrix< T > &dm, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols) override
Definition: petsc_matrix.C:778
virtual void set(const numeric_index_type i, const numeric_index_type j, const T value) override
SparseMatrix interface to PETSc Mat.
bool initialized()
Definition: libmesh.C:258
virtual void _get_submatrix(SparseMatrix< T > &submatrix, const std::vector< numeric_index_type > &rows, const std::vector< numeric_index_type > &cols, const bool reuse_submatrix) const override
Definition: petsc_matrix.C:849
virtual void zero_rows(std::vector< numeric_index_type > &rows, T diag_value=0.0) override
Definition: petsc_matrix.C:496
unsigned int n() const
A matrix object used for finite element assembly and numerics.
Definition: dense_matrix.h:54
virtual void flush() override
Definition: petsc_matrix.C:977
virtual void close() override
Definition: petsc_matrix.C:957
virtual Real linfty_norm() const override
Definition: petsc_matrix.C:569
virtual numeric_index_type n() const =0