1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2020 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 #include "libmesh/parallel.h"
33 #include "libmesh/utility.h"
34
35
36 #ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
37
38 namespace
39 {
40 using namespace libMesh;
41
42 // historic libMesh n_nz & n_oz arrays are set up for PETSc's AIJ format.
43 // however, when the blocksize is >1, we need to transform these into
44 // their BAIJ counterparts.
45 inline
transform_preallocation_arrays(const PetscInt blocksize,const std::vector<numeric_index_type> & n_nz,const std::vector<numeric_index_type> & n_oz,std::vector<numeric_index_type> & b_n_nz,std::vector<numeric_index_type> & b_n_oz)46 void transform_preallocation_arrays (const PetscInt blocksize,
47 const std::vector<numeric_index_type> & n_nz,
48 const std::vector<numeric_index_type> & n_oz,
49 std::vector<numeric_index_type> & b_n_nz,
50 std::vector<numeric_index_type> & b_n_oz)
51 {
52 libmesh_assert_equal_to (n_nz.size(), n_oz.size());
53 libmesh_assert_equal_to (n_nz.size()%blocksize, 0);
54
55 b_n_nz.clear(); /**/ b_n_nz.reserve(n_nz.size()/blocksize);
56 b_n_oz.clear(); /**/ b_n_oz.reserve(n_oz.size()/blocksize);
57
58 for (std::size_t nn=0, nnzs=n_nz.size(); nn<nnzs; nn += blocksize)
59 {
60 b_n_nz.push_back (n_nz[nn]/blocksize);
61 b_n_oz.push_back (n_oz[nn]/blocksize);
62 }
63 }
64 }
65
66 #endif
67
68
69
70 namespace libMesh
71 {
72
73
74 //-----------------------------------------------------------------------
75 // PetscMatrix members
76
77
78 // Constructor
79 template <typename T>
PetscMatrix(const Parallel::Communicator & comm_in)80 PetscMatrix<T>::PetscMatrix(const Parallel::Communicator & comm_in) :
81 SparseMatrix<T>(comm_in),
82 _destroy_mat_on_exit(true),
83 _mat_type(AIJ)
84 {}
85
86
87
88 // Constructor taking an existing Mat but not the responsibility
89 // for destroying it
90 template <typename T>
PetscMatrix(Mat mat_in,const Parallel::Communicator & comm_in)91 PetscMatrix<T>::PetscMatrix(Mat mat_in,
92 const Parallel::Communicator & comm_in) :
93 SparseMatrix<T>(comm_in),
94 _destroy_mat_on_exit(false),
95 _mat_type(AIJ)
96 {
97 this->_mat = mat_in;
98 this->_is_initialized = true;
99 }
100
101
102
103 // Destructor
104 template <typename T>
~PetscMatrix()105 PetscMatrix<T>::~PetscMatrix()
106 {
107 this->clear();
108 }
109
110 template <typename T>
set_matrix_type(PetscMatrixType mat_type)111 void PetscMatrix<T>::set_matrix_type(PetscMatrixType mat_type)
112 {
113 _mat_type = mat_type;
114 }
115
116 template <typename T>
init(const numeric_index_type m_in,const numeric_index_type n_in,const numeric_index_type m_l,const numeric_index_type n_l,const numeric_index_type nnz,const numeric_index_type noz,const numeric_index_type blocksize_in)117 void PetscMatrix<T>::init (const numeric_index_type m_in,
118 const numeric_index_type n_in,
119 const numeric_index_type m_l,
120 const numeric_index_type n_l,
121 const numeric_index_type nnz,
122 const numeric_index_type noz,
123 const numeric_index_type blocksize_in)
124 {
125 // So compilers don't warn when !LIBMESH_ENABLE_BLOCKED_STORAGE
126 libmesh_ignore(blocksize_in);
127
128 // Clear initialized matrices
129 if (this->initialized())
130 this->clear();
131
132 this->_is_initialized = true;
133
134
135 PetscErrorCode ierr = 0;
136 PetscInt m_global = static_cast<PetscInt>(m_in);
137 PetscInt n_global = static_cast<PetscInt>(n_in);
138 PetscInt m_local = static_cast<PetscInt>(m_l);
139 PetscInt n_local = static_cast<PetscInt>(n_l);
140 PetscInt n_nz = static_cast<PetscInt>(nnz);
141 PetscInt n_oz = static_cast<PetscInt>(noz);
142
143 ierr = MatCreate(this->comm().get(), &_mat);
144 LIBMESH_CHKERR(ierr);
145 ierr = MatSetSizes(_mat, m_local, n_local, m_global, n_global);
146 LIBMESH_CHKERR(ierr);
147 PetscInt blocksize = static_cast<PetscInt>(blocksize_in);
148 ierr = MatSetBlockSize(_mat,blocksize);
149 LIBMESH_CHKERR(ierr);
150
151 #ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
152 if (blocksize > 1)
153 {
154 // specified blocksize, bs>1.
155 // double check sizes.
156 libmesh_assert_equal_to (m_local % blocksize, 0);
157 libmesh_assert_equal_to (n_local % blocksize, 0);
158 libmesh_assert_equal_to (m_global % blocksize, 0);
159 libmesh_assert_equal_to (n_global % blocksize, 0);
160 libmesh_assert_equal_to (n_nz % blocksize, 0);
161 libmesh_assert_equal_to (n_oz % blocksize, 0);
162
163 ierr = MatSetType(_mat, MATBAIJ); // Automatically chooses seqbaij or mpibaij
164 LIBMESH_CHKERR(ierr);
165 ierr = MatSeqBAIJSetPreallocation(_mat, blocksize, n_nz/blocksize, NULL);
166 LIBMESH_CHKERR(ierr);
167 ierr = MatMPIBAIJSetPreallocation(_mat, blocksize,
168 n_nz/blocksize, NULL,
169 n_oz/blocksize, NULL);
170 LIBMESH_CHKERR(ierr);
171 }
172 else
173 #endif
174 {
175 switch (_mat_type) {
176 case AIJ:
177 ierr = MatSetType(_mat, MATAIJ); // Automatically chooses seqaij or mpiaij
178 LIBMESH_CHKERR(ierr);
179 ierr = MatSeqAIJSetPreallocation(_mat, n_nz, NULL);
180 LIBMESH_CHKERR(ierr);
181 ierr = MatMPIAIJSetPreallocation(_mat, n_nz, NULL, n_oz, NULL);
182 LIBMESH_CHKERR(ierr);
183 break;
184
185 case HYPRE:
186 #if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
187 ierr = MatSetType(_mat, MATHYPRE);
188 LIBMESH_CHKERR(ierr);
189 ierr = MatHYPRESetPreallocation(_mat, n_nz, NULL, n_oz, NULL);
190 LIBMESH_CHKERR(ierr);
191 #else
192 libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
193 #endif
194 break;
195
196 default: libmesh_error_msg("Unsupported petsc matrix type");
197 }
198 }
199
200 // Make it an error for PETSc to allocate new nonzero entries during assembly
201 ierr = MatSetOption(_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
202 LIBMESH_CHKERR(ierr);
203
204 // Is prefix information available somewhere? Perhaps pass in the system name?
205 ierr = MatSetOptionsPrefix(_mat, "");
206 LIBMESH_CHKERR(ierr);
207 ierr = MatSetFromOptions(_mat);
208 LIBMESH_CHKERR(ierr);
209
210 this->zero ();
211 }
212
213
214 template <typename T>
init(const numeric_index_type m_in,const numeric_index_type n_in,const numeric_index_type m_l,const numeric_index_type n_l,const std::vector<numeric_index_type> & n_nz,const std::vector<numeric_index_type> & n_oz,const numeric_index_type blocksize_in)215 void PetscMatrix<T>::init (const numeric_index_type m_in,
216 const numeric_index_type n_in,
217 const numeric_index_type m_l,
218 const numeric_index_type n_l,
219 const std::vector<numeric_index_type> & n_nz,
220 const std::vector<numeric_index_type> & n_oz,
221 const numeric_index_type blocksize_in)
222 {
223 // So compilers don't warn when !LIBMESH_ENABLE_BLOCKED_STORAGE
224 libmesh_ignore(blocksize_in);
225
226 // Clear initialized matrices
227 if (this->initialized())
228 this->clear();
229
230 this->_is_initialized = true;
231
232 // Make sure the sparsity pattern isn't empty unless the matrix is 0x0
233 libmesh_assert_equal_to (n_nz.size(), m_l);
234 libmesh_assert_equal_to (n_oz.size(), m_l);
235
236 PetscErrorCode ierr = 0;
237 PetscInt m_global = static_cast<PetscInt>(m_in);
238 PetscInt n_global = static_cast<PetscInt>(n_in);
239 PetscInt m_local = static_cast<PetscInt>(m_l);
240 PetscInt n_local = static_cast<PetscInt>(n_l);
241
242 ierr = MatCreate(this->comm().get(), &_mat);
243 LIBMESH_CHKERR(ierr);
244 ierr = MatSetSizes(_mat, m_local, n_local, m_global, n_global);
245 LIBMESH_CHKERR(ierr);
246 PetscInt blocksize = static_cast<PetscInt>(blocksize_in);
247 ierr = MatSetBlockSize(_mat,blocksize);
248 LIBMESH_CHKERR(ierr);
249
250 #ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
251 if (blocksize > 1)
252 {
253 // specified blocksize, bs>1.
254 // double check sizes.
255 libmesh_assert_equal_to (m_local % blocksize, 0);
256 libmesh_assert_equal_to (n_local % blocksize, 0);
257 libmesh_assert_equal_to (m_global % blocksize, 0);
258 libmesh_assert_equal_to (n_global % blocksize, 0);
259
260 ierr = MatSetType(_mat, MATBAIJ); // Automatically chooses seqbaij or mpibaij
261 LIBMESH_CHKERR(ierr);
262
263 // transform the per-entry n_nz and n_oz arrays into their block counterparts.
264 std::vector<numeric_index_type> b_n_nz, b_n_oz;
265
266 transform_preallocation_arrays (blocksize,
267 n_nz, n_oz,
268 b_n_nz, b_n_oz);
269
270 ierr = MatSeqBAIJSetPreallocation (_mat,
271 blocksize,
272 0,
273 numeric_petsc_cast(b_n_nz.empty() ? nullptr : b_n_nz.data()));
274 LIBMESH_CHKERR(ierr);
275
276 ierr = MatMPIBAIJSetPreallocation (_mat,
277 blocksize,
278 0,
279 numeric_petsc_cast(b_n_nz.empty() ? nullptr : b_n_nz.data()),
280 0,
281 numeric_petsc_cast(b_n_oz.empty() ? nullptr : b_n_oz.data()));
282 LIBMESH_CHKERR(ierr);
283 }
284 else
285 #endif
286 {
287 switch (_mat_type) {
288 case AIJ:
289 ierr = MatSetType(_mat, MATAIJ); // Automatically chooses seqaij or mpiaij
290 LIBMESH_CHKERR(ierr);
291 ierr = MatSeqAIJSetPreallocation (_mat,
292 0,
293 numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()));
294 LIBMESH_CHKERR(ierr);
295 ierr = MatMPIAIJSetPreallocation (_mat,
296 0,
297 numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
298 0,
299 numeric_petsc_cast(n_oz.empty() ? nullptr : n_oz.data()));
300 break;
301
302 case HYPRE:
303 #if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
304 ierr = MatSetType(_mat, MATHYPRE);
305 LIBMESH_CHKERR(ierr);
306 ierr = MatHYPRESetPreallocation (_mat,
307 0,
308 numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
309 0,
310 numeric_petsc_cast(n_oz.empty() ? nullptr : n_oz.data()));
311 LIBMESH_CHKERR(ierr);
312 #else
313 libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
314 #endif
315 break;
316
317 default: libmesh_error_msg("Unsupported petsc matrix type");
318 }
319
320 }
321
322 // Make it an error for PETSc to allocate new nonzero entries during assembly
323 ierr = MatSetOption(_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
324 LIBMESH_CHKERR(ierr);
325
326 // Is prefix information available somewhere? Perhaps pass in the system name?
327 ierr = MatSetOptionsPrefix(_mat, "");
328 LIBMESH_CHKERR(ierr);
329 ierr = MatSetFromOptions(_mat);
330 LIBMESH_CHKERR(ierr);
331
332
333 this->zero();
334 }
335
336
337 template <typename T>
init(const ParallelType)338 void PetscMatrix<T>::init (const ParallelType)
339 {
340 libmesh_assert(this->_dof_map);
341
342 // Clear initialized matrices
343 if (this->initialized())
344 this->clear();
345
346 this->_is_initialized = true;
347
348
349 const numeric_index_type my_m = this->_dof_map->n_dofs();
350 const numeric_index_type my_n = my_m;
351 const numeric_index_type n_l = this->_dof_map->n_dofs_on_processor(this->processor_id());
352 const numeric_index_type m_l = n_l;
353
354
355 const std::vector<numeric_index_type> & n_nz = this->_dof_map->get_n_nz();
356 const std::vector<numeric_index_type> & n_oz = this->_dof_map->get_n_oz();
357
358 // Make sure the sparsity pattern isn't empty unless the matrix is 0x0
359 libmesh_assert_equal_to (n_nz.size(), m_l);
360 libmesh_assert_equal_to (n_oz.size(), m_l);
361
362 PetscErrorCode ierr = 0;
363 PetscInt m_global = static_cast<PetscInt>(my_m);
364 PetscInt n_global = static_cast<PetscInt>(my_n);
365 PetscInt m_local = static_cast<PetscInt>(m_l);
366 PetscInt n_local = static_cast<PetscInt>(n_l);
367
368 ierr = MatCreate(this->comm().get(), &_mat);
369 LIBMESH_CHKERR(ierr);
370 ierr = MatSetSizes(_mat, m_local, n_local, m_global, n_global);
371 LIBMESH_CHKERR(ierr);
372 PetscInt blocksize = static_cast<PetscInt>(this->_dof_map->block_size());
373 ierr = MatSetBlockSize(_mat,blocksize);
374 LIBMESH_CHKERR(ierr);
375
376 #ifdef LIBMESH_ENABLE_BLOCKED_STORAGE
377 if (blocksize > 1)
378 {
379 // specified blocksize, bs>1.
380 // double check sizes.
381 libmesh_assert_equal_to (m_local % blocksize, 0);
382 libmesh_assert_equal_to (n_local % blocksize, 0);
383 libmesh_assert_equal_to (m_global % blocksize, 0);
384 libmesh_assert_equal_to (n_global % blocksize, 0);
385
386 ierr = MatSetType(_mat, MATBAIJ); // Automatically chooses seqbaij or mpibaij
387 LIBMESH_CHKERR(ierr);
388
389 // transform the per-entry n_nz and n_oz arrays into their block counterparts.
390 std::vector<numeric_index_type> b_n_nz, b_n_oz;
391
392 transform_preallocation_arrays (blocksize,
393 n_nz, n_oz,
394 b_n_nz, b_n_oz);
395
396 ierr = MatSeqBAIJSetPreallocation (_mat,
397 blocksize,
398 0,
399 numeric_petsc_cast(b_n_nz.empty() ? nullptr : b_n_nz.data()));
400 LIBMESH_CHKERR(ierr);
401
402 ierr = MatMPIBAIJSetPreallocation (_mat,
403 blocksize,
404 0,
405 numeric_petsc_cast(b_n_nz.empty() ? nullptr : b_n_nz.data()),
406 0,
407 numeric_petsc_cast(b_n_oz.empty() ? nullptr : b_n_oz.data()));
408 LIBMESH_CHKERR(ierr);
409 }
410 else
411 #endif
412 {
413 switch (_mat_type) {
414 case AIJ:
415 ierr = MatSetType(_mat, MATAIJ); // Automatically chooses seqaij or mpiaij
416 LIBMESH_CHKERR(ierr);
417 ierr = MatSeqAIJSetPreallocation (_mat,
418 0,
419 numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()));
420 LIBMESH_CHKERR(ierr);
421 ierr = MatMPIAIJSetPreallocation (_mat,
422 0,
423 numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
424 0,
425 numeric_petsc_cast(n_oz.empty() ? nullptr : n_oz.data()));
426 break;
427
428 case HYPRE:
429 #if !PETSC_VERSION_LESS_THAN(3,9,4) && LIBMESH_HAVE_PETSC_HYPRE
430 ierr = MatSetType(_mat, MATHYPRE);
431 LIBMESH_CHKERR(ierr);
432 ierr = MatHYPRESetPreallocation (_mat,
433 0,
434 numeric_petsc_cast(n_nz.empty() ? nullptr : n_nz.data()),
435 0,
436 numeric_petsc_cast(n_oz.empty() ? nullptr : n_oz.data()));
437 LIBMESH_CHKERR(ierr);
438 #else
439 libmesh_error_msg("PETSc 3.9.4 or higher with hypre is required for MatHypre");
440 #endif
441 break;
442
443 default: libmesh_error_msg("Unsupported petsc matrix type");
444 }
445 }
446
447 // Make it an error for PETSc to allocate new nonzero entries during assembly
448 ierr = MatSetOption(_mat, MAT_NEW_NONZERO_ALLOCATION_ERR, PETSC_TRUE);
449 LIBMESH_CHKERR(ierr);
450
451 // Is prefix information available somewhere? Perhaps pass in the system name?
452 ierr = MatSetOptionsPrefix(_mat, "");
453 LIBMESH_CHKERR(ierr);
454 ierr = MatSetFromOptions(_mat);
455 LIBMESH_CHKERR(ierr);
456
457 this->zero();
458 }
459
460
461 template <typename T>
update_preallocation_and_zero()462 void PetscMatrix<T>::update_preallocation_and_zero ()
463 {
464 libmesh_not_implemented();
465 }
466
467 template <typename T>
reset_preallocation()468 void PetscMatrix<T>::reset_preallocation()
469 {
470 #if !PETSC_VERSION_LESS_THAN(3,9,0)
471 libmesh_assert (this->initialized());
472
473 auto ierr = MatResetPreallocation(_mat);
474 LIBMESH_CHKERR(ierr);
475 #else
476 libmesh_warning("Your version of PETSc doesn't support resetting of "
477 "preallocation, so we will use your most recent sparsity "
478 "pattern. This may result in a degradation of performance\n");
479 #endif
480 }
481
482 template <typename T>
zero()483 void PetscMatrix<T>::zero ()
484 {
485 libmesh_assert (this->initialized());
486
487 semiparallel_only();
488
489 PetscErrorCode ierr=0;
490
491 PetscInt m_l, n_l;
492
493 ierr = MatGetLocalSize(_mat,&m_l,&n_l);
494 LIBMESH_CHKERR(ierr);
495
496 if (n_l)
497 {
498 ierr = MatZeroEntries(_mat);
499 LIBMESH_CHKERR(ierr);
500 }
501 }
502
503 template <typename T>
zero_rows(std::vector<numeric_index_type> & rows,T diag_value)504 void PetscMatrix<T>::zero_rows (std::vector<numeric_index_type> & rows, T diag_value)
505 {
506 libmesh_assert (this->initialized());
507
508 semiparallel_only();
509
510 PetscErrorCode ierr=0;
511
512 // As of petsc-dev at the time of 3.1.0, MatZeroRows now takes two additional
513 // optional arguments. The optional arguments (x,b) can be used to specify the
514 // solutions for the zeroed rows (x) and right hand side (b) to update.
515 // Could be useful for setting boundary conditions...
516 if (!rows.empty())
517 ierr = MatZeroRows(_mat, cast_int<PetscInt>(rows.size()),
518 numeric_petsc_cast(rows.data()), PS(diag_value),
519 NULL, NULL);
520 else
521 ierr = MatZeroRows(_mat, 0, NULL, PS(diag_value), NULL, NULL);
522
523 LIBMESH_CHKERR(ierr);
524 }
525
526 template <typename T>
zero_clone()527 std::unique_ptr<SparseMatrix<T>> PetscMatrix<T>::zero_clone () const
528 {
529 libmesh_error_msg_if(!this->closed(), "Matrix must be closed before it can be cloned!");
530
531 // Copy the nonzero pattern only
532 Mat copy;
533 PetscErrorCode ierr = MatDuplicate(_mat, MAT_DO_NOT_COPY_VALUES, ©);
534 LIBMESH_CHKERR(ierr);
535
536 // Call wrapping PetscMatrix constructor, have it take over
537 // ownership.
538 auto ret = libmesh_make_unique<PetscMatrix<T>>(copy, this->comm());
539 ret->set_destroy_mat_on_exit(true);
540
541 // Work around an issue on older compilers. We are able to simply
542 // "return ret;" on newer compilers
543 return std::unique_ptr<SparseMatrix<T>>(ret.release());
544 }
545
546
547
548 template <typename T>
clone()549 std::unique_ptr<SparseMatrix<T>> PetscMatrix<T>::clone () const
550 {
551 libmesh_error_msg_if(!this->closed(), "Matrix must be closed before it can be cloned!");
552
553 // Copy the nonzero pattern and numerical values
554 Mat copy;
555 PetscErrorCode ierr = MatDuplicate(_mat, MAT_COPY_VALUES, ©);
556 LIBMESH_CHKERR(ierr);
557
558 // Call wrapping PetscMatrix constructor, have it take over
559 // ownership.
560 auto ret = libmesh_make_unique<PetscMatrix<T>>(copy, this->comm());
561 ret->set_destroy_mat_on_exit(true);
562
563 // Work around an issue on older compilers. We are able to simply
564 // "return ret;" on newer compilers
565 return std::unique_ptr<SparseMatrix<T>>(ret.release());
566 }
567
568 template <typename T>
clear()569 void PetscMatrix<T>::clear ()
570 {
571 PetscErrorCode ierr=0;
572
573 if ((this->initialized()) && (this->_destroy_mat_on_exit))
574 {
575 semiparallel_only();
576
577 ierr = MatDestroy (&_mat);
578 LIBMESH_CHKERR(ierr);
579
580 this->_is_initialized = false;
581 }
582 }
583
584
585
586 template <typename T>
l1_norm()587 Real PetscMatrix<T>::l1_norm () const
588 {
589 libmesh_assert (this->initialized());
590
591 semiparallel_only();
592
593 PetscErrorCode ierr=0;
594 PetscReal petsc_value;
595 Real value;
596
597 libmesh_assert (this->closed());
598
599 ierr = MatNorm(_mat, NORM_1, &petsc_value);
600 LIBMESH_CHKERR(ierr);
601
602 value = static_cast<Real>(petsc_value);
603
604 return value;
605 }
606
607
608
609 template <typename T>
linfty_norm()610 Real PetscMatrix<T>::linfty_norm () const
611 {
612 libmesh_assert (this->initialized());
613
614 semiparallel_only();
615
616 PetscErrorCode ierr=0;
617 PetscReal petsc_value;
618 Real value;
619
620 libmesh_assert (this->closed());
621
622 ierr = MatNorm(_mat, NORM_INFINITY, &petsc_value);
623 LIBMESH_CHKERR(ierr);
624
625 value = static_cast<Real>(petsc_value);
626
627 return value;
628 }
629
630
631
632 template <typename T>
print_matlab(const std::string & name)633 void PetscMatrix<T>::print_matlab (const std::string & name) const
634 {
635 libmesh_assert (this->initialized());
636
637 semiparallel_only();
638
639 if (!this->closed())
640 {
641 libmesh_deprecated();
642 libmesh_warning("The matrix must be assembled before calling PetscMatrix::print_matlab().\n"
643 "Please update your code, as this warning will become an error in a future release.");
644 const_cast<PetscMatrix<T> *>(this)->close();
645 }
646
647 PetscErrorCode ierr=0;
648 PetscViewer petsc_viewer;
649
650
651 ierr = PetscViewerCreate (this->comm().get(),
652 &petsc_viewer);
653 LIBMESH_CHKERR(ierr);
654
655 /**
656 * Create an ASCII file containing the matrix
657 * if a filename was provided.
658 */
659 if (name != "")
660 {
661 ierr = PetscViewerASCIIOpen( this->comm().get(),
662 name.c_str(),
663 &petsc_viewer);
664 LIBMESH_CHKERR(ierr);
665
666 #if PETSC_VERSION_LESS_THAN(3,7,0)
667 ierr = PetscViewerSetFormat (petsc_viewer,
668 PETSC_VIEWER_ASCII_MATLAB);
669 #else
670 ierr = PetscViewerPushFormat (petsc_viewer,
671 PETSC_VIEWER_ASCII_MATLAB);
672 #endif
673
674 LIBMESH_CHKERR(ierr);
675
676 ierr = MatView (_mat, petsc_viewer);
677 LIBMESH_CHKERR(ierr);
678 }
679
680 /**
681 * Otherwise the matrix will be dumped to the screen.
682 */
683 else
684 {
685 #if PETSC_VERSION_LESS_THAN(3,7,0)
686 ierr = PetscViewerSetFormat (PETSC_VIEWER_STDOUT_WORLD,
687 PETSC_VIEWER_ASCII_MATLAB);
688 #else
689 ierr = PetscViewerPushFormat (PETSC_VIEWER_STDOUT_WORLD,
690 PETSC_VIEWER_ASCII_MATLAB);
691 #endif
692
693 LIBMESH_CHKERR(ierr);
694
695 ierr = MatView (_mat, PETSC_VIEWER_STDOUT_WORLD);
696 LIBMESH_CHKERR(ierr);
697 }
698
699
700 /**
701 * Destroy the viewer.
702 */
703 ierr = PetscViewerDestroy (&petsc_viewer);
704 LIBMESH_CHKERR(ierr);
705 }
706
707
708
709
710
711 template <typename T>
print_personal(std::ostream & os)712 void PetscMatrix<T>::print_personal(std::ostream & os) const
713 {
714 libmesh_assert (this->initialized());
715
716 // Routine must be called in parallel on parallel matrices
717 // and serial on serial matrices.
718 semiparallel_only();
719
720 // #ifndef NDEBUG
721 // if (os != std::cout)
722 // libMesh::err << "Warning! PETSc can only print to std::cout!" << std::endl;
723 // #endif
724
725 // Matrix must be in an assembled state to be printed
726 if (!this->closed())
727 {
728 libmesh_deprecated();
729 libmesh_warning("The matrix must be assembled before calling PetscMatrix::print_personal().\n"
730 "Please update your code, as this warning will become an error in a future release.");
731 const_cast<PetscMatrix<T> *>(this)->close();
732 }
733
734 PetscErrorCode ierr=0;
735
736 // Print to screen if ostream is stdout
737 if (os.rdbuf() == std::cout.rdbuf())
738 {
739 ierr = MatView(_mat, PETSC_VIEWER_STDOUT_SELF);
740 LIBMESH_CHKERR(ierr);
741 }
742
743 // Otherwise, print to the requested file, in a roundabout way...
744 else
745 {
746 // We will create a temporary filename, and file, for PETSc to
747 // write to.
748 std::string temp_filename;
749
750 {
751 // Template for temporary filename
752 char c[] = "temp_petsc_matrix.XXXXXX";
753
754 // Generate temporary, unique filename only on processor 0. We will
755 // use this filename for PetscViewerASCIIOpen, before copying it into
756 // the user's stream
757 if (this->processor_id() == 0)
758 {
759 int fd = mkstemp(c);
760
761 // Check to see that mkstemp did not fail.
762 libmesh_error_msg_if(fd == -1, "mkstemp failed in PetscMatrix::print_personal()");
763
764 // mkstemp returns a file descriptor for an open file,
765 // so let's close it before we hand it to PETSc!
766 ::close (fd);
767 }
768
769 // Store temporary filename as string, makes it easier to broadcast
770 temp_filename = c;
771 }
772
773 // Now broadcast the filename from processor 0 to all processors.
774 this->comm().broadcast(temp_filename);
775
776 // PetscViewer object for passing to MatView
777 PetscViewer petsc_viewer;
778
779 // This PETSc function only takes a string and handles the opening/closing
780 // of the file internally. Since print_personal() takes a reference to
781 // an ostream, we have to do an extra step... print_personal() should probably
782 // have a version that takes a string to get rid of this problem.
783 ierr = PetscViewerASCIIOpen( this->comm().get(),
784 temp_filename.c_str(),
785 &petsc_viewer);
786 LIBMESH_CHKERR(ierr);
787
788 // Probably don't need to set the format if it's default...
789 // ierr = PetscViewerSetFormat (petsc_viewer,
790 // PETSC_VIEWER_DEFAULT);
791 // LIBMESH_CHKERR(ierr);
792
793 // Finally print the matrix using the viewer
794 ierr = MatView (_mat, petsc_viewer);
795 LIBMESH_CHKERR(ierr);
796
797 if (this->processor_id() == 0)
798 {
799 // Now the inefficient bit: open temp_filename as an ostream and copy the contents
800 // into the user's desired ostream. We can't just do a direct file copy, we don't even have the filename!
801 std::ifstream input_stream(temp_filename.c_str());
802 os << input_stream.rdbuf(); // The "most elegant" way to copy one stream into another.
803 // os.close(); // close not defined in ostream
804
805 // Now remove the temporary file
806 input_stream.close();
807 std::remove(temp_filename.c_str());
808 }
809 }
810 }
811
812
813
814
815
816
817 template <typename T>
add_matrix(const DenseMatrix<T> & dm,const std::vector<numeric_index_type> & rows,const std::vector<numeric_index_type> & cols)818 void PetscMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
819 const std::vector<numeric_index_type> & rows,
820 const std::vector<numeric_index_type> & cols)
821 {
822 libmesh_assert (this->initialized());
823
824 const numeric_index_type n_rows = dm.m();
825 const numeric_index_type n_cols = dm.n();
826
827 libmesh_assert_equal_to (rows.size(), n_rows);
828 libmesh_assert_equal_to (cols.size(), n_cols);
829
830 PetscErrorCode ierr=0;
831 ierr = MatSetValues(_mat,
832 n_rows, numeric_petsc_cast(rows.data()),
833 n_cols, numeric_petsc_cast(cols.data()),
834 pPS(const_cast<T*>(dm.get_values().data())),
835 ADD_VALUES);
836 LIBMESH_CHKERR(ierr);
837 }
838
839
840
841
842
843
844 template <typename T>
add_block_matrix(const DenseMatrix<T> & dm,const std::vector<numeric_index_type> & brows,const std::vector<numeric_index_type> & bcols)845 void PetscMatrix<T>::add_block_matrix(const DenseMatrix<T> & dm,
846 const std::vector<numeric_index_type> & brows,
847 const std::vector<numeric_index_type> & bcols)
848 {
849 libmesh_assert (this->initialized());
850
851 const numeric_index_type n_brows =
852 cast_int<numeric_index_type>(brows.size());
853 const numeric_index_type n_bcols =
854 cast_int<numeric_index_type>(bcols.size());
855
856 PetscErrorCode ierr=0;
857
858 #ifndef NDEBUG
859 const numeric_index_type n_rows =
860 cast_int<numeric_index_type>(dm.m());
861 const numeric_index_type n_cols =
862 cast_int<numeric_index_type>(dm.n());
863 const numeric_index_type blocksize = n_rows / n_brows;
864
865 libmesh_assert_equal_to (n_cols / n_bcols, blocksize);
866 libmesh_assert_equal_to (blocksize*n_brows, n_rows);
867 libmesh_assert_equal_to (blocksize*n_bcols, n_cols);
868
869 PetscInt petsc_blocksize;
870 ierr = MatGetBlockSize(_mat, &petsc_blocksize);
871 LIBMESH_CHKERR(ierr);
872 libmesh_assert_equal_to (blocksize, static_cast<numeric_index_type>(petsc_blocksize));
873 #endif
874
875 // These casts are required for PETSc <= 2.1.5
876 ierr = MatSetValuesBlocked(_mat,
877 n_brows, numeric_petsc_cast(brows.data()),
878 n_bcols, numeric_petsc_cast(bcols.data()),
879 pPS(const_cast<T*>(dm.get_values().data())),
880 ADD_VALUES);
881 LIBMESH_CHKERR(ierr);
882 }
883
884
885
886
887
888 template <typename T>
_get_submatrix(SparseMatrix<T> & submatrix,const std::vector<numeric_index_type> & rows,const std::vector<numeric_index_type> & cols,const bool reuse_submatrix)889 void PetscMatrix<T>::_get_submatrix(SparseMatrix<T> & submatrix,
890 const std::vector<numeric_index_type> & rows,
891 const std::vector<numeric_index_type> & cols,
892 const bool reuse_submatrix) const
893 {
894 if (!this->closed())
895 {
896 libmesh_deprecated();
897 libmesh_warning("The matrix must be assembled before calling PetscMatrix::create_submatrix().\n"
898 "Please update your code, as this warning will become an error in a future release.");
899 const_cast<PetscMatrix<T> *>(this)->close();
900 }
901
902 // Make sure the SparseMatrix passed in is really a PetscMatrix
903 PetscMatrix<T> * petsc_submatrix = cast_ptr<PetscMatrix<T> *>(&submatrix);
904
905 // If we're not reusing submatrix and submatrix is already initialized
906 // then we need to clear it, otherwise we get a memory leak.
907 if (!reuse_submatrix && submatrix.initialized())
908 submatrix.clear();
909
910 // Construct row and column index sets.
911 PetscErrorCode ierr=0;
912 IS isrow, iscol;
913
914 ierr = ISCreateGeneral(this->comm().get(),
915 cast_int<PetscInt>(rows.size()),
916 numeric_petsc_cast(rows.data()),
917 PETSC_USE_POINTER,
918 &isrow); LIBMESH_CHKERR(ierr);
919
920 ierr = ISCreateGeneral(this->comm().get(),
921 cast_int<PetscInt>(cols.size()),
922 numeric_petsc_cast(cols.data()),
923 PETSC_USE_POINTER,
924 &iscol); LIBMESH_CHKERR(ierr);
925
926 // Extract submatrix
927 ierr = LibMeshCreateSubMatrix(_mat,
928 isrow,
929 iscol,
930 (reuse_submatrix ? MAT_REUSE_MATRIX : MAT_INITIAL_MATRIX),
931 &(petsc_submatrix->_mat)); LIBMESH_CHKERR(ierr);
932
933 // Specify that the new submatrix is initialized and close it.
934 petsc_submatrix->_is_initialized = true;
935 petsc_submatrix->close();
936
937 // Clean up PETSc data structures
938 ierr = ISDestroy(&isrow); LIBMESH_CHKERR(ierr);
939 ierr = ISDestroy(&iscol); LIBMESH_CHKERR(ierr);
940 }
941
942
943
944 template <typename T>
get_diagonal(NumericVector<T> & dest)945 void PetscMatrix<T>::get_diagonal (NumericVector<T> & dest) const
946 {
947 // Make sure the NumericVector passed in is really a PetscVector
948 PetscVector<T> & petsc_dest = cast_ref<PetscVector<T> &>(dest);
949
950 // Needs a const_cast since PETSc does not work with const.
951 PetscErrorCode ierr =
952 MatGetDiagonal(const_cast<PetscMatrix<T> *>(this)->mat(),petsc_dest.vec()); LIBMESH_CHKERR(ierr);
953 }
954
955
956
957 template <typename T>
get_transpose(SparseMatrix<T> & dest)958 void PetscMatrix<T>::get_transpose (SparseMatrix<T> & dest) const
959 {
960 // Make sure the SparseMatrix passed in is really a PetscMatrix
961 PetscMatrix<T> & petsc_dest = cast_ref<PetscMatrix<T> &>(dest);
962
963 // If we aren't reusing the matrix then need to clear dest,
964 // otherwise we get a memory leak
965 if (&petsc_dest != this)
966 dest.clear();
967
968 PetscErrorCode ierr;
969 if (&petsc_dest == this)
970 // The MAT_REUSE_MATRIX flag was replaced by MAT_INPLACE_MATRIX
971 // in PETSc 3.7.0
972 #if PETSC_VERSION_LESS_THAN(3,7,0)
973 ierr = MatTranspose(_mat,MAT_REUSE_MATRIX,&petsc_dest._mat);
974 #else
975 ierr = MatTranspose(_mat, MAT_INPLACE_MATRIX, &petsc_dest._mat);
976 #endif
977 else
978 ierr = MatTranspose(_mat,MAT_INITIAL_MATRIX,&petsc_dest._mat);
979 LIBMESH_CHKERR(ierr);
980
981 // Specify that the transposed matrix is initialized and close it.
982 petsc_dest._is_initialized = true;
983 petsc_dest.close();
984 }
985
986
987
988
989
990 template <typename T>
close()991 void PetscMatrix<T>::close ()
992 {
993 semiparallel_only();
994
995 // BSK - 1/19/2004
996 // strictly this check should be OK, but it seems to
997 // fail on matrix-free matrices. Do they falsely
998 // state they are assembled? Check with the developers...
999 // if (this->closed())
1000 // return;
1001
1002 PetscErrorCode ierr=0;
1003
1004 ierr = MatAssemblyBegin (_mat, MAT_FINAL_ASSEMBLY);
1005 LIBMESH_CHKERR(ierr);
1006 ierr = MatAssemblyEnd (_mat, MAT_FINAL_ASSEMBLY);
1007 LIBMESH_CHKERR(ierr);
1008 }
1009
1010 template <typename T>
flush()1011 void PetscMatrix<T>::flush ()
1012 {
1013 semiparallel_only();
1014
1015 PetscErrorCode ierr=0;
1016
1017 ierr = MatAssemblyBegin (_mat, MAT_FLUSH_ASSEMBLY);
1018 LIBMESH_CHKERR(ierr);
1019 ierr = MatAssemblyEnd (_mat, MAT_FLUSH_ASSEMBLY);
1020 LIBMESH_CHKERR(ierr);
1021 }
1022
1023
1024
1025 template <typename T>
m()1026 numeric_index_type PetscMatrix<T>::m () const
1027 {
1028 libmesh_assert (this->initialized());
1029
1030 PetscInt petsc_m=0, petsc_n=0;
1031 PetscErrorCode ierr=0;
1032
1033 ierr = MatGetSize (_mat, &petsc_m, &petsc_n);
1034 LIBMESH_CHKERR(ierr);
1035
1036 return static_cast<numeric_index_type>(petsc_m);
1037 }
1038
1039 template <typename T>
local_m()1040 numeric_index_type PetscMatrix<T>::local_m () const
1041 {
1042 libmesh_assert (this->initialized());
1043
1044 PetscInt m = 0;
1045
1046 auto ierr = MatGetLocalSize (_mat, &m, NULL);
1047 LIBMESH_CHKERR(ierr);
1048
1049 return static_cast<numeric_index_type>(m);
1050 }
1051
1052 template <typename T>
n()1053 numeric_index_type PetscMatrix<T>::n () const
1054 {
1055 libmesh_assert (this->initialized());
1056
1057 PetscInt petsc_m=0, petsc_n=0;
1058 PetscErrorCode ierr=0;
1059
1060 ierr = MatGetSize (_mat, &petsc_m, &petsc_n);
1061 LIBMESH_CHKERR(ierr);
1062
1063 return static_cast<numeric_index_type>(petsc_n);
1064 }
1065
1066 template <typename T>
local_n()1067 numeric_index_type PetscMatrix<T>::local_n () const
1068 {
1069 libmesh_assert (this->initialized());
1070
1071 PetscInt n = 0;
1072
1073 auto ierr = MatGetLocalSize (_mat, NULL, &n);
1074 LIBMESH_CHKERR(ierr);
1075
1076 return static_cast<numeric_index_type>(n);
1077 }
1078
1079 template <typename T>
get_local_size(numeric_index_type & m,numeric_index_type & n)1080 void PetscMatrix<T>::get_local_size (numeric_index_type & m,
1081 numeric_index_type & n) const
1082 {
1083 libmesh_assert (this->initialized());
1084
1085 PetscInt petsc_m = 0, petsc_n = 0;
1086
1087 auto ierr = MatGetLocalSize (_mat, &petsc_m, &petsc_n);
1088 LIBMESH_CHKERR(ierr);
1089
1090 m = static_cast<numeric_index_type>(petsc_m);
1091 n = static_cast<numeric_index_type>(petsc_n);
1092 }
1093
1094 template <typename T>
row_start()1095 numeric_index_type PetscMatrix<T>::row_start () const
1096 {
1097 libmesh_assert (this->initialized());
1098
1099 PetscInt start=0, stop=0;
1100 PetscErrorCode ierr=0;
1101
1102 ierr = MatGetOwnershipRange(_mat, &start, &stop);
1103 LIBMESH_CHKERR(ierr);
1104
1105 return static_cast<numeric_index_type>(start);
1106 }
1107
1108
1109
1110 template <typename T>
row_stop()1111 numeric_index_type PetscMatrix<T>::row_stop () const
1112 {
1113 libmesh_assert (this->initialized());
1114
1115 PetscInt start=0, stop=0;
1116 PetscErrorCode ierr=0;
1117
1118 ierr = MatGetOwnershipRange(_mat, &start, &stop);
1119 LIBMESH_CHKERR(ierr);
1120
1121 return static_cast<numeric_index_type>(stop);
1122 }
1123
1124
1125
1126 template <typename T>
set(const numeric_index_type i,const numeric_index_type j,const T value)1127 void PetscMatrix<T>::set (const numeric_index_type i,
1128 const numeric_index_type j,
1129 const T value)
1130 {
1131 libmesh_assert (this->initialized());
1132
1133 PetscErrorCode ierr=0;
1134 PetscInt i_val=i, j_val=j;
1135
1136 PetscScalar petsc_value = static_cast<PetscScalar>(value);
1137 ierr = MatSetValues(_mat, 1, &i_val, 1, &j_val,
1138 &petsc_value, INSERT_VALUES);
1139 LIBMESH_CHKERR(ierr);
1140 }
1141
1142
1143
1144 template <typename T>
add(const numeric_index_type i,const numeric_index_type j,const T value)1145 void PetscMatrix<T>::add (const numeric_index_type i,
1146 const numeric_index_type j,
1147 const T value)
1148 {
1149 libmesh_assert (this->initialized());
1150
1151 PetscErrorCode ierr=0;
1152 PetscInt i_val=i, j_val=j;
1153
1154 PetscScalar petsc_value = static_cast<PetscScalar>(value);
1155 ierr = MatSetValues(_mat, 1, &i_val, 1, &j_val,
1156 &petsc_value, ADD_VALUES);
1157 LIBMESH_CHKERR(ierr);
1158 }
1159
1160
1161
1162 template <typename T>
add_matrix(const DenseMatrix<T> & dm,const std::vector<numeric_index_type> & dof_indices)1163 void PetscMatrix<T>::add_matrix(const DenseMatrix<T> & dm,
1164 const std::vector<numeric_index_type> & dof_indices)
1165 {
1166 this->add_matrix (dm, dof_indices, dof_indices);
1167 }
1168
1169
1170
1171
1172
1173
1174
1175 template <typename T>
add(const T a_in,const SparseMatrix<T> & X_in)1176 void PetscMatrix<T>::add (const T a_in, const SparseMatrix<T> & X_in)
1177 {
1178 libmesh_assert (this->initialized());
1179
1180 // sanity check. but this cannot avoid
1181 // crash due to incompatible sparsity structure...
1182 libmesh_assert_equal_to (this->m(), X_in.m());
1183 libmesh_assert_equal_to (this->n(), X_in.n());
1184
1185 PetscScalar a = static_cast<PetscScalar> (a_in);
1186 const PetscMatrix<T> * X = cast_ptr<const PetscMatrix<T> *> (&X_in);
1187
1188 libmesh_assert (X);
1189
1190 PetscErrorCode ierr=0;
1191
1192 // the matrix from which we copy the values has to be assembled/closed
1193 libmesh_assert(X->closed());
1194
1195 semiparallel_only();
1196
1197 ierr = MatAXPY(_mat, a, X->_mat, DIFFERENT_NONZERO_PATTERN);
1198 LIBMESH_CHKERR(ierr);
1199 }
1200
1201
1202 template <typename T>
matrix_matrix_mult(SparseMatrix<T> & X_in,SparseMatrix<T> & Y_out)1203 void PetscMatrix<T>::matrix_matrix_mult (SparseMatrix<T> & X_in, SparseMatrix<T> & Y_out)
1204 {
1205 libmesh_assert (this->initialized());
1206
1207 // sanity check.
1208 libmesh_assert_equal_to (this->n(), X_in.m());
1209 libmesh_assert_equal_to (this->m(), Y_out.m());
1210 libmesh_assert_equal_to (X_in.n(), Y_out.n());
1211
1212 const PetscMatrix<T> * X = cast_ptr<const PetscMatrix<T> *> (&X_in);
1213 PetscMatrix<T> * Y = cast_ptr<PetscMatrix<T> *> (&Y_out);
1214
1215 PetscErrorCode ierr = 0;
1216
1217 // the matrix from which we copy the values has to be assembled/closed
1218 libmesh_assert(X->closed());
1219
1220 semiparallel_only();
1221
1222 ierr = MatMatMult(_mat, X->_mat, MAT_INITIAL_MATRIX, PETSC_DEFAULT, &Y->_mat);
1223 LIBMESH_CHKERR(ierr);
1224 }
1225
1226 template <typename T>
1227 void
add_sparse_matrix(const SparseMatrix<T> & spm,const std::map<numeric_index_type,numeric_index_type> & row_ltog,const std::map<numeric_index_type,numeric_index_type> & col_ltog,const T scalar)1228 PetscMatrix<T>::add_sparse_matrix (const SparseMatrix<T> & spm,
1229 const std::map<numeric_index_type, numeric_index_type> & row_ltog,
1230 const std::map<numeric_index_type, numeric_index_type> & col_ltog,
1231 const T scalar)
1232 {
1233 // size of spm is usually greater than row_ltog and col_ltog in parallel as the indices are owned by the processor
1234 // also, we should allow adding certain parts of spm to _mat
1235 libmesh_assert_greater_equal(spm.m(), row_ltog.size());
1236 libmesh_assert_greater_equal(spm.n(), col_ltog.size());
1237
1238 // make sure matrix has larger size than spm
1239 libmesh_assert_greater_equal(this->m(), spm.m());
1240 libmesh_assert_greater_equal(this->n(), spm.n());
1241
1242 if (!this->closed())
1243 this->close();
1244
1245 PetscErrorCode ierr = 0;
1246
1247 auto pscm = cast_ptr<const PetscMatrix<T> *>(&spm);
1248
1249 PetscInt ncols = 0;
1250
1251 const PetscInt * lcols;
1252 const PetscScalar * vals;
1253
1254 std::vector<PetscInt> gcols;
1255 std::vector<PetscScalar> values;
1256
1257 for (auto ltog : row_ltog)
1258 {
1259 PetscInt grow[] = {static_cast<PetscInt>(ltog.second)}; // global row index
1260
1261 ierr = MatGetRow(pscm->_mat, static_cast<PetscInt>(ltog.first), &ncols, &lcols, &vals);
1262 LIBMESH_CHKERR(ierr);
1263
1264 // get global indices (gcols) from lcols, increment values = vals*scalar
1265 gcols.resize(ncols);
1266 values.resize(ncols);
1267 for (auto i : index_range(gcols))
1268 {
1269 gcols[i] = libmesh_map_find(col_ltog, lcols[i]);
1270 values[i] = scalar * vals[i];
1271 }
1272
1273 ierr = MatSetValues(_mat, 1, grow, ncols, gcols.data(), values.data(), ADD_VALUES);
1274 LIBMESH_CHKERR(ierr);
1275 ierr = MatRestoreRow(pscm->_mat, static_cast<PetscInt>(ltog.first), &ncols, &lcols, &vals);
1276 LIBMESH_CHKERR(ierr);
1277 }
1278 // Note: We are not closing the matrix because it is expensive to do so when adding multiple sparse matrices.
1279 // Remember to manually close the matrix once all changes to the matrix have been made.
1280 }
1281
1282 template <typename T>
operator()1283 T PetscMatrix<T>::operator () (const numeric_index_type i_in,
1284 const numeric_index_type j_in) const
1285 {
1286 libmesh_assert (this->initialized());
1287
1288 // PETSc 2.2.1 & newer
1289 const PetscScalar * petsc_row;
1290 const PetscInt * petsc_cols;
1291
1292 // If the entry is not in the sparse matrix, it is 0.
1293 T value=0.;
1294
1295 PetscErrorCode
1296 ierr=0;
1297 PetscInt
1298 ncols=0,
1299 i_val=static_cast<PetscInt>(i_in),
1300 j_val=static_cast<PetscInt>(j_in);
1301
1302
1303 // the matrix needs to be closed for this to work
1304 // this->close();
1305 // but closing it is a semiparallel operation; we want operator()
1306 // to run on one processor.
1307 libmesh_assert(this->closed());
1308
1309 ierr = MatGetRow(_mat, i_val, &ncols, &petsc_cols, &petsc_row);
1310 LIBMESH_CHKERR(ierr);
1311
1312 // Perform a binary search to find the contiguous index in
1313 // petsc_cols (resp. petsc_row) corresponding to global index j_val
1314 std::pair<const PetscInt *, const PetscInt *> p =
1315 std::equal_range (petsc_cols, petsc_cols + ncols, j_val);
1316
1317 // Found an entry for j_val
1318 if (p.first != p.second)
1319 {
1320 // The entry in the contiguous row corresponding
1321 // to the j_val column of interest
1322 const std::size_t j =
1323 std::distance (const_cast<PetscInt *>(petsc_cols),
1324 const_cast<PetscInt *>(p.first));
1325
1326 libmesh_assert_less (j, ncols);
1327 libmesh_assert_equal_to (petsc_cols[j], j_val);
1328
1329 value = static_cast<T> (petsc_row[j]);
1330 }
1331
1332 ierr = MatRestoreRow(_mat, i_val,
1333 &ncols, &petsc_cols, &petsc_row);
1334 LIBMESH_CHKERR(ierr);
1335
1336 return value;
1337 }
1338
1339 template <typename T>
get_row(numeric_index_type i_in,std::vector<numeric_index_type> & indices,std::vector<T> & values)1340 void PetscMatrix<T>::get_row (numeric_index_type i_in,
1341 std::vector<numeric_index_type> & indices,
1342 std::vector<T> & values) const
1343 {
1344 libmesh_assert (this->initialized());
1345
1346 const PetscScalar * petsc_row;
1347 const PetscInt * petsc_cols;
1348
1349 PetscErrorCode ierr=0;
1350 PetscInt
1351 ncols=0,
1352 i_val = static_cast<PetscInt>(i_in);
1353
1354 // the matrix needs to be closed for this to work
1355 // this->close();
1356 // but closing it is a semiparallel operation; we want operator()
1357 // to run on one processor.
1358 libmesh_assert(this->closed());
1359
1360 // PETSc makes no effort at being thread safe. Helgrind complains about
1361 // possible data races even just in PetscFunctionBegin (due to things
1362 // like stack counter incrementing). Perhaps we could ignore
1363 // this, but there are legitimate data races for Mat data members like
1364 // mat->getrowactive between MatGetRow and MatRestoreRow. Moreover,
1365 // there could be a write into mat->rowvalues during MatGetRow from
1366 // one thread while we are attempting to read from mat->rowvalues
1367 // (through petsc_cols) during data copy in another thread. So
1368 // the safe thing to do is to lock the whole method
1369
1370 #ifdef LIBMESH_HAVE_CXX11_THREAD
1371 std::lock_guard<std::mutex>
1372 #else
1373 Threads::spin_mutex::scoped_lock
1374 #endif
1375 lock(_petsc_matrix_mutex);
1376
1377 ierr = MatGetRow(_mat, i_val, &ncols, &petsc_cols, &petsc_row);
1378 LIBMESH_CHKERR(ierr);
1379
1380 // Copy the data
1381 indices.resize(static_cast<std::size_t>(ncols));
1382 values.resize(static_cast<std::size_t>(ncols));
1383
1384 for (auto i : index_range(indices))
1385 {
1386 indices[i] = static_cast<numeric_index_type>(petsc_cols[i]);
1387 values[i] = static_cast<T>(petsc_row[i]);
1388 }
1389
1390 ierr = MatRestoreRow(_mat, i_val,
1391 &ncols, &petsc_cols, &petsc_row);
1392 LIBMESH_CHKERR(ierr);
1393 }
1394
1395
1396 template <typename T>
closed()1397 bool PetscMatrix<T>::closed() const
1398 {
1399 libmesh_assert (this->initialized());
1400
1401 PetscErrorCode ierr=0;
1402 PetscBool assembled;
1403
1404 ierr = MatAssembled(_mat, &assembled);
1405 LIBMESH_CHKERR(ierr);
1406
1407 return (assembled == PETSC_TRUE);
1408 }
1409
1410 template <typename T>
set_destroy_mat_on_exit(bool destroy)1411 void PetscMatrix<T>::set_destroy_mat_on_exit(bool destroy)
1412 {
1413 this->_destroy_mat_on_exit = destroy;
1414 }
1415
1416
1417 template <typename T>
swap(PetscMatrix<T> & m_in)1418 void PetscMatrix<T>::swap(PetscMatrix<T> & m_in)
1419 {
1420 std::swap(_mat, m_in._mat);
1421 std::swap(_destroy_mat_on_exit, m_in._destroy_mat_on_exit);
1422 }
1423
1424
1425
1426 //------------------------------------------------------------------
1427 // Explicit instantiations
1428 template class PetscMatrix<Number>;
1429
1430 } // namespace libMesh
1431
1432
1433 #endif // #ifdef LIBMESH_HAVE_PETSC
1434