Line data Source code
1 : // The libMesh Finite Element Library.
2 : // Copyright (C) 2002-2026 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_DENSE_MATRIX_H
21 : #define LIBMESH_DENSE_MATRIX_H
22 :
23 : // Local Includes
24 : #include "libmesh/libmesh_common.h"
25 : #include "libmesh/dense_matrix_base.h"
26 : #include "libmesh/int_range.h"
27 : #include "libmesh/parallel_algorithms.h"
28 :
29 : // For the definition of PetscBLASInt.
30 : #if (LIBMESH_HAVE_PETSC)
31 : # include "libmesh/petsc_macro.h"
32 : # ifdef I
33 : # define LIBMESH_SAW_I
34 : # endif
35 :
36 : #include "libmesh/ignore_warnings.h"
37 : # include <petscsys.h>
38 : #include "libmesh/restore_warnings.h"
39 :
40 : # ifndef LIBMESH_SAW_I
41 : # undef I // Avoid complex.h contamination
42 : # endif
43 : #endif
44 :
45 : // C++ includes
46 : #include <algorithm>
47 : #include <initializer_list>
48 : #include <vector>
49 :
50 : #ifdef LIBMESH_HAVE_METAPHYSICL
51 : #include "metaphysicl/dualnumber_decl.h"
52 : #include "metaphysicl/raw_type.h"
53 : #endif
54 :
55 : namespace libMesh
56 : {
57 :
58 : // Forward Declarations
59 : template <typename T> class DenseVector;
60 :
61 : /**
62 : * Defines a dense matrix for use in Finite Element-type computations.
63 : * Useful for storing element stiffness matrices before summation into
64 : * a global matrix. All overridden virtual functions are documented
65 : * in dense_matrix_base.h.
66 : *
67 : * \author Benjamin S. Kirk
68 : * \date 2002
69 : * \brief A matrix object used for finite element assembly and numerics.
70 : */
71 : template<typename T>
72 315200 : class DenseMatrix : public DenseMatrixBase<T>
73 : {
74 : public:
75 :
76 : /**
77 : * Constructor. Creates a dense matrix of dimension \p m by \p n.
78 : */
79 : DenseMatrix(const unsigned int new_m=0,
80 : const unsigned int new_n=0);
81 :
82 : /**
83 : * Constructor taking the number of rows, columns, and an
84 : * initializer_list, which must be of length nrow * ncol, of
85 : * row-major values to initialize the DenseMatrix with.
86 : */
87 : template <typename T2>
88 : DenseMatrix(unsigned int nrow,
89 : unsigned int ncol,
90 : std::initializer_list<T2> init_list);
91 :
92 : /**
93 : * The 5 special functions can be defaulted for this class, as it
94 : * does not manage any memory itself.
95 : */
96 0 : DenseMatrix (DenseMatrix &&) = default;
97 881148 : DenseMatrix (const DenseMatrix &) = default;
98 262326 : DenseMatrix & operator= (const DenseMatrix &) = default;
99 : DenseMatrix & operator= (DenseMatrix &&) = default;
100 104461424 : virtual ~DenseMatrix() = default;
101 :
102 : /**
103 : * Sets all elements of the matrix to 0 and resets any decomposition
104 : * flag which may have been previously set. This allows e.g. a new
105 : * LU decomposition to be computed while reusing the same storage.
106 : */
107 : virtual void zero() override final;
108 :
109 : /**
110 : * Get submatrix with the smallest row and column indices and the submatrix size.
111 : */
112 : DenseMatrix sub_matrix(unsigned int row_id, unsigned int row_size,
113 : unsigned int col_id, unsigned int col_size) const;
114 :
115 : /**
116 : * \returns The \p (i,j) element of the matrix.
117 : */
118 : T operator() (const unsigned int i,
119 : const unsigned int j) const;
120 :
121 : /**
122 : * \returns The \p (i,j) element of the matrix as a writable reference.
123 : */
124 : T & operator() (const unsigned int i,
125 : const unsigned int j);
126 :
127 54989648 : virtual T el(const unsigned int i,
128 : const unsigned int j) const override final
129 54989648 : { return (*this)(i,j); }
130 :
131 22183512 : virtual T & el(const unsigned int i,
132 : const unsigned int j) override final
133 22183512 : { return (*this)(i,j); }
134 :
135 : virtual void left_multiply (const DenseMatrixBase<T> & M2) override final;
136 :
137 : /**
138 : * Left multiplies by the matrix \p M2 of different type
139 : */
140 : template <typename T2>
141 : void left_multiply (const DenseMatrixBase<T2> & M2);
142 :
143 : virtual void right_multiply (const DenseMatrixBase<T> & M2) override final;
144 :
145 : /**
146 : * Right multiplies by the matrix \p M2 of different type
147 : */
148 : template <typename T2>
149 : void right_multiply (const DenseMatrixBase<T2> & M2);
150 :
151 : /**
152 : * Performs the matrix-vector multiplication,
153 : * \p dest := (*this) * \p arg.
154 : */
155 : void vector_mult (DenseVector<T> & dest,
156 : const DenseVector<T> & arg) const;
157 :
158 : /**
159 : * Performs the matrix-vector multiplication,
160 : * \p dest := (*this) * \p arg
161 : * on mixed types
162 : */
163 : template <typename T2>
164 : void vector_mult (DenseVector<typename CompareTypes<T,T2>::supertype> & dest,
165 : const DenseVector<T2> & arg) const;
166 :
167 : /**
168 : * Performs the matrix-vector multiplication,
169 : * \p dest := (*this)^T * \p arg.
170 : */
171 : void vector_mult_transpose (DenseVector<T> & dest,
172 : const DenseVector<T> & arg) const;
173 :
174 : /**
175 : * Performs the matrix-vector multiplication,
176 : * \p dest := (*this)^T * \p arg.
177 : * on mixed types
178 : */
179 : template <typename T2>
180 : void vector_mult_transpose (DenseVector<typename CompareTypes<T,T2>::supertype> & dest,
181 : const DenseVector<T2> & arg) const;
182 :
183 : /**
184 : * Performs the scaled matrix-vector multiplication,
185 : * \p dest += \p factor * (*this) * \p arg.
186 : */
187 : void vector_mult_add (DenseVector<T> & dest,
188 : const T factor,
189 : const DenseVector<T> & arg) const;
190 :
191 : /**
192 : * Performs the scaled matrix-vector multiplication,
193 : * \p dest += \p factor * (*this) * \p arg.
194 : * on mixed types
195 : */
196 : template <typename T2, typename T3>
197 : void vector_mult_add (DenseVector<typename CompareTypes<T, typename CompareTypes<T2,T3>::supertype>::supertype> & dest,
198 : const T2 factor,
199 : const DenseVector<T3> & arg) const;
200 :
201 : /**
202 : * Put the \p sub_m x \p sub_n principal submatrix into \p dest.
203 : */
204 : void get_principal_submatrix (unsigned int sub_m, unsigned int sub_n, DenseMatrix<T> & dest) const;
205 :
206 : /**
207 : * Put the \p sub_m x \p sub_m principal submatrix into \p dest.
208 : */
209 : void get_principal_submatrix (unsigned int sub_m, DenseMatrix<T> & dest) const;
210 :
211 : /**
212 : * Computes the outer (dyadic) product of two vectors and stores in (*this).
213 : *
214 : * The outer product of two real-valued vectors \f$\mathbf{a}\f$ and \f$\mathbf{b}\f$ is
215 : * \f[
216 : * (\mathbf{a}\mathbf{b}^T)_{i,j} = \mathbf{a}_i \mathbf{b}_j .
217 : * \f]
218 : * The outer product of two complex-valued vectors \f$\mathbf{a}\f$ and \f$\mathbf{b}\f$ is
219 : * \f[
220 : * (\mathbf{a}\mathbf{b}^H)_{i,j} = \mathbf{a}_i \mathbf{b}^*_j ,
221 : * \f]
222 : * where \f$H\f$ denotes the conjugate transpose of the vector and \f$*\f$
223 : * denotes the complex conjugate.
224 : *
225 : * \param[in] a Vector whose entries correspond to rows in the product matrix.
226 : * \param[in] b Vector whose entries correspond to columns in the product matrix.
227 : */
228 : void outer_product(const DenseVector<T> & a, const DenseVector<T> & b);
229 :
230 : /**
231 : * Assignment-from-other-matrix-type operator.
232 : *
233 : * Copies the dense matrix of type T2 into the present matrix. This
234 : * is useful for copying real matrices into complex ones for further
235 : * operations.
236 : *
237 : * \returns A reference to *this.
238 : */
239 : template <typename T2>
240 : DenseMatrix<T> & operator = (const DenseMatrix<T2> & other_matrix);
241 :
242 : /**
243 : * STL-like swap method
244 : */
245 : void swap(DenseMatrix<T> & other_matrix);
246 :
247 : /**
248 : * Resizes the matrix to the specified size and calls zero(). Will
249 : * never free memory, but may allocate more. Note: when the matrix
250 : * is zero()'d, any decomposition (LU, Cholesky, etc.) is also
251 : * cleared, forcing a new decomposition to be computed the next time
252 : * e.g. lu_solve() is called.
253 : */
254 : void resize(const unsigned int new_m,
255 : const unsigned int new_n);
256 :
257 : /**
258 : * Multiplies every element in the matrix by \p factor.
259 : */
260 : void scale (const T factor);
261 :
262 : /**
263 : * Multiplies every element in the column \p col matrix by \p factor.
264 : */
265 : void scale_column (const unsigned int col, const T factor);
266 :
267 : /**
268 : * Multiplies every element in the matrix by \p factor.
269 : *
270 : * \returns A reference to *this.
271 : */
272 : DenseMatrix<T> & operator *= (const T factor);
273 :
274 : /**
275 : * Adds \p factor times \p mat to this matrix.
276 : *
277 : * \returns A reference to *this.
278 : */
279 : template<typename T2, typename T3>
280 : typename std::enable_if<
281 : ScalarTraits<T2>::value, void >::type add (const T2 factor,
282 : const DenseMatrix<T3> & mat);
283 :
284 : /**
285 : * \returns \p true if \p mat is exactly equal to this matrix, \p false otherwise.
286 : */
287 : bool operator== (const DenseMatrix<T> & mat) const;
288 :
289 : /**
290 : * \returns \p true if \p mat is not exactly equal to this matrix, false otherwise.
291 : */
292 : bool operator!= (const DenseMatrix<T> & mat) const;
293 :
294 : /**
295 : * Adds \p mat to this matrix.
296 : *
297 : * \returns A reference to *this.
298 : */
299 : DenseMatrix<T> & operator+= (const DenseMatrix<T> & mat);
300 :
301 : /**
302 : * Subtracts \p mat from this matrix.
303 : *
304 : * \returns A reference to *this.
305 : */
306 : DenseMatrix<T> & operator-= (const DenseMatrix<T> & mat);
307 :
308 : /**
309 : * \returns The minimum entry in the matrix, or the minimum real
310 : * part in the case of complex numbers.
311 : */
312 : auto min () const -> decltype(libmesh_real(T(0)));
313 :
314 : /**
315 : * \returns The maximum entry in the matrix, or the maximum real
316 : * part in the case of complex numbers.
317 : */
318 : auto max () const -> decltype(libmesh_real(T(0)));
319 :
320 : /**
321 : * \returns The l1-norm of the matrix, that is, the max column sum:
322 : *
323 : * \f$ |M|_1 = max_{all columns j} \sum_{all rows i} |M_ij| \f$,
324 : *
325 : * This is the natural matrix norm that is compatible to the l1-norm
326 : * for vectors, i.e. \f$ |Mv|_1 \leq |M|_1 |v|_1 \f$.
327 : */
328 : auto l1_norm () const;
329 :
330 : /**
331 : * \returns The linfty-norm of the matrix, that is, the max row sum:
332 : *
333 : * \f$ |M|_\infty = max_{all rows i} \sum_{all columns j} |M_ij| \f$,
334 : *
335 : * This is the natural matrix norm that is compatible to the
336 : * linfty-norm of vectors, i.e. \f$ |Mv|_\infty \leq |M|_\infty |v|_\infty \f$.
337 : */
338 : auto linfty_norm () const;
339 :
340 : /**
341 : * Left multiplies by the transpose of the matrix \p A.
342 : */
343 : void left_multiply_transpose (const DenseMatrix<T> & A);
344 :
345 : /**
346 : * Left multiplies by the transpose of the matrix \p A which
347 : * contains a different numerical type.
348 : */
349 : template <typename T2>
350 : void left_multiply_transpose (const DenseMatrix<T2> & A);
351 :
352 :
353 : /**
354 : * Right multiplies by the transpose of the matrix \p A
355 : */
356 : void right_multiply_transpose (const DenseMatrix<T> & A);
357 :
358 : /**
359 : * Right multiplies by the transpose of the matrix \p A which
360 : * contains a different numerical type.
361 : */
362 : template <typename T2>
363 : void right_multiply_transpose (const DenseMatrix<T2> & A);
364 :
365 : /**
366 : * \returns The \p (i,j) element of the transposed matrix.
367 : */
368 : T transpose (const unsigned int i,
369 : const unsigned int j) const;
370 :
371 : /**
372 : * Put the tranposed matrix into \p dest.
373 : */
374 : void get_transpose(DenseMatrix<T> & dest) const;
375 :
376 : /**
377 : * \returns A reference to the underlying data storage vector.
378 : *
379 : * This should be used with caution (i.e. one should not change the
380 : * size of the vector, etc.) but is useful for interoperating with
381 : * low level BLAS routines which expect a simple array.
382 : */
383 2326073 : std::vector<T> & get_values() { return _val; }
384 :
385 : /**
386 : * \returns A constant reference to the underlying data storage vector.
387 : */
388 3300132 : const std::vector<T> & get_values() const { return _val; }
389 :
390 : /**
391 : * Condense-out the \p (i,j) entry of the matrix, forcing
392 : * it to take on the value \p val. This is useful in numerical
393 : * simulations for applying boundary conditions. Preserves the
394 : * symmetry of the matrix.
395 : */
396 0 : void condense(const unsigned int i,
397 : const unsigned int j,
398 : const T val,
399 : DenseVector<T> & rhs)
400 0 : { DenseMatrixBase<T>::condense (i, j, val, rhs); }
401 :
402 : /**
403 : * Solve the system Ax=b given the input vector b. Partial pivoting
404 : * is performed by default in order to keep the algorithm stable to
405 : * the effects of round-off error.
406 : *
407 : * Important note: once you call lu_solve(), you must _not_ modify
408 : * the entries of the matrix via calls to operator(i,j) and call
409 : * lu_solve() again without first calling either zero() or resize(),
410 : * otherwise the code will skip computing the decomposition of the
411 : * matrix and go directly to the back substitution step. This is
412 : * done on purpose for efficiency, so that the same LU decomposition
413 : * can be used with multiple right-hand sides, but it does also make
414 : * it possible to "shoot yourself in the foot", so be careful!
415 : */
416 : void lu_solve (const DenseVector<T> & b,
417 : DenseVector<T> & x);
418 :
419 : /**
420 : * For symmetric positive definite (SPD) matrices. A Cholesky factorization
421 : * of A such that A = L L^T is about twice as fast as a standard LU
422 : * factorization. Therefore you can use this method if you know a-priori
423 : * that the matrix is SPD. If the matrix is not SPD, an error is generated.
424 : * One nice property of Cholesky decompositions is that they do not require
425 : * pivoting for stability.
426 : *
427 : * Important note: once you call cholesky_solve(), you must _not_
428 : * modify the entries of the matrix via calls to operator(i,j) and
429 : * call cholesky_solve() again without first calling either zero()
430 : * or resize(), otherwise the code will skip computing the
431 : * decomposition of the matrix and go directly to the back
432 : * substitution step. This is done on purpose for efficiency, so
433 : * that the same decomposition can be used with multiple right-hand
434 : * sides, but it does also make it possible to "shoot yourself in
435 : * the foot", so be careful!
436 : *
437 : * \note This method may also be used when A is real-valued and x
438 : * and b are complex-valued.
439 : */
440 : template <typename T2>
441 : void cholesky_solve(const DenseVector<T2> & b,
442 : DenseVector<T2> & x);
443 :
444 : /**
445 : * Compute the singular value decomposition of the matrix.
446 : * On exit, sigma holds all of the singular values (in
447 : * descending order).
448 : *
449 : * The implementation uses PETSc's interface to BLAS/LAPACK.
450 : * If this is not available, this function throws an error.
451 : */
452 : void svd(DenseVector<Real> & sigma);
453 :
454 : /**
455 : * Compute the "reduced" singular value decomposition of the matrix.
456 : * On exit, sigma holds all of the singular values (in
457 : * descending order), U holds the left singular vectors,
458 : * and VT holds the transpose of the right singular vectors.
459 : * In the reduced SVD, U has min(m,n) columns and VT has
460 : * min(m,n) rows. (In the "full" SVD, U and VT would be square.)
461 : *
462 : * The implementation uses PETSc's interface to BLAS/LAPACK.
463 : * If this is not available, this function throws an error.
464 : */
465 : void svd(DenseVector<Real> & sigma,
466 : DenseMatrix<Number> & U,
467 : DenseMatrix<Number> & VT);
468 :
469 : /**
470 : * Solve the system of equations \f$ A x = rhs \f$ for \f$ x \f$ in the
471 : * least-squares sense. \f$ A \f$ may be non-square and/or rank-deficient.
472 : * You can control which singular values are treated as zero by
473 : * changing the "rcond" parameter. Singular values S(i) for which
474 : * S(i) <= rcond*S(1) are treated as zero for purposes of the solve.
475 : * Passing a negative number for rcond forces a "machine precision"
476 : * value to be used instead.
477 : *
478 : * This function is marked const, since due to various
479 : * implementation details, we do not need to modify the contents of
480 : * A in order to compute the SVD (a copy is made internally
481 : * instead).
482 : *
483 : * Requires PETSc >= 3.1 since this was the first version to provide
484 : * the LAPACKgelss_ wrapper.
485 : */
486 : void svd_solve(const DenseVector<T> & rhs,
487 : DenseVector<T> & x,
488 : Real rcond=std::numeric_limits<Real>::epsilon()) const;
489 :
490 : /**
491 : * Compute the eigenvalues (both real and imaginary parts) of a general matrix.
492 : *
493 : * Warning: the contents of \p *this are overwritten by this function!
494 : *
495 : * The implementation requires the LAPACKgeev_ function which is wrapped by PETSc.
496 : */
497 : void evd(DenseVector<T> & lambda_real,
498 : DenseVector<T> & lambda_imag);
499 :
500 : /**
501 : * Compute the eigenvalues (both real and imaginary parts) and left
502 : * eigenvectors of a general matrix, \f$ A \f$.
503 : *
504 : * Warning: the contents of \p *this are overwritten by this function!
505 : *
506 : * The left eigenvector \f$ u_j \f$ of \f$ A \f$ satisfies:
507 : * \f$ u_j^H A = lambda_j u_j^H \f$
508 : * where \f$ u_j^H \f$ denotes the conjugate-transpose of \f$ u_j \f$.
509 : *
510 : * If the j-th and (j+1)-st eigenvalues form a complex conjugate
511 : * pair, then the j-th and (j+1)-st columns of VL "share" their
512 : * real-valued storage in the following way:
513 : * u_j = VL(:,j) + i*VL(:,j+1) and
514 : * u_{j+1} = VL(:,j) - i*VL(:,j+1).
515 : *
516 : * The implementation requires the LAPACKgeev_ routine which is provided by PETSc.
517 : */
518 : void evd_left(DenseVector<T> & lambda_real,
519 : DenseVector<T> & lambda_imag,
520 : DenseMatrix<T> & VL);
521 :
522 : /**
523 : * Compute the eigenvalues (both real and imaginary parts) and right
524 : * eigenvectors of a general matrix, \f$ A \f$.
525 : *
526 : * Warning: the contents of \p *this are overwritten by this function!
527 : *
528 : * The right eigenvector \f$ v_j \f$ of \f$ A \f$ satisfies:
529 : * \f$ A v_j = lambda_j v_j \f$
530 : * where \f$ lambda_j \f$ is its corresponding eigenvalue.
531 : *
532 : * \note If the j-th and (j+1)-st eigenvalues form a complex
533 : * conjugate pair, then the j-th and (j+1)-st columns of VR "share"
534 : * their real-valued storage in the following way:
535 : * v_j = VR(:,j) + i*VR(:,j+1) and
536 : * v_{j+1} = VR(:,j) - i*VR(:,j+1).
537 : *
538 : * The implementation requires the LAPACKgeev_ routine which is provided by PETSc.
539 : */
540 : void evd_right(DenseVector<T> & lambda_real,
541 : DenseVector<T> & lambda_imag,
542 : DenseMatrix<T> & VR);
543 :
544 : /**
545 : * Compute the eigenvalues (both real and imaginary parts) as well as the left
546 : * and right eigenvectors of a general matrix.
547 : *
548 : * Warning: the contents of \p *this are overwritten by this function!
549 : *
550 : * See the documentation of the \p evd_left() and \p evd_right()
551 : * functions for more information. The implementation requires the
552 : * LAPACKgeev_ routine which is provided by PETSc.
553 : */
554 : void evd_left_and_right(DenseVector<T> & lambda_real,
555 : DenseVector<T> & lambda_imag,
556 : DenseMatrix<T> & VL,
557 : DenseMatrix<T> & VR);
558 :
559 : /**
560 : * \returns The determinant of the matrix.
561 : *
562 : * \note Implemented by computing an LU decomposition and then
563 : * taking the product of the diagonal terms. Therefore this is a
564 : * non-const method which modifies the entries of the matrix.
565 : */
566 : T det();
567 :
568 : /**
569 : * Returns true iff every entry is finite.
570 : */
571 : friend bool isfinite (const DenseMatrix<T> & var)
572 : {
573 : using std::isfinite;
574 : using libMesh::isfinite; // for T==complex
575 : for (const T & v : var._val)
576 : if (!isfinite(v))
577 : return false;
578 : return true;
579 : }
580 :
581 : /**
582 : * Returns true iff no entry is NaN and any entry is infinite.
583 : *
584 : * This is arguably inconsistent with our std::complex overload (and
585 : * the C99 Annex G recommendations for _Complex, and C++
586 : * std::complex arithmetic), which treats mixed (inf,NaN) pairs as
587 : * infinite, but this is probably safer for users.
588 : */
589 : friend bool isinf (const DenseMatrix<T> & var)
590 : {
591 : using std::isinf;
592 : using libMesh::isinf; // for T==complex
593 : using std::isnan;
594 : using libMesh::isnan;
595 : bool has_inf = false;
596 : for (const T & v : var._val)
597 : {
598 : // NaN anywhere makes us NaN, not inf
599 : if (isnan(v))
600 : return false;
601 : has_inf = has_inf || isinf(v);
602 : }
603 : return has_inf;
604 : }
605 :
606 : /**
607 : * Returns true iff any entry is NaN.
608 : *
609 : * This is arguably inconsistent with our std::complex overload (and
610 : * the C99 Annex G recommendations for _Complex, and C++
611 : * std::complex arithmetic), which treats mixed (inf,NaN) pairs as
612 : * infinite, but this is probably safer for users.
613 : */
614 : friend bool isnan (const DenseMatrix<T> & var)
615 : {
616 : using std::isnan;
617 : using libMesh::isnan; // for T==complex
618 : for (const T & v : var._val)
619 : if (isnan(v))
620 : return true;
621 : return false;
622 : }
623 :
624 : /**
625 : * Computes the inverse of the dense matrix (assuming it is invertible)
626 : * by first computing the LU decomposition and then performing multiple
627 : * back substitution steps. Follows the algorithm from Numerical Recipes
628 : * in C that is available on the web.
629 : *
630 : * This routine is commented out since it is not really a memory- or
631 : * computationally- efficient implementation. Also, you typically
632 : * don't need the actual inverse for anything, and can use something
633 : * like lu_solve() instead.
634 : */
635 : // void inverse();
636 :
637 : /**
638 : * Run-time selectable option to turn on/off BLAS support.
639 : * This was primarily used for testing purposes, and could be
640 : * removed...
641 : */
642 : bool use_blas_lapack;
643 :
644 : /**
645 : * Helper structure for determining whether to use blas_lapack
646 : */
647 : struct UseBlasLapack
648 : {
649 : static const bool value = false;
650 : };
651 :
652 : private:
653 :
654 : /**
655 : * The actual data values, stored as a 1D array.
656 : */
657 : std::vector<T> _val;
658 :
659 : /**
660 : * Form the LU decomposition of the matrix. This function
661 : * is private since it is only called as part of the implementation
662 : * of the lu_solve(...) function.
663 : */
664 : void _lu_decompose ();
665 :
666 : /**
667 : * Solves the system Ax=b through back substitution. This function
668 : * is private since it is only called as part of the implementation
669 : * of the lu_solve(...) function.
670 : */
671 : void _lu_back_substitute (const DenseVector<T> & b,
672 : DenseVector<T> & x) const;
673 :
674 : /**
675 : * Decomposes a symmetric positive definite matrix into a
676 : * product of two lower triangular matrices according to
677 : * A = LL^T.
678 : *
679 : * \note This program generates an error if the matrix is not SPD.
680 : */
681 : void _cholesky_decompose();
682 :
683 : /**
684 : * Solves the equation Ax=b for the unknown value x and rhs
685 : * b based on the Cholesky factorization of A.
686 : *
687 : * \note This method may be used when A is real-valued and b and x
688 : * are complex-valued.
689 : */
690 : template <typename T2>
691 : void _cholesky_back_substitute(const DenseVector<T2> & b,
692 : DenseVector<T2> & x) const;
693 :
694 : /**
695 : * The decomposition schemes above change the entries of the matrix
696 : * A. It is therefore an error to call A.lu_solve() and subsequently
697 : * call A.cholesky_solve() since the result will probably not match
698 : * any desired outcome. This typedef keeps track of which decomposition
699 : * has been called for this matrix.
700 : */
701 : enum DecompositionType {LU=0, CHOLESKY=1, LU_BLAS_LAPACK, NONE};
702 :
703 : /**
704 : * This flag keeps track of which type of decomposition has been
705 : * performed on the matrix.
706 : */
707 : DecompositionType _decomposition_type;
708 :
709 : /**
710 : * Enumeration used to determine the behavior of the _multiply_blas
711 : * function.
712 : */
713 : enum _BLAS_Multiply_Flag {
714 : LEFT_MULTIPLY = 0,
715 : RIGHT_MULTIPLY,
716 : LEFT_MULTIPLY_TRANSPOSE,
717 : RIGHT_MULTIPLY_TRANSPOSE
718 : };
719 :
720 : /**
721 : * The _multiply_blas function computes A <- op(A) * op(B) using
722 : * BLAS gemm function. Used in the right_multiply(),
723 : * left_multiply(), right_multiply_transpose(), and
724 : * left_multiply_transpose() routines.
725 : * [ Implementation in dense_matrix_blas_lapack.C ]
726 : */
727 : void _multiply_blas(const DenseMatrixBase<T> & other,
728 : _BLAS_Multiply_Flag flag);
729 :
730 : /**
731 : * Computes an LU factorization of the matrix using the
732 : * Lapack routine "getrf". This routine should only be
733 : * used by the "use_blas_lapack" branch of the lu_solve()
734 : * function. After the call to this function, the matrix
735 : * is replaced by its factorized version, and the
736 : * DecompositionType is set to LU_BLAS_LAPACK.
737 : * [ Implementation in dense_matrix_blas_lapack.C ]
738 : */
739 : void _lu_decompose_lapack();
740 :
741 : /**
742 : * Computes an SVD of the matrix using the
743 : * Lapack routine "getsvd".
744 : * [ Implementation in dense_matrix_blas_lapack.C ]
745 : */
746 : void _svd_lapack(DenseVector<Real> & sigma);
747 :
748 : /**
749 : * Computes a "reduced" SVD of the matrix using the
750 : * Lapack routine "getsvd".
751 : * [ Implementation in dense_matrix_blas_lapack.C ]
752 : */
753 : void _svd_lapack(DenseVector<Real> & sigma,
754 : DenseMatrix<Number> & U,
755 : DenseMatrix<Number> & VT);
756 :
757 : /**
758 : * Called by svd_solve(rhs).
759 : */
760 : void _svd_solve_lapack(const DenseVector<T> & rhs,
761 : DenseVector<T> & x,
762 : Real rcond) const;
763 :
764 : /**
765 : * Helper function that actually performs the SVD.
766 : * [ Implementation in dense_matrix_blas_lapack.C ]
767 : */
768 : void _svd_helper (char JOBU,
769 : char JOBVT,
770 : std::vector<Real> & sigma_val,
771 : std::vector<Number> & U_val,
772 : std::vector<Number> & VT_val);
773 :
774 : /**
775 : * Computes the eigenvalues of the matrix using the Lapack routine
776 : * "DGEEV". If VR and/or VL are not nullptr, then the matrix of right
777 : * and/or left eigenvectors is also computed and returned by this
778 : * function.
779 : *
780 : * [ Implementation in dense_matrix_blas_lapack.C ]
781 : */
782 : void _evd_lapack(DenseVector<T> & lambda_real,
783 : DenseVector<T> & lambda_imag,
784 : DenseMatrix<T> * VL = nullptr,
785 : DenseMatrix<T> * VR = nullptr);
786 :
787 : /**
788 : * Array used to store pivot indices. May be used by whatever
789 : * factorization is currently active, clients of the class should
790 : * not rely on it for any reason.
791 : */
792 : #if (LIBMESH_HAVE_PETSC && LIBMESH_USE_REAL_NUMBERS)
793 : typedef PetscBLASInt pivot_index_t;
794 : #else
795 : typedef int pivot_index_t;
796 : #endif
797 : std::vector<pivot_index_t> _pivots;
798 :
799 : /**
800 : * Companion function to _lu_decompose_lapack(). Do not use
801 : * directly, called through the public lu_solve() interface.
802 : * This function is logically const in that it does not modify
803 : * the matrix, but since we are just calling LAPACK routines,
804 : * it's less const_cast hassle to just declare the function
805 : * non-const.
806 : * [ Implementation in dense_matrix_blas_lapack.C ]
807 : */
808 : void _lu_back_substitute_lapack (const DenseVector<T> & b,
809 : DenseVector<T> & x);
810 :
811 : /**
812 : * Uses the BLAS GEMV function (through PETSc) to compute
813 : *
814 : * dest := alpha*A*arg + beta*dest
815 : *
816 : * where alpha and beta are scalars, A is this matrix, and
817 : * arg and dest are input vectors of appropriate size. If
818 : * trans is true, the transpose matvec is computed instead.
819 : * By default, trans==false.
820 : *
821 : * [ Implementation in dense_matrix_blas_lapack.C ]
822 : */
823 : void _matvec_blas(T alpha, T beta,
824 : DenseVector<T> & dest,
825 : const DenseVector<T> & arg,
826 : bool trans=false) const;
827 :
828 : /**
829 : * Left multiplies by the transpose of the matrix \p A which
830 : * may contain a different numerical type.
831 : */
832 : template <typename T2>
833 : void _left_multiply_transpose (const DenseMatrix<T2> & A);
834 :
835 : /**
836 : * Right multiplies by the transpose of the matrix \p A which
837 : * may contain a different numerical type.
838 : */
839 : template <typename T2>
840 : void _right_multiply_transpose (const DenseMatrix<T2> & A);
841 : };
842 :
843 :
844 :
845 :
846 :
847 : // ------------------------------------------------------------
848 : /**
849 : * Provide Typedefs for dense matrices
850 : */
851 : namespace DenseMatrices
852 : {
853 :
854 : /**
855 : * Convenient definition of a real-only
856 : * dense matrix.
857 : */
858 : typedef DenseMatrix<Real> RealDenseMatrix;
859 :
860 : /**
861 : * This typedef may be either a real-only matrix, or a truly complex
862 : * matrix, depending on how \p Number was defined in \p
863 : * libmesh_common.h. Also, be aware of the fact that \p
864 : * DenseMatrix<T> is likely to be more efficient for real than for
865 : * complex data.
866 : */
867 : typedef DenseMatrix<Complex> ComplexDenseMatrix;
868 :
869 : }
870 :
871 :
872 :
873 : using namespace DenseMatrices;
874 :
875 : // The PETSc Lapack wrappers are only for PetscScalar, therefore we
876 : // can't e.g. get a Lapack version of DenseMatrix<Real>::lu_solve()
877 : // when libmesh/PETSc are compiled with complex numbers.
878 : #if defined(LIBMESH_HAVE_PETSC) && \
879 : defined(LIBMESH_USE_REAL_NUMBERS) && \
880 : defined(LIBMESH_DEFAULT_DOUBLE_PRECISION)
881 : template <>
882 : struct DenseMatrix<double>::UseBlasLapack
883 : {
884 : static const bool value = true;
885 : };
886 : #endif
887 :
888 :
889 : // ------------------------------------------------------------
890 : // Dense Matrix member functions
891 : template<typename T>
892 : inline
893 108454285 : DenseMatrix<T>::DenseMatrix(const unsigned int new_m,
894 : const unsigned int new_n) :
895 : DenseMatrixBase<T>(new_m,new_n),
896 92818226 : use_blas_lapack(DenseMatrix<T>::UseBlasLapack::value),
897 : _val(),
898 119191868 : _decomposition_type(NONE)
899 : {
900 97716702 : this->resize(new_m,new_n);
901 108454285 : }
902 :
903 : template <typename T>
904 : template <typename T2>
905 : DenseMatrix<T>::DenseMatrix(unsigned int nrow,
906 : unsigned int ncol,
907 : std::initializer_list<T2> init_list) :
908 : DenseMatrixBase<T>(nrow, ncol),
909 : use_blas_lapack(DenseMatrix<T>::UseBlasLapack::value),
910 : _val(init_list.begin(), init_list.end()),
911 : _decomposition_type(NONE)
912 : {
913 : // Make sure the user passed us an amount of data which is
914 : // consistent with the size of the matrix.
915 : libmesh_assert_equal_to(nrow * ncol, init_list.size());
916 : }
917 :
918 :
919 :
920 : template<typename T>
921 : inline
922 5991596 : void DenseMatrix<T>::swap(DenseMatrix<T> & other_matrix)
923 : {
924 : using std::swap;
925 594836 : swap(this->_m, other_matrix._m);
926 594836 : swap(this->_n, other_matrix._n);
927 594836 : _val.swap(other_matrix._val);
928 6586432 : DecompositionType _temp = _decomposition_type;
929 6586432 : _decomposition_type = other_matrix._decomposition_type;
930 6586432 : other_matrix._decomposition_type = _temp;
931 5991596 : }
932 :
933 :
934 : template <typename T>
935 : template <typename T2>
936 : inline
937 : DenseMatrix<T> &
938 : DenseMatrix<T>::operator=(const DenseMatrix<T2> & mat)
939 : {
940 : unsigned int mat_m = mat.m(), mat_n = mat.n();
941 : this->resize(mat_m, mat_n);
942 : for (unsigned int i=0; i<mat_m; i++)
943 : for (unsigned int j=0; j<mat_n; j++)
944 : (*this)(i,j) = mat(i,j);
945 :
946 : return *this;
947 : }
948 :
949 :
950 :
951 : template<typename T>
952 : inline
953 138896310 : void DenseMatrix<T>::resize(const unsigned int new_m,
954 : const unsigned int new_n)
955 : {
956 153650441 : _val.resize(new_m*new_n);
957 :
958 153650441 : this->_m = new_m;
959 153650441 : this->_n = new_n;
960 :
961 : // zero and set decomposition_type to NONE
962 8458623 : this->zero();
963 138896310 : }
964 :
965 :
966 :
967 : template<typename T>
968 : inline
969 8672739 : void DenseMatrix<T>::zero()
970 : {
971 154332520 : _decomposition_type = NONE;
972 :
973 8672739 : std::fill (_val.begin(), _val.end(), static_cast<T>(0));
974 8672739 : }
975 :
976 :
977 :
978 : template<typename T>
979 : inline
980 0 : DenseMatrix<T> DenseMatrix<T>::sub_matrix(unsigned int row_id, unsigned int row_size,
981 : unsigned int col_id, unsigned int col_size) const
982 : {
983 0 : libmesh_assert_less (row_id + row_size - 1, this->_m);
984 0 : libmesh_assert_less (col_id + col_size - 1, this->_n);
985 :
986 0 : DenseMatrix<T> sub;
987 0 : sub._m = row_size;
988 0 : sub._n = col_size;
989 0 : sub._val.resize(row_size * col_size);
990 :
991 0 : unsigned int end_col = this->_n - col_size - col_id;
992 0 : unsigned int p = row_id * this->_n;
993 0 : unsigned int q = 0;
994 0 : for (unsigned int i=0; i<row_size; i++)
995 : {
996 : // skip the beginning columns
997 0 : p += col_id;
998 0 : for (unsigned int j=0; j<col_size; j++)
999 0 : sub._val[q++] = _val[p++];
1000 : // skip the rest columns
1001 0 : p += end_col;
1002 : }
1003 :
1004 0 : return sub;
1005 0 : }
1006 :
1007 :
1008 :
1009 : template<typename T>
1010 : inline
1011 7902664 : T DenseMatrix<T>::operator () (const unsigned int i,
1012 : const unsigned int j) const
1013 : {
1014 7902664 : libmesh_assert_less (i*j, _val.size());
1015 7902664 : libmesh_assert_less (i, this->_m);
1016 7902664 : libmesh_assert_less (j, this->_n);
1017 :
1018 :
1019 : // return _val[(i) + (this->_m)*(j)]; // col-major
1020 884595700 : return _val[(i)*(this->_n) + (j)]; // row-major
1021 : }
1022 :
1023 :
1024 :
1025 : template<typename T>
1026 : inline
1027 0 : T & DenseMatrix<T>::operator () (const unsigned int i,
1028 : const unsigned int j)
1029 : {
1030 0 : libmesh_assert_less (i*j, _val.size());
1031 0 : libmesh_assert_less (i, this->_m);
1032 0 : libmesh_assert_less (j, this->_n);
1033 :
1034 : //return _val[(i) + (this->_m)*(j)]; // col-major
1035 8563029181 : return _val[(i)*(this->_n) + (j)]; // row-major
1036 : }
1037 :
1038 :
1039 :
1040 :
1041 :
1042 : template<typename T>
1043 : inline
1044 122111 : void DenseMatrix<T>::scale (const T factor)
1045 : {
1046 2005351276 : for (auto & v : _val)
1047 1904676072 : v *= factor;
1048 122111 : }
1049 :
1050 :
1051 : template<typename T>
1052 : inline
1053 0 : void DenseMatrix<T>::scale_column (const unsigned int col, const T factor)
1054 : {
1055 0 : for (auto i : make_range(this->m()))
1056 0 : (*this)(i, col) *= factor;
1057 0 : }
1058 :
1059 :
1060 :
1061 : template<typename T>
1062 : inline
1063 128111 : DenseMatrix<T> & DenseMatrix<T>::operator *= (const T factor)
1064 : {
1065 128111 : this->scale(factor);
1066 128111 : return *this;
1067 : }
1068 :
1069 :
1070 :
1071 : template<typename T>
1072 : template<typename T2, typename T3>
1073 : inline
1074 : typename std::enable_if<
1075 : ScalarTraits<T2>::value, void >::type
1076 1353372 : DenseMatrix<T>::add (const T2 factor,
1077 : const DenseMatrix<T3> & mat)
1078 : {
1079 114036 : libmesh_assert_equal_to (this->m(), mat.m());
1080 114036 : libmesh_assert_equal_to (this->n(), mat.n());
1081 :
1082 11507940 : for (auto i : make_range(this->m()))
1083 117040026 : for (auto j : make_range(this->n()))
1084 109118208 : (*this)(i,j) += factor * mat(i,j);
1085 1353372 : }
1086 :
1087 :
1088 :
1089 : template<typename T>
1090 : inline
1091 547008 : bool DenseMatrix<T>::operator == (const DenseMatrix<T> & mat) const
1092 : {
1093 91652784 : for (auto i : index_range(_val))
1094 91105776 : if (_val[i] != mat._val[i])
1095 0 : return false;
1096 :
1097 547008 : return true;
1098 : }
1099 :
1100 :
1101 :
1102 : template<typename T>
1103 : inline
1104 0 : bool DenseMatrix<T>::operator != (const DenseMatrix<T> & mat) const
1105 : {
1106 0 : for (auto i : index_range(_val))
1107 0 : if (_val[i] != mat._val[i])
1108 0 : return true;
1109 :
1110 0 : return false;
1111 : }
1112 :
1113 :
1114 :
1115 : template<typename T>
1116 : inline
1117 6263411 : DenseMatrix<T> & DenseMatrix<T>::operator += (const DenseMatrix<T> & mat)
1118 : {
1119 1899817523 : for (auto i : index_range(_val))
1120 2119545086 : _val[i] += mat._val[i];
1121 :
1122 6263411 : return *this;
1123 : }
1124 :
1125 :
1126 :
1127 : template<typename T>
1128 : inline
1129 0 : DenseMatrix<T> & DenseMatrix<T>::operator -= (const DenseMatrix<T> & mat)
1130 : {
1131 0 : for (auto i : index_range(_val))
1132 0 : _val[i] -= mat._val[i];
1133 :
1134 0 : return *this;
1135 : }
1136 :
1137 :
1138 :
1139 : template<typename T>
1140 : inline
1141 0 : auto DenseMatrix<T>::min () const -> decltype(libmesh_real(T(0)))
1142 : {
1143 0 : libmesh_assert (this->_m);
1144 0 : libmesh_assert (this->_n);
1145 : typedef decltype(libmesh_real(T(0))) realfromT;
1146 : return libmesh_transform_reduce
1147 0 : (_val.begin(), _val.end(), std::numeric_limits<realfromT>::max(),
1148 0 : [](const auto & a, const auto & b){using std::min; return min(a,b);},
1149 0 : [](const T & v){return libmesh_real(v);});
1150 : }
1151 :
1152 :
1153 :
1154 : template<typename T>
1155 : inline
1156 0 : auto DenseMatrix<T>::max () const -> decltype(libmesh_real(T(0)))
1157 : {
1158 0 : libmesh_assert (this->_m);
1159 0 : libmesh_assert (this->_n);
1160 : typedef decltype(libmesh_real(T(0))) realfromT;
1161 : return libmesh_transform_reduce
1162 0 : (_val.begin(), _val.end(), std::numeric_limits<realfromT>::lowest(),
1163 0 : [](const auto & a, const auto & b){using std::max; return max(a,b);},
1164 0 : [](const T & v){return libmesh_real(v);});
1165 : }
1166 :
1167 :
1168 :
1169 : template<typename T>
1170 : inline
1171 27228 : auto DenseMatrix<T>::l1_norm () const
1172 : {
1173 27228 : libmesh_assert (this->_m);
1174 27228 : libmesh_assert (this->_n);
1175 :
1176 : using std::abs;
1177 27228 : auto columnsum = abs(T(0));
1178 300732 : for (unsigned int i=0; i!=this->_m; i++)
1179 : {
1180 273504 : columnsum += abs((*this)(i,0));
1181 : }
1182 27228 : auto my_max = columnsum;
1183 273504 : for (unsigned int j=1; j!=this->_n; j++)
1184 : {
1185 246276 : columnsum = 0.;
1186 3149772 : for (unsigned int i=0; i!=this->_m; i++)
1187 : {
1188 2903496 : columnsum += abs((*this)(i,j));
1189 : }
1190 246276 : my_max = (my_max > columnsum? my_max : columnsum);
1191 : }
1192 27228 : return my_max;
1193 : }
1194 :
1195 :
1196 :
1197 : template<typename T>
1198 : inline
1199 0 : auto DenseMatrix<T>::linfty_norm () const
1200 : {
1201 0 : libmesh_assert (this->_m);
1202 0 : libmesh_assert (this->_n);
1203 : using std::abs;
1204 :
1205 0 : auto rowsum = abs(T(0));
1206 0 : for (unsigned int j=0; j!=this->_n; j++)
1207 : {
1208 0 : rowsum += abs((*this)(0,j));
1209 : }
1210 0 : auto my_max = rowsum;
1211 0 : for (unsigned int i=1; i!=this->_m; i++)
1212 : {
1213 0 : rowsum = 0.;
1214 0 : for (unsigned int j=0; j!=this->_n; j++)
1215 : {
1216 0 : rowsum += abs((*this)(i,j));
1217 : }
1218 0 : my_max = (my_max > rowsum? my_max : rowsum);
1219 : }
1220 0 : return my_max;
1221 : }
1222 :
1223 :
1224 :
1225 : template<typename T>
1226 : inline
1227 0 : T DenseMatrix<T>::transpose (const unsigned int i,
1228 : const unsigned int j) const
1229 : {
1230 : // Implement in terms of operator()
1231 0 : return (*this)(j,i);
1232 : }
1233 :
1234 :
1235 :
1236 :
1237 :
1238 : // template<typename T>
1239 : // inline
1240 : // void DenseMatrix<T>::condense(const unsigned int iv,
1241 : // const unsigned int jv,
1242 : // const T val,
1243 : // DenseVector<T> & rhs)
1244 : // {
1245 : // libmesh_assert_equal_to (this->_m, rhs.size());
1246 : // libmesh_assert_equal_to (iv, jv);
1247 :
1248 :
1249 : // // move the known value into the RHS
1250 : // // and zero the column
1251 : // for (auto i : make_range(this->m()))
1252 : // {
1253 : // rhs(i) -= ((*this)(i,jv))*val;
1254 : // (*this)(i,jv) = 0.;
1255 : // }
1256 :
1257 : // // zero the row
1258 : // for (auto j : make_range(this->n()))
1259 : // (*this)(iv,j) = 0.;
1260 :
1261 : // (*this)(iv,jv) = 1.;
1262 : // rhs(iv) = val;
1263 :
1264 : // }
1265 :
1266 :
1267 : } // namespace libMesh
1268 :
1269 : #ifdef LIBMESH_HAVE_METAPHYSICL
1270 : namespace MetaPhysicL
1271 : {
1272 : template <typename T>
1273 : struct RawType<libMesh::DenseMatrix<T>>
1274 : {
1275 : typedef libMesh::DenseMatrix<typename RawType<T>::value_type> value_type;
1276 :
1277 : static value_type value (const libMesh::DenseMatrix<T> & in)
1278 : {
1279 : const auto m = in.m(), n = in.n();
1280 : value_type ret(m, n);
1281 : for (unsigned int i = 0; i < m; ++i)
1282 : for (unsigned int j = 0; j < n; ++j)
1283 : ret(i,j) = raw_value(in(i,j));
1284 :
1285 : return ret;
1286 : }
1287 : };
1288 : }
1289 : #endif
1290 :
1291 :
1292 : #endif // LIBMESH_DENSE_MATRIX_H
|
