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