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

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

Hosted By:
SourceForge.net Logo