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