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_SPARSE_MATRIX_H
21 : #define LIBMESH_SPARSE_MATRIX_H
22 :
23 :
24 : // Local includes
25 : #include "libmesh/libmesh.h"
26 : #include "libmesh/libmesh_common.h"
27 : #include "libmesh/reference_counted_object.h"
28 : #include "libmesh/id_types.h"
29 : #include "libmesh/parallel_object.h"
30 : #include "libmesh/enum_parallel_type.h" // PARALLEL
31 : #include "libmesh/enum_matrix_build_type.h" // AUTOMATIC
32 : #include "libmesh/enum_solver_package.h"
33 :
34 : // C++ includes
35 : #include <cstddef>
36 : #include <iomanip>
37 : #include <vector>
38 : #include <memory>
39 :
40 : namespace libMesh
41 : {
42 :
43 : // forward declarations
44 : template <typename T> class SparseMatrix;
45 : template <typename T> class DenseMatrix;
46 : class DofMap;
47 : namespace SparsityPattern {
48 : class Build;
49 : class Graph;
50 : }
51 : template <typename T> class NumericVector;
52 :
53 : // This template helper function must be declared before it
54 : // can be defined below.
55 : template <typename T>
56 : std::ostream & operator << (std::ostream & os, const SparseMatrix<T> & m);
57 :
58 :
59 : /**
60 : * Generic sparse matrix. This class contains pure virtual members
61 : * that must be overridden in derived classes. Using a common base
62 : * class allows for uniform access to sparse matrices from various
63 : * different solver packages in different formats.
64 : *
65 : * \author Benjamin S. Kirk
66 : * \date 2003
67 : */
68 : template <typename T>
69 : class SparseMatrix : public ReferenceCountedObject<SparseMatrix<T>>,
70 : public ParallelObject
71 : {
72 : public:
73 : /**
74 : * Constructor; initializes the matrix to be empty, without any
75 : * structure, i.e. the matrix is not usable at all. This
76 : * constructor is therefore only useful for matrices which are
77 : * members of a class. All other matrices should be created at a
78 : * point in the data flow where all necessary information is
79 : * available.
80 : *
81 : * You have to initialize the matrix before usage with
82 : * \p init(...).
83 : */
84 : explicit
85 : SparseMatrix (const Parallel::Communicator & comm);
86 :
87 : /**
88 : * This _looks_ like a copy assignment operator, but note that,
89 : * unlike normal copy assignment operators, it is pure virtual. This
90 : * function should be overridden in derived classes so that they can
91 : * be copied correctly via references to the base class. This design
92 : * usually isn't a good idea in general, but in this context it
93 : * works because we usually don't have a mix of different kinds of
94 : * SparseMatrix active in the library at a single time.
95 : *
96 : * \returns A reference to *this as the base type.
97 : */
98 0 : virtual SparseMatrix<T> & operator= (const SparseMatrix<T> &)
99 : {
100 0 : libmesh_not_implemented();
101 : return *this;
102 : }
103 :
104 : /**
105 : * These 3 special functions can be defaulted for this class, as it
106 : * does not manage any memory itself.
107 : */
108 : SparseMatrix (SparseMatrix &&) = default;
109 284 : SparseMatrix (const SparseMatrix &) = default;
110 : SparseMatrix & operator= (SparseMatrix &&) = default;
111 :
112 : /**
113 : * While this class doesn't manage any memory, the derived class might and
114 : * users may be deleting through a pointer to this base class
115 : */
116 2674 : virtual ~SparseMatrix() = default;
117 :
118 : /**
119 : * Builds a \p SparseMatrix<T> using the linear solver package specified by
120 : * \p solver_package
121 : */
122 : static std::unique_ptr<SparseMatrix<T>>
123 : build(const Parallel::Communicator & comm,
124 : const SolverPackage solver_package = libMesh::default_solver_package(),
125 : const MatrixBuildType matrix_build_type = MatrixBuildType::AUTOMATIC);
126 :
127 : virtual SolverPackage solver_package() = 0;
128 :
129 : /**
130 : * \returns \p true if the matrix has been initialized,
131 : * \p false otherwise.
132 : */
133 3865791 : virtual bool initialized() const { return _is_initialized; }
134 :
135 : /**
136 : * Set a pointer to the \p DofMap to use. If a separate sparsity
137 : * pattern is not being used, use the one from the DofMap.
138 : *
139 : * The lifetime of \p dof_map must exceed the lifetime of \p this.
140 : */
141 : void attach_dof_map (const DofMap & dof_map);
142 :
143 : /**
144 : * Set a pointer to a sparsity pattern to use. Useful in cases
145 : * where a matrix requires a wider (or for efficiency narrower)
146 : * pattern than most matrices in the system, or in cases where no
147 : * system sparsity pattern is being calculated by the DofMap.
148 : *
149 : * The lifetime of \p sp must exceed the lifetime of \p this.
150 : */
151 : void attach_sparsity_pattern (const SparsityPattern::Build & sp);
152 :
153 : /**
154 : * \returns \p true if this sparse matrix format needs to be fed the
155 : * graph of the sparse matrix.
156 : *
157 : * This is true for \p LaspackMatrix, but not \p PetscMatrixBase subclasses. In
158 : * the case where the full graph is not required, we can efficiently
159 : * approximate it to provide a good estimate of the required size of
160 : * the sparse matrix.
161 : */
162 40374 : virtual bool need_full_sparsity_pattern() const
163 40374 : { return false; }
164 :
165 : /**
166 : * @returns Whether this matrix needs the sparsity pattern computed by the \p DofMap
167 : */
168 14396 : virtual bool require_sparsity_pattern() const { return !this->use_hash_table(); }
169 :
170 : /**
171 : * Updates the matrix sparsity pattern. When your \p SparseMatrix<T>
172 : * implementation does not need this data, simply do not override
173 : * this method.
174 : */
175 0 : virtual void update_sparsity_pattern (const SparsityPattern::Graph &) {}
176 :
177 : /**
178 : * Initialize SparseMatrix with the specified sizes.
179 : *
180 : * \param m The global number of rows.
181 : * \param n The global number of columns.
182 : * \param m_l The local number of rows.
183 : * \param n_l The local number of columns.
184 : * \param nnz The number of on-diagonal nonzeros per row (defaults to 30).
185 : * \param noz The number of off-diagonal nonzeros per row (defaults to 10).
186 : * \param blocksize Optional value indicating dense coupled blocks
187 : * for systems with multiple variables all of the same type.
188 : */
189 : virtual void init (const numeric_index_type m,
190 : const numeric_index_type n,
191 : const numeric_index_type m_l,
192 : const numeric_index_type n_l,
193 : const numeric_index_type nnz=30,
194 : const numeric_index_type noz=10,
195 : const numeric_index_type blocksize=1) = 0;
196 :
197 : /**
198 : * Initialize this matrix using the sparsity structure computed by \p dof_map.
199 : * @param type The serial/parallel/ghosted type of the matrix
200 : */
201 : virtual void init (ParallelType type = PARALLEL) = 0;
202 :
203 : /**
204 : * Restores the \p SparseMatrix<T> to a pristine state.
205 : */
206 : virtual void clear () = 0;
207 :
208 : /**
209 : * Set all entries to 0.
210 : */
211 : virtual void zero () = 0;
212 :
213 : /**
214 : * \returns A smart pointer to a copy of this matrix with the same
215 : * type, size, and partitioning, but with all zero entries.
216 : *
217 : * \note This must be overridden in the derived classes.
218 : */
219 : virtual std::unique_ptr<SparseMatrix<T>> zero_clone () const = 0;
220 :
221 : /**
222 : * \returns A smart pointer to a copy of this matrix.
223 : *
224 : * \note This must be overridden in the derived classes.
225 : */
226 : virtual std::unique_ptr<SparseMatrix<T>> clone () const = 0;
227 :
228 : /**
229 : * Sets all row entries to 0 then puts \p diag_value in the diagonal entry.
230 : */
231 : virtual void zero_rows (std::vector<numeric_index_type> & rows, T diag_value = 0.0);
232 :
233 : /**
234 : * Calls the SparseMatrix's internal assembly routines, ensuring
235 : * that the values are consistent across processors.
236 : */
237 : virtual void close () = 0;
238 :
239 : /**
240 : * For PETSc matrix , this function is similar to close but without shrinking memory.
241 : * This is useful when we want to switch between ADD_VALUES and INSERT_VALUES.
242 : * close should be called before using the matrix.
243 : */
244 0 : virtual void flush () { close(); }
245 :
246 : /**
247 : * \returns The row-dimension of the matrix.
248 : */
249 : virtual numeric_index_type m () const = 0;
250 :
251 : /**
252 : * Get the number of rows owned by this process
253 : */
254 8 : virtual numeric_index_type local_m () const { return row_stop() - row_start(); }
255 :
256 : /**
257 : * Get the number of columns owned by this process
258 : */
259 16 : virtual numeric_index_type local_n () const { return col_stop() - col_start(); }
260 :
261 : /**
262 : * \returns The column-dimension of the matrix.
263 : */
264 : virtual numeric_index_type n () const = 0;
265 :
266 : /**
267 : * \returns The index of the first matrix row stored on this
268 : * processor.
269 : */
270 : virtual numeric_index_type row_start () const = 0;
271 :
272 : /**
273 : * \returns The index of the last matrix row (+1) stored on this
274 : * processor.
275 : */
276 : virtual numeric_index_type row_stop () const = 0;
277 :
278 : /**
279 : * \returns The index of the first matrix column owned by this
280 : * processor.
281 : */
282 : virtual numeric_index_type col_start () const = 0;
283 :
284 : /**
285 : * \returns The index of the last matrix column (+1) owned by this
286 : * processor.
287 : */
288 : virtual numeric_index_type col_stop () const = 0;
289 :
290 : /**
291 : * Set the element \p (i,j) to \p value. Throws an error if the
292 : * entry does not exist. Zero values can be "stored" in non-existent
293 : * fields.
294 : */
295 : virtual void set (const numeric_index_type i,
296 : const numeric_index_type j,
297 : const T value) = 0;
298 :
299 : /**
300 : * Add \p value to the element \p (i,j). Throws an error if the
301 : * entry does not exist. Zero values can be "added" to non-existent
302 : * entries.
303 : */
304 : virtual void add (const numeric_index_type i,
305 : const numeric_index_type j,
306 : const T value) = 0;
307 :
308 : /**
309 : * Add the full matrix \p dm to the SparseMatrix. This is useful
310 : * for adding an element matrix at assembly time.
311 : */
312 : virtual void add_matrix (const DenseMatrix<T> & dm,
313 : const std::vector<numeric_index_type> & rows,
314 : const std::vector<numeric_index_type> & cols) = 0;
315 :
316 : /**
317 : * Same as \p add_matrix, but assumes the row and column maps are the same.
318 : * Thus the matrix \p dm must be square.
319 : */
320 : virtual void add_matrix (const DenseMatrix<T> & dm,
321 : const std::vector<numeric_index_type> & dof_indices) = 0;
322 :
323 : /**
324 : * Add the full matrix \p dm to the SparseMatrix. This is useful
325 : * for adding an element matrix at assembly time. The matrix is
326 : * assumed blocked, and \p brow, \p bcol correspond to the *block*
327 : * row and column indices.
328 : */
329 : virtual void add_block_matrix (const DenseMatrix<T> & dm,
330 : const std::vector<numeric_index_type> & brows,
331 : const std::vector<numeric_index_type> & bcols);
332 :
333 : /**
334 : * Same as \p add_block_matrix(), but assumes the row and column
335 : * maps are the same. Thus the matrix \p dm must be square.
336 : */
337 0 : virtual void add_block_matrix (const DenseMatrix<T> & dm,
338 : const std::vector<numeric_index_type> & dof_indices)
339 0 : { this->add_block_matrix (dm, dof_indices, dof_indices); }
340 :
341 : /**
342 : * Compute \f$ A \leftarrow A + a*X \f$ for scalar \p a, matrix \p X.
343 : */
344 : virtual void add (const T a, const SparseMatrix<T> & X) = 0;
345 :
346 : /**
347 : * Compute Y = A*X for matrix \p X.
348 : */
349 0 : virtual void matrix_matrix_mult (SparseMatrix<T> & /*X*/, SparseMatrix<T> & /*Y*/, bool /*reuse*/)
350 0 : { libmesh_not_implemented(); }
351 :
352 : /**
353 : * Add \p scalar* \p spm to the rows and cols of this matrix (A):
354 : * A(rows[i], cols[j]) += scalar * spm(i,j)
355 : */
356 0 : virtual void add_sparse_matrix (const SparseMatrix<T> & /*spm*/,
357 : const std::map<numeric_index_type, numeric_index_type> & /*rows*/,
358 : const std::map<numeric_index_type, numeric_index_type> & /*cols*/,
359 : const T /*scalar*/)
360 0 : { libmesh_not_implemented(); }
361 :
362 : /**
363 : * \returns A copy of matrix entry \p (i,j).
364 : *
365 : * \note This may be an expensive operation, and you should always
366 : * be careful where you call this function.
367 : */
368 : virtual T operator () (const numeric_index_type i,
369 : const numeric_index_type j) const = 0;
370 :
371 : /**
372 : * \returns The \f$ \ell_1 \f$-norm of the matrix, that is the max column sum:
373 : * \f$ |M|_1 = \max_{j} \sum_{i} |M_{ij}| \f$
374 : *
375 : * This is the natural matrix norm that is compatible with the
376 : * \f$ \ell_1 \f$-norm for vectors, i.e. \f$ |Mv|_1 \leq |M|_1 |v|_1 \f$.
377 : * (cf. Haemmerlin-Hoffmann : Numerische Mathematik)
378 : */
379 : virtual Real l1_norm () const = 0;
380 :
381 : /**
382 : * \returns The \f$ \ell_{\infty} \f$-norm of the matrix, that is the max row sum:
383 : *
384 : * \f$ |M|_{\infty} = \max_{i} \sum_{j} |M_{ij}| \f$
385 : *
386 : * This is the natural matrix norm that is compatible to the
387 : * \f$ \ell_{\infty} \f$-norm of vectors, i.e.
388 : * \f$ |Mv|_{\infty} \leq |M|_{\infty} |v|_{\infty} \f$.
389 : * (cf. Haemmerlin-Hoffmann : Numerische Mathematik)
390 : */
391 : virtual Real linfty_norm () const = 0;
392 :
393 : /**
394 : * \returns The l1_norm() of the difference of \p this and \p other_mat
395 : */
396 : Real l1_norm_diff (const SparseMatrix<T> & other_mat) const;
397 :
398 : /**
399 : * \returns \p true if the matrix has been assembled.
400 : */
401 : virtual bool closed() const = 0;
402 :
403 : /**
404 : * \returns the global number of non-zero entries in the matrix
405 : * sparsity pattern
406 : */
407 : virtual std::size_t n_nonzeros() const;
408 :
409 : /**
410 : * Print the contents of the matrix to the screen in a uniform
411 : * style, regardless of matrix/solver package being used.
412 : */
413 : void print(std::ostream & os=libMesh::out, const bool sparse=false) const;
414 :
415 : /**
416 : * Same as the print method above, but allows you to print to a
417 : * stream in the standard syntax.
418 : *
419 : * \code
420 : * template <typename U>
421 : * friend std::ostream & operator << (std::ostream & os, const SparseMatrix<U> & m);
422 : * \endcode
423 : *
424 : * \note The above syntax, which does not require any
425 : * prior declaration of operator<<, declares *any* instantiation of
426 : * SparseMatrix<X> is friend to *any* instantiation of
427 : * operator<<(ostream &, SparseMatrix<Y> &). It would not happen in
428 : * practice, but in principle it means that SparseMatrix<Complex>
429 : * would be friend to operator<<(ostream &, SparseMatrix<Real>).
430 : *
431 : * \note The form below, which requires a previous
432 : * declaration of the operator<<(stream &, SparseMatrix<T> &) function
433 : * (see top of this file), means that any instantiation of
434 : * SparseMatrix<T> is friend to the specialization
435 : * operator<<(ostream &, SparseMatrix<T> &), but e.g. SparseMatrix<U>
436 : * is *not* friend to the same function. So this is slightly
437 : * different to the form above...
438 : *
439 : * This method seems to be the "preferred" technique, see
440 : * http://www.parashift.com/c++-faq-lite/template-friends.html
441 : */
442 : friend std::ostream & operator << <>(std::ostream & os, const SparseMatrix<T> & m);
443 :
444 : /**
445 : * Print the contents of the matrix to the screen in a
446 : * package-personalized style, if available.
447 : */
448 : virtual void print_personal(std::ostream & os=libMesh::out) const = 0;
449 :
450 : /**
451 : * Print the contents of the matrix in Matlab's sparse matrix
452 : * format. Optionally prints the matrix to the file named \p name.
453 : * If \p name is not specified it is dumped to the screen.
454 : */
455 : virtual void print_matlab(const std::string & /*name*/ = "") const;
456 :
457 : /**
458 : * Write the contents of the matrix to a file in PETSc's binary
459 : * sparse matrix format.
460 : */
461 : virtual void print_petsc_binary(const std::string & filename);
462 :
463 : /**
464 : * Write the contents of the matrix to a file in PETSc's HDF5
465 : * sparse matrix format.
466 : */
467 : virtual void print_petsc_hdf5(const std::string & filename);
468 :
469 : /**
470 : * Read the contents of the matrix from a file, with the file format
471 : * inferred from the extension of \p filename.
472 : */
473 : virtual void read(const std::string & filename);
474 :
475 : /**
476 : * Read the contents of the matrix from a file, with the HDF5 sparse
477 : * matrix format used by CoreForm, expecing sparse matrix data in
478 : * the group given by \p groupname.
479 : *
480 : * This will be initialized with the sparsity from the file,
481 : * linearly partitioned onto the number of processors available
482 : * unless \p this matrix is pre-sized and pre-partitionsed.
483 : */
484 : virtual void read_coreform_hdf5(const std::string & filename,
485 : const std::string & groupname = "extraction");
486 :
487 : /**
488 : * Read the contents of the matrix from the Matlab-script sparse
489 : * matrix format used by PETSc.
490 : *
491 : * If the size and sparsity of the matrix in \p filename appears
492 : * consistent with the existing sparsity of \p this then the
493 : * existing parallel decomposition and sparsity will be retained.
494 : * If not, then \p this will be initialized with the sparsity from
495 : * the file, linearly partitioned onto the number of processors
496 : * available.
497 : */
498 : virtual void read_matlab(const std::string & filename);
499 :
500 : /**
501 : * Read the contents of the matrix from a file in PETSc's binary
502 : * sparse matrix format.
503 : */
504 : virtual void read_petsc_binary(const std::string & filename);
505 :
506 : /**
507 : * Read the contents of the matrix from a file in PETSc's HDF5
508 : * sparse matrix format.
509 : */
510 : virtual void read_petsc_hdf5(const std::string & filename);
511 :
512 : /**
513 : * This function creates a matrix called "submatrix" which is defined
514 : * by the row and column indices given in the "rows" and "cols" entries.
515 : * Currently this operation is only defined for the PetscMatrixBase subclasses.
516 : * Note: The \p rows and \p cols vectors need to be sorted;
517 : * Use the nosort version below if \p rows and \p cols vectors are not sorted;
518 : * The \p rows and \p cols only contain indices that are owned by this processor.
519 : */
520 594 : virtual void create_submatrix(SparseMatrix<T> & submatrix,
521 : const std::vector<numeric_index_type> & rows,
522 : const std::vector<numeric_index_type> & cols) const
523 : {
524 594 : this->_get_submatrix(submatrix,
525 : rows,
526 : cols,
527 : false); // false means DO NOT REUSE submatrix
528 594 : }
529 :
530 : /**
531 : * Similar to the above function, this function creates a \p
532 : * submatrix which is defined by the indices given in the \p rows
533 : * and \p cols vectors.
534 : * Note: Both \p rows and \p cols can be unsorted;
535 : * Use the above function for better efficiency if your indices are sorted;
536 : * \p rows and \p cols can contain indices that are owned by other processors.
537 : */
538 0 : virtual void create_submatrix_nosort(SparseMatrix<T> & /*submatrix*/,
539 : const std::vector<numeric_index_type> & /*rows*/,
540 : const std::vector<numeric_index_type> & /*cols*/) const
541 : {
542 0 : libmesh_not_implemented();
543 : }
544 :
545 : /**
546 : * This function is similar to the one above, but it allows you to reuse
547 : * the existing sparsity pattern of "submatrix" instead of reallocating
548 : * it again. This should hopefully be more efficient if you are frequently
549 : * extracting submatrices of the same size.
550 : */
551 0 : virtual void reinit_submatrix(SparseMatrix<T> & submatrix,
552 : const std::vector<numeric_index_type> & rows,
553 : const std::vector<numeric_index_type> & cols) const
554 : {
555 0 : this->_get_submatrix(submatrix,
556 : rows,
557 : cols,
558 : true); // true means REUSE submatrix
559 0 : }
560 :
561 : /**
562 : * Multiplies the matrix by the NumericVector \p arg and stores the
563 : * result in NumericVector \p dest.
564 : */
565 : void vector_mult (NumericVector<T> & dest,
566 : const NumericVector<T> & arg) const;
567 :
568 : /**
569 : * Multiplies the matrix by the NumericVector \p arg and adds the
570 : * result to the NumericVector \p dest.
571 : */
572 : void vector_mult_add (NumericVector<T> & dest,
573 : const NumericVector<T> & arg) const;
574 :
575 : /**
576 : * Copies the diagonal part of the matrix into \p dest.
577 : */
578 : virtual void get_diagonal (NumericVector<T> & dest) const = 0;
579 :
580 : /**
581 : * Copies the transpose of the matrix into \p dest, which may be
582 : * *this.
583 : */
584 : virtual void get_transpose (SparseMatrix<T> & dest) const = 0;
585 :
586 : /**
587 : * Get a row from the matrix
588 : * @param i The matrix row to get
589 : * @param indices A container that will be filled with the column indices
590 : * corresponding to (possibly) non-zero values
591 : * @param values A container holding the column values
592 : */
593 : virtual void get_row(numeric_index_type i,
594 : std::vector<numeric_index_type> & indices,
595 : std::vector<T> & values) const = 0;
596 :
597 : /**
598 : * Scales all elements of this matrix by \p scale
599 : */
600 : virtual void scale(const T scale);
601 :
602 : /**
603 : * @returns Whether the matrix supports hash table assembly
604 : */
605 0 : virtual bool supports_hash_table() const { return false; }
606 :
607 : /**
608 : * Sets whether to use hash table assembly. This will error if the passed-in value is true and the
609 : * matrix type does not support hash tables. Hash table or hash map assembly means storing maps
610 : * from i-j locations in the matrix to values. Because it is a hash map as opposed to a contiguous
611 : * array of data, no preallocation is required to use it
612 : */
613 : void use_hash_table(bool use_hash);
614 :
615 : /**
616 : * @returns Whether this matrix is using hash table assembly. Hash table or hash map assembly
617 : * means storing maps from i-j locations in the matrix to values. Because it is a hash map as
618 : * opposed to a contiguous array of data, no preallocation is required to use it
619 : */
620 32566 : bool use_hash_table() const { return _use_hash_table; }
621 :
622 : /**
623 : * Reset the memory storage of the matrix. Unlike \p clear(), this does not destroy the matrix but
624 : * rather will reset the matrix to use the original preallocation or when using hash table matrix
625 : * assembly (see \p use_hash_table()) will reset (clear) the hash table used for assembly. In the
626 : * words of the \p MatResetPreallocation documentation in PETSc, 'current values in the matrix are
627 : * lost in this call', so a user can expect to have back their original sparsity pattern in a
628 : * zeroed state
629 : */
630 0 : virtual void restore_original_nonzero_pattern() { libmesh_not_implemented(); }
631 :
632 : protected:
633 : /**
634 : * Protected implementation of the create_submatrix and reinit_submatrix
635 : * routines.
636 : *
637 : * \note This function must be overridden in derived classes for it
638 : * to work properly!
639 : */
640 0 : virtual void _get_submatrix(SparseMatrix<T> & /*submatrix*/,
641 : const std::vector<numeric_index_type> & /*rows*/,
642 : const std::vector<numeric_index_type> & /*cols*/,
643 : const bool /*reuse_submatrix*/) const
644 : {
645 0 : libmesh_not_implemented();
646 : }
647 :
648 : /**
649 : * The \p DofMap object associated with this object. May be queried
650 : * for degree-of-freedom counts on processors.
651 : */
652 : DofMap const * _dof_map;
653 :
654 : /**
655 : * The \p sparsity pattern associated with this object. Should be
656 : * queried for entry counts (or with need_full_sparsity_pattern,
657 : * patterns) when needed.
658 : */
659 : SparsityPattern::Build const * _sp;
660 :
661 : /**
662 : * Flag indicating whether or not the matrix has been initialized.
663 : */
664 : bool _is_initialized;
665 :
666 : /**
667 : * Flag indicating whether the matrix is assembled using a hash table
668 : */
669 : bool _use_hash_table;
670 : };
671 :
672 :
673 :
674 : //-----------------------------------------------------------------------
675 : // SparseMatrix inline members
676 : template <typename T>
677 : void
678 19374 : SparseMatrix<T>::use_hash_table(const bool use_hash)
679 : {
680 19374 : libmesh_error_msg_if(use_hash && !this->supports_hash_table(),
681 : "This matrix class does not support hash table assembly");
682 19374 : this->_use_hash_table = use_hash;
683 19374 : }
684 :
685 : // For SGI MIPSpro this implementation must occur after
686 : // the full specialization of the print() member.
687 : //
688 : // It's generally easier to define these friend functions in the header
689 : // file.
690 : template <typename T>
691 0 : std::ostream & operator << (std::ostream & os, const SparseMatrix<T> & m)
692 : {
693 0 : m.print(os);
694 0 : return os;
695 : }
696 :
697 : template <typename T>
698 : auto
699 : l1_norm(const SparseMatrix<T> & mat)
700 : {
701 : return mat.l1_norm();
702 : }
703 :
704 : template <typename T>
705 : auto
706 : l1_norm_diff(const SparseMatrix<T> & mat1, const SparseMatrix<T> & mat2)
707 : {
708 : return mat1.l1_norm_diff(mat2);
709 : }
710 :
711 : } // namespace libMesh
712 :
713 :
714 : #endif // LIBMESH_SPARSE_MATRIX_H
|
