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 : #ifndef LIBMESH_PETSC_MATRIX_H
21 : #define LIBMESH_PETSC_MATRIX_H
22 :
23 : #include "libmesh/libmesh_common.h"
24 :
25 : #ifdef LIBMESH_HAVE_PETSC
26 :
27 : // Local includes
28 : #include "libmesh/petsc_matrix_base.h"
29 : #include "libmesh/petsc_macro.h"
30 : #include "libmesh/petsc_solver_exception.h"
31 :
32 : // C++ includes
33 : #include <algorithm>
34 : #ifdef LIBMESH_HAVE_CXX11_THREAD
35 : #include <atomic>
36 : #include <mutex>
37 : #endif
38 :
39 : class PetscMatrixTest;
40 :
41 : namespace libMesh
42 : {
43 :
44 : // Forward Declarations
45 : template <typename T> class DenseMatrix;
46 :
47 : enum PetscMatrixType : int {
48 : AIJ=0,
49 : HYPRE};
50 :
51 : /**
52 : * This class provides a nice interface to the PETSc C-based AIJ data
53 : * structures for parallel, sparse matrices. All overridden virtual
54 : * functions are documented in sparse_matrix.h.
55 : *
56 : * \author Benjamin S. Kirk
57 : * \date 2002
58 : * \brief SparseMatrix interface to PETSc Mat.
59 : */
60 : template <typename T>
61 3056 : class PetscMatrix final : public PetscMatrixBase<T>
62 : {
63 : public:
64 : /**
65 : * Constructor; initializes the matrix to be empty, without any
66 : * structure, i.e. the matrix is not usable at all. This
67 : * constructor is therefore only useful for matrices which are
68 : * members of a class. All other matrices should be created at a
69 : * point in the data flow where all necessary information is
70 : * available.
71 : *
72 : * You have to initialize the matrix before usage with \p init(...).
73 : */
74 : explicit
75 : PetscMatrix (const Parallel::Communicator & comm_in);
76 :
77 : /**
78 : * Constructor. Creates a PetscMatrix assuming you already have a
79 : * valid Mat object. In this case, m may not be destroyed by the
80 : * PetscMatrix destructor when this object goes out of scope. This
81 : * allows ownership of m to remain with the original creator, and to
82 : * simply provide additional functionality with the PetscMatrix.
83 : */
84 : explicit
85 : PetscMatrix (Mat m,
86 : const Parallel::Communicator & comm_in,
87 : bool destroy_on_exit = false);
88 :
89 : /**
90 : * Constructor. Creates and initializes a PetscMatrix with the given
91 : * structure. See \p init(...) for a description of the parameters.
92 : */
93 : explicit
94 : PetscMatrix (const Parallel::Communicator & comm_in,
95 : const numeric_index_type m,
96 : const numeric_index_type n,
97 : const numeric_index_type m_l,
98 : const numeric_index_type n_l,
99 : const numeric_index_type n_nz=30,
100 : const numeric_index_type n_oz=10,
101 : const numeric_index_type blocksize=1);
102 :
103 : /**
104 : * This class manages a C-style struct (Mat) manually, so we
105 : * don't want to allow any automatic copy/move functions to be
106 : * generated, and we can't default the destructor.
107 : */
108 : PetscMatrix (PetscMatrix &&) = delete;
109 : PetscMatrix (const PetscMatrix &) = delete;
110 : PetscMatrix & operator= (PetscMatrix &&) = delete;
111 : virtual ~PetscMatrix ();
112 :
113 : PetscMatrix & operator= (const PetscMatrix &);
114 : virtual SparseMatrix<T> & operator= (const SparseMatrix<T> & v) override;
115 :
116 : /**
117 : * Initialize a PETSc matrix.
118 : *
119 : * \param m The global number of rows.
120 : * \param n The global number of columns.
121 : * \param m_l The local number of rows.
122 : * \param n_l The local number of columns.
123 : * \param n_nz The number of nonzeros in each row of the DIAGONAL portion of the local submatrix.
124 : * \param n_oz The number of nonzeros in each row of the OFF-DIAGONAL portion of the local submatrix.
125 : * \param blocksize Optional value indicating dense coupled blocks for systems with multiple variables all of the same type.
126 : */
127 : virtual void init (const numeric_index_type m,
128 : const numeric_index_type n,
129 : const numeric_index_type m_l,
130 : const numeric_index_type n_l,
131 : const numeric_index_type n_nz=30,
132 : const numeric_index_type n_oz=10,
133 : const numeric_index_type blocksize=1) override;
134 :
135 : /**
136 : * Initialize a PETSc matrix.
137 : *
138 : * \param m The global number of rows.
139 : * \param n The global number of columns.
140 : * \param m_l The local number of rows.
141 : * \param n_l The local number of columns.
142 : * \param n_nz array containing the number of nonzeros in each row of the DIAGONAL portion of the local submatrix.
143 : * \param n_oz Array containing the number of nonzeros in each row of the OFF-DIAGONAL portion of the local submatrix.
144 : * \param blocksize Optional value indicating dense coupled blocks for systems with multiple variables all of the same type.
145 : */
146 : void init (const numeric_index_type m,
147 : const numeric_index_type n,
148 : const numeric_index_type m_l,
149 : const numeric_index_type n_l,
150 : const std::vector<numeric_index_type> & n_nz,
151 : const std::vector<numeric_index_type> & n_oz,
152 : const numeric_index_type blocksize=1);
153 :
154 : virtual void init (ParallelType = PARALLEL) override;
155 :
156 : /**
157 : * Update the sparsity pattern based on \p dof_map, and set the matrix
158 : * to zero. This is useful in cases where the sparsity pattern changes
159 : * during a computation.
160 : */
161 : void update_preallocation_and_zero();
162 :
163 : /**
164 : * Reset matrix to use the original nonzero pattern provided by users
165 : */
166 : void reset_preallocation();
167 :
168 : virtual void zero () override;
169 :
170 : virtual std::unique_ptr<SparseMatrix<T>> zero_clone () const override;
171 :
172 : virtual std::unique_ptr<SparseMatrix<T>> clone () const override;
173 :
174 : virtual void zero_rows (std::vector<numeric_index_type> & rows, T diag_value = 0.0) override;
175 :
176 : virtual void flush () override;
177 :
178 : virtual void set (const numeric_index_type i,
179 : const numeric_index_type j,
180 : const T value) override;
181 :
182 : virtual void add (const numeric_index_type i,
183 : const numeric_index_type j,
184 : const T value) override;
185 :
186 : virtual void add_matrix (const DenseMatrix<T> & dm,
187 : const std::vector<numeric_index_type> & rows,
188 : const std::vector<numeric_index_type> & cols) override;
189 :
190 : virtual void add_matrix (const DenseMatrix<T> & dm,
191 : const std::vector<numeric_index_type> & dof_indices) override;
192 :
193 : virtual void add_block_matrix (const DenseMatrix<T> & dm,
194 : const std::vector<numeric_index_type> & brows,
195 : const std::vector<numeric_index_type> & bcols) override;
196 :
197 0 : virtual void add_block_matrix (const DenseMatrix<T> & dm,
198 : const std::vector<numeric_index_type> & dof_indices) override
199 0 : { this->add_block_matrix (dm, dof_indices, dof_indices); }
200 :
201 : /**
202 : * Compute A += a*X for scalar \p a, matrix \p X.
203 : *
204 : * \note The matrices A and X need to have the same nonzero pattern,
205 : * otherwise \p PETSc will crash!
206 : *
207 : * \note It is advisable to not only allocate appropriate memory
208 : * with \p init(), but also explicitly zero the terms of \p this
209 : * whenever you add a non-zero value to \p X.
210 : *
211 : * \note \p X will be closed, if not already done, before performing
212 : * any work.
213 : */
214 : virtual void add (const T a, const SparseMatrix<T> & X) override;
215 :
216 : /**
217 : * Compute Y = A*X for matrix \p X.
218 : * Set reuse = true if this->_mat and X have the same nonzero pattern as before,
219 : * and Y is obtained from a previous call to this function with reuse = false
220 : */
221 : virtual void matrix_matrix_mult (SparseMatrix<T> & X, SparseMatrix<T> & Y, bool reuse = false) override;
222 :
223 : /**
224 : * Add \p scalar* \p spm to the rows and cols of this matrix (A):
225 : * A(rows[i], cols[j]) += scalar * spm(i,j)
226 : */
227 : virtual
228 : void add_sparse_matrix (const SparseMatrix<T> & spm,
229 : const std::map<numeric_index_type,numeric_index_type> & row_ltog,
230 : const std::map<numeric_index_type,numeric_index_type> & col_ltog,
231 : const T scalar) override;
232 :
233 : virtual T operator () (const numeric_index_type i,
234 : const numeric_index_type j) const override;
235 :
236 : virtual Real l1_norm () const override;
237 :
238 : virtual Real frobenius_norm () const;
239 :
240 : virtual Real linfty_norm () const override;
241 :
242 : /**
243 : * Print the contents of the matrix to the screen with the PETSc
244 : * viewer. This function only allows printing to standard out since
245 : * we have limited ourselves to one PETSc implementation for
246 : * writing.
247 : */
248 : virtual void print_personal(std::ostream & os=libMesh::out) const override;
249 :
250 : virtual void print_matlab(const std::string & name = "") const override;
251 :
252 : virtual void print_petsc_binary(const std::string & filename) override;
253 :
254 : virtual void print_petsc_hdf5(const std::string & filename) override;
255 :
256 : virtual void read_petsc_binary(const std::string & filename) override;
257 :
258 : virtual void read_petsc_hdf5(const std::string & filename) override;
259 :
260 : virtual void get_diagonal (NumericVector<T> & dest) const override;
261 :
262 : virtual void get_transpose (SparseMatrix<T> & dest) const override;
263 :
264 : virtual
265 : void get_row(numeric_index_type i,
266 : std::vector<numeric_index_type> & indices,
267 : std::vector<T> & values) const override;
268 :
269 : /**
270 : * Similar to the `create_submatrix` function, this function creates a \p
271 : * submatrix which is defined by the indices given in the \p rows
272 : * and \p cols vectors.
273 : * Note: Both \p rows and \p cols can be unsorted;
274 : * Use the above function for better efficiency if your indices are sorted;
275 : * \p rows and \p cols can contain indices that are owned by other processors.
276 : */
277 : virtual void create_submatrix_nosort(SparseMatrix<T> & submatrix,
278 : const std::vector<numeric_index_type> & rows,
279 : const std::vector<numeric_index_type> & cols) const override;
280 :
281 : virtual void scale(const T scale) override;
282 :
283 : #if PETSC_RELEASE_GREATER_EQUALS(3, 23, 0)
284 : /**
285 : * Creates a copy of the current hash table matrix and then performs assembly. This is very useful
286 : * in cases where you are not done filling this matrix but want to be able to read the current
287 : * state of it
288 : */
289 : std::unique_ptr<PetscMatrix<T>> copy_from_hash();
290 : #endif
291 :
292 : virtual bool supports_hash_table() const override;
293 :
294 : virtual void restore_original_nonzero_pattern() override;
295 :
296 : protected:
297 : /**
298 : * Perform matrix initialization steps sans preallocation
299 : * @param m The global number of rows
300 : * @param n The global number of columns
301 : * @param m_l The local number of rows
302 : * @param n_l The local number of columns
303 : * @param blocksize The matrix block size
304 : */
305 : void init_without_preallocation (numeric_index_type m,
306 : numeric_index_type n,
307 : numeric_index_type m_l,
308 : numeric_index_type n_l,
309 : numeric_index_type blocksize);
310 :
311 : /*
312 : * Performs matrix preallcation
313 : * \param m_l The local number of rows.
314 : * \param n_nz array containing the number of nonzeros in each row of the DIAGONAL portion of the local submatrix.
315 : * \param n_oz Array containing the number of nonzeros in each row of the OFF-DIAGONAL portion of the local submatrix.
316 : * \param blocksize Optional value indicating dense coupled blocks for systems with multiple variables all of the same */
317 : void preallocate(numeric_index_type m_l,
318 : const std::vector<numeric_index_type> & n_nz,
319 : const std::vector<numeric_index_type> & n_oz,
320 : numeric_index_type blocksize);
321 :
322 : /**
323 : * Finish up the initialization process. This method does a few things which include
324 : * - Setting the option to make new nonzeroes an error (otherwise users will just have a silent
325 : (often huge) performance penalty
326 : * - Marking the matrix as initialized
327 : */
328 : void finish_initialization();
329 :
330 : /**
331 : * This function either creates or re-initializes a matrix called \p
332 : * submatrix which is defined by the indices given in the \p rows
333 : * and \p cols vectors.
334 : *
335 : * This function is implemented in terms of MatGetSubMatrix(). The
336 : * \p reuse_submatrix parameter determines whether or not PETSc will
337 : * treat \p submatrix as one which has already been used (had memory
338 : * allocated) or as a new matrix.
339 : */
340 : virtual void _get_submatrix(SparseMatrix<T> & submatrix,
341 : const std::vector<numeric_index_type> & rows,
342 : const std::vector<numeric_index_type> & cols,
343 : const bool reuse_submatrix) const override;
344 :
345 : void _petsc_viewer(const std::string & filename,
346 : PetscViewerType viewertype,
347 : PetscFileMode filemode);
348 :
349 : PetscMatrixType _mat_type;
350 :
351 : private:
352 : #ifdef LIBMESH_HAVE_CXX11_THREAD
353 : mutable std::mutex _petsc_matrix_mutex;
354 : #else
355 : mutable Threads::spin_mutex _petsc_matrix_mutex;
356 : #endif
357 :
358 : /**
359 : * \returns A norm of the matrix, where the type of norm to compute is
360 : * determined by the template parameter N of the PETSc-defined type NormType.
361 : * The valid template arguments are NORM_1, NORM_FROBENIUS and NORM_INFINITY,
362 : * as used to define l1_norm(), frobenius_norm() and linfty_norm().
363 : */
364 : template <NormType N> Real norm () const;
365 :
366 : friend class ::PetscMatrixTest;
367 : };
368 :
369 : } // namespace libMesh
370 :
371 : #endif // #ifdef LIBMESH_HAVE_PETSC
372 : #endif // LIBMESH_PETSC_MATRIX_H
|
