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