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 : #ifndef LIBMESH_DENSE_MATRIX_IMPL_H
20 : #define LIBMESH_DENSE_MATRIX_IMPL_H
21 :
22 : // C++ Includes
23 : #include <cstdlib> // *must* precede <cmath> for proper std:abs() on PGI, Sun Studio CC
24 : #include <cmath> // for sqrt
25 :
26 : // Local Includes
27 : #include "libmesh/dense_matrix.h"
28 : #include "libmesh/dense_vector.h"
29 : #include "libmesh/int_range.h"
30 : #include "libmesh/libmesh.h"
31 :
32 : #ifdef LIBMESH_HAVE_METAPHYSICL
33 : #include "metaphysicl/dualnumber_decl.h"
34 : #endif
35 :
36 : namespace libMesh
37 : {
38 :
39 :
40 :
41 : // ------------------------------------------------------------
42 : // Dense Matrix member functions
43 :
44 : template<typename T>
45 0 : void DenseMatrix<T>::left_multiply (const DenseMatrixBase<T> & M2)
46 : {
47 0 : if (this->use_blas_lapack)
48 0 : this->_multiply_blas(M2, LEFT_MULTIPLY);
49 : else
50 : {
51 : // (*this) <- M2 * (*this)
52 : // Where:
53 : // (*this) = (m x n),
54 : // M2 = (m x p),
55 : // M3 = (p x n)
56 :
57 : // M3 is a copy of *this before it gets resize()d
58 0 : DenseMatrix<T> M3(*this);
59 :
60 : // Resize *this so that the result can fit
61 0 : this->resize (M2.m(), M3.n());
62 :
63 : // Call the multiply function in the base class
64 0 : this->multiply(*this, M2, M3);
65 0 : }
66 0 : }
67 :
68 :
69 :
70 : template<typename T>
71 : template<typename T2>
72 : void DenseMatrix<T>::left_multiply (const DenseMatrixBase<T2> & M2)
73 : {
74 : // (*this) <- M2 * (*this)
75 : // Where:
76 : // (*this) = (m x n),
77 : // M2 = (m x p),
78 : // M3 = (p x n)
79 :
80 : // M3 is a copy of *this before it gets resize()d
81 : DenseMatrix<T> M3(*this);
82 :
83 : // Resize *this so that the result can fit
84 : this->resize (M2.m(), M3.n());
85 :
86 : // Call the multiply function in the base class
87 : this->multiply(*this, M2, M3);
88 : }
89 :
90 :
91 :
92 : template<typename T>
93 1462935 : void DenseMatrix<T>::left_multiply_transpose(const DenseMatrix<T> & A)
94 : {
95 1462935 : if (this->use_blas_lapack)
96 1372844 : this->_multiply_blas(A, LEFT_MULTIPLY_TRANSPOSE);
97 : else
98 90091 : this->_left_multiply_transpose(A);
99 1462935 : }
100 :
101 :
102 : template<typename T>
103 : template<typename T2>
104 : void DenseMatrix<T>::left_multiply_transpose(const DenseMatrix<T2> & A)
105 : {
106 : this->_left_multiply_transpose(A);
107 : }
108 :
109 :
110 : template<typename T>
111 : template<typename T2>
112 90091 : void DenseMatrix<T>::_left_multiply_transpose(const DenseMatrix<T2> & A)
113 : {
114 : //Check to see if we are doing (A^T)*A
115 90091 : if (this == &A)
116 : {
117 : //libmesh_here();
118 0 : DenseMatrix<T> B(*this);
119 :
120 : // Simple but inefficient way
121 : // return this->left_multiply_transpose(B);
122 :
123 : // More efficient, but more code way
124 : // If A is mxn, the result will be a square matrix of Size n x n.
125 0 : const unsigned int n_rows = A.m();
126 0 : const unsigned int n_cols = A.n();
127 :
128 : // resize() *this and also zero out all entries.
129 0 : this->resize(n_cols,n_cols);
130 :
131 : // Compute the lower-triangular part
132 0 : for (unsigned int i=0; i<n_cols; ++i)
133 0 : for (unsigned int j=0; j<=i; ++j)
134 0 : for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
135 0 : (*this)(i,j) += B(k,i)*B(k,j);
136 :
137 : // Copy lower-triangular part into upper-triangular part
138 0 : for (unsigned int i=0; i<n_cols; ++i)
139 0 : for (unsigned int j=i+1; j<n_cols; ++j)
140 0 : (*this)(i,j) = (*this)(j,i);
141 0 : }
142 :
143 : else
144 : {
145 90091 : DenseMatrix<T> B(*this);
146 :
147 90091 : this->resize (A.n(), B.n());
148 :
149 0 : libmesh_assert_equal_to (A.m(), B.m());
150 0 : libmesh_assert_equal_to (this->m(), A.n());
151 0 : libmesh_assert_equal_to (this->n(), B.n());
152 :
153 0 : const unsigned int m_s = A.n();
154 0 : const unsigned int p_s = A.m();
155 0 : const unsigned int n_s = this->n();
156 :
157 : // Do it this way because there is a
158 : // decent chance (at least for constraint matrices)
159 : // that A.transpose(i,k) = 0.
160 1235752 : for (unsigned int i=0; i<m_s; i++)
161 25054174 : for (unsigned int k=0; k<p_s; k++)
162 0 : if (A.transpose(i,k) != 0.)
163 23205925 : for (unsigned int j=0; j<n_s; j++)
164 0 : (*this)(i,j) += A.transpose(i,k)*B(k,j);
165 90091 : }
166 90091 : }
167 :
168 :
169 :
170 : template<typename T>
171 3121863 : void DenseMatrix<T>::right_multiply (const DenseMatrixBase<T> & M3)
172 : {
173 3121863 : if (this->use_blas_lapack)
174 2847448 : this->_multiply_blas(M3, RIGHT_MULTIPLY);
175 : else
176 : {
177 : // (*this) <- (*this) * M3
178 : // Where:
179 : // (*this) = (m x n),
180 : // M2 = (m x p),
181 : // M3 = (p x n)
182 :
183 : // M2 is a copy of *this before it gets resize()d
184 274415 : DenseMatrix<T> M2(*this);
185 :
186 : // Resize *this so that the result can fit
187 274415 : this->resize (M2.m(), M3.n());
188 :
189 274415 : this->multiply(*this, M2, M3);
190 274415 : }
191 3121863 : }
192 :
193 :
194 :
195 : template<typename T>
196 : template<typename T2>
197 : void DenseMatrix<T>::right_multiply (const DenseMatrixBase<T2> & M3)
198 : {
199 : // (*this) <- M3 * (*this)
200 : // Where:
201 : // (*this) = (m x n),
202 : // M2 = (m x p),
203 : // M3 = (p x n)
204 :
205 : // M2 is a copy of *this before it gets resize()d
206 : DenseMatrix<T> M2(*this);
207 :
208 : // Resize *this so that the result can fit
209 : this->resize (M2.m(), M3.n());
210 :
211 : this->multiply(*this, M2, M3);
212 : }
213 :
214 :
215 :
216 :
217 : template<typename T>
218 1658880 : void DenseMatrix<T>::right_multiply_transpose (const DenseMatrix<T> & B)
219 : {
220 1658880 : if (this->use_blas_lapack)
221 1474560 : this->_multiply_blas(B, RIGHT_MULTIPLY_TRANSPOSE);
222 : else
223 184320 : this->_right_multiply_transpose(B);
224 1658880 : }
225 :
226 :
227 :
228 : template<typename T>
229 : template<typename T2>
230 : void DenseMatrix<T>::right_multiply_transpose (const DenseMatrix<T2> & B)
231 : {
232 : this->_right_multiply_transpose(B);
233 : }
234 :
235 :
236 :
237 : template<typename T>
238 : template<typename T2>
239 184320 : void DenseMatrix<T>::_right_multiply_transpose (const DenseMatrix<T2> & B)
240 : {
241 : //Check to see if we are doing B*(B^T)
242 184320 : if (this == &B)
243 : {
244 : //libmesh_here();
245 0 : DenseMatrix<T> A(*this);
246 :
247 : // Simple but inefficient way
248 : // return this->right_multiply_transpose(A);
249 :
250 : // More efficient, more code
251 : // If B is mxn, the result will be a square matrix of Size m x m.
252 0 : const unsigned int n_rows = B.m();
253 0 : const unsigned int n_cols = B.n();
254 :
255 : // resize() *this and also zero out all entries.
256 0 : this->resize(n_rows,n_rows);
257 :
258 : // Compute the lower-triangular part
259 0 : for (unsigned int i=0; i<n_rows; ++i)
260 0 : for (unsigned int j=0; j<=i; ++j)
261 0 : for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
262 0 : (*this)(i,j) += A(i,k)*A(j,k);
263 :
264 : // Copy lower-triangular part into upper-triangular part
265 0 : for (unsigned int i=0; i<n_rows; ++i)
266 0 : for (unsigned int j=i+1; j<n_rows; ++j)
267 0 : (*this)(i,j) = (*this)(j,i);
268 0 : }
269 :
270 : else
271 : {
272 184320 : DenseMatrix<T> A(*this);
273 :
274 184320 : this->resize (A.m(), B.m());
275 :
276 0 : libmesh_assert_equal_to (A.n(), B.n());
277 0 : libmesh_assert_equal_to (this->m(), A.m());
278 0 : libmesh_assert_equal_to (this->n(), B.m());
279 :
280 0 : const unsigned int m_s = A.m();
281 0 : const unsigned int p_s = A.n();
282 0 : const unsigned int n_s = this->n();
283 :
284 : // Do it this way because there is a
285 : // decent chance (at least for constraint matrices)
286 : // that B.transpose(k,j) = 0.
287 737280 : for (unsigned int j=0; j<n_s; j++)
288 3870720 : for (unsigned int k=0; k<p_s; k++)
289 3317760 : if (B.transpose(k,j) != 0.)
290 6635520 : for (unsigned int i=0; i<m_s; i++)
291 4976640 : (*this)(i,j) += A(i,k)*B.transpose(k,j);
292 184320 : }
293 184320 : }
294 :
295 :
296 :
297 :
298 : template<typename T>
299 27557811 : void DenseMatrix<T>::vector_mult (DenseVector<T> & dest,
300 : const DenseVector<T> & arg) const
301 : {
302 : // Make sure the input sizes are compatible
303 2290782 : libmesh_assert_equal_to (this->n(), arg.size());
304 :
305 : // Resize and clear dest.
306 : // Note: DenseVector::resize() also zeros the vector.
307 27557811 : dest.resize(this->m());
308 :
309 : // Short-circuit if the matrix is empty
310 27557811 : if(this->m() == 0 || this->n() == 0)
311 10000 : return;
312 :
313 27447811 : if (this->use_blas_lapack)
314 25632945 : this->_matvec_blas(1., 0., dest, arg);
315 : else
316 : {
317 0 : const unsigned int n_rows = this->m();
318 0 : const unsigned int n_cols = this->n();
319 :
320 6282190 : for (unsigned int i=0; i<n_rows; i++)
321 19415104 : for (unsigned int j=0; j<n_cols; j++)
322 10648800 : dest(i) += (*this)(i,j)*arg(j);
323 : }
324 : }
325 :
326 :
327 :
328 : template<typename T>
329 : template<typename T2>
330 0 : void DenseMatrix<T>::vector_mult (DenseVector<typename CompareTypes<T,T2>::supertype> & dest,
331 : const DenseVector<T2> & arg) const
332 : {
333 : // Make sure the input sizes are compatible
334 : libmesh_assert_equal_to (this->n(), arg.size());
335 :
336 : // Resize and clear dest.
337 : // Note: DenseVector::resize() also zeros the vector.
338 0 : dest.resize(this->m());
339 :
340 : // Short-circuit if the matrix is empty
341 0 : if (this->m() == 0 || this->n() == 0)
342 : return;
343 :
344 : const unsigned int n_rows = this->m();
345 : const unsigned int n_cols = this->n();
346 :
347 0 : for (unsigned int i=0; i<n_rows; i++)
348 0 : for (unsigned int j=0; j<n_cols; j++)
349 : dest(i) += (*this)(i,j)*arg(j);
350 : }
351 :
352 :
353 :
354 : template<typename T>
355 2869406 : void DenseMatrix<T>::vector_mult_transpose (DenseVector<T> & dest,
356 : const DenseVector<T> & arg) const
357 : {
358 : // Make sure the input sizes are compatible
359 265395 : libmesh_assert_equal_to (this->m(), arg.size());
360 :
361 : // Resize and clear dest.
362 : // Note: DenseVector::resize() also zeros the vector.
363 2869406 : dest.resize(this->n());
364 :
365 : // Short-circuit if the matrix is empty
366 2869406 : if (this->m() == 0)
367 0 : return;
368 :
369 2869406 : if (this->use_blas_lapack)
370 : {
371 2729923 : this->_matvec_blas(1., 0., dest, arg, /*trans=*/true);
372 : }
373 : else
374 : {
375 0 : const unsigned int n_rows = this->m();
376 0 : const unsigned int n_cols = this->n();
377 :
378 : // WORKS
379 : // for (unsigned int j=0; j<n_cols; j++)
380 : // for (unsigned int i=0; i<n_rows; i++)
381 : // dest(j) += (*this)(i,j)*arg(i);
382 :
383 : // ALSO WORKS, (i,j) just swapped
384 2446076 : for (unsigned int i=0; i<n_cols; i++)
385 55990312 : for (unsigned int j=0; j<n_rows; j++)
386 0 : dest(i) += (*this)(j,i)*arg(j);
387 : }
388 : }
389 :
390 :
391 :
392 : template<typename T>
393 : template<typename T2>
394 0 : void DenseMatrix<T>::vector_mult_transpose (DenseVector<typename CompareTypes<T,T2>::supertype> & dest,
395 : const DenseVector<T2> & arg) const
396 : {
397 : // Make sure the input sizes are compatible
398 : libmesh_assert_equal_to (this->m(), arg.size());
399 :
400 : // Resize and clear dest.
401 : // Note: DenseVector::resize() also zeros the vector.
402 0 : dest.resize(this->n());
403 :
404 : // Short-circuit if the matrix is empty
405 0 : if (this->m() == 0)
406 : return;
407 :
408 : const unsigned int n_rows = this->m();
409 : const unsigned int n_cols = this->n();
410 :
411 : // WORKS
412 : // for (unsigned int j=0; j<n_cols; j++)
413 : // for (unsigned int i=0; i<n_rows; i++)
414 : // dest(j) += (*this)(i,j)*arg(i);
415 :
416 : // ALSO WORKS, (i,j) just swapped
417 0 : for (unsigned int i=0; i<n_cols; i++)
418 0 : for (unsigned int j=0; j<n_rows; j++)
419 : dest(i) += (*this)(j,i)*arg(j);
420 : }
421 :
422 :
423 :
424 : template<typename T>
425 0 : void DenseMatrix<T>::vector_mult_add (DenseVector<T> & dest,
426 : const T factor,
427 : const DenseVector<T> & arg) const
428 : {
429 : // Short-circuit if the matrix is empty
430 0 : if (this->m() == 0)
431 : {
432 0 : dest.resize(0);
433 0 : return;
434 : }
435 :
436 0 : if (this->use_blas_lapack)
437 0 : this->_matvec_blas(factor, 1., dest, arg);
438 : else
439 : {
440 0 : DenseVector<T> temp(arg.size());
441 0 : this->vector_mult(temp, arg);
442 0 : dest.add(factor, temp);
443 : }
444 : }
445 :
446 :
447 :
448 : template<typename T>
449 : template<typename T2, typename T3>
450 0 : void DenseMatrix<T>::vector_mult_add (DenseVector<typename CompareTypes<T, typename CompareTypes<T2,T3>::supertype>::supertype> & dest,
451 : const T2 factor,
452 : const DenseVector<T3> & arg) const
453 : {
454 : // Short-circuit if the matrix is empty
455 0 : if (this->m() == 0)
456 : {
457 0 : dest.resize(0);
458 0 : return;
459 : }
460 :
461 : DenseVector<typename CompareTypes<T,T3>::supertype>
462 0 : temp(arg.size());
463 0 : this->vector_mult(temp, arg);
464 0 : dest.add(factor, temp);
465 : }
466 :
467 :
468 :
469 : template<typename T>
470 1244422 : void DenseMatrix<T>::get_principal_submatrix (unsigned int sub_m,
471 : unsigned int sub_n,
472 : DenseMatrix<T> & dest) const
473 : {
474 104066 : libmesh_assert( (sub_m <= this->m()) && (sub_n <= this->n()) );
475 :
476 1140356 : dest.resize(sub_m, sub_n);
477 10361690 : for (unsigned int i=0; i<sub_m; i++)
478 102478626 : for (unsigned int j=0; j<sub_n; j++)
479 93361358 : dest(i,j) = (*this)(i,j);
480 1244422 : }
481 :
482 :
483 :
484 : template<typename T>
485 70 : void DenseMatrix<T>::outer_product (const DenseVector<T> & a,
486 : const DenseVector<T> & b)
487 : {
488 2 : const unsigned int m = a.size();
489 2 : const unsigned int n = b.size();
490 :
491 68 : this->resize(m, n);
492 210 : for (unsigned int i = 0; i < m; ++i)
493 560 : for (unsigned int j = 0; j < n; ++j)
494 444 : (*this)(i,j) = a(i) * libmesh_conj(b(j));
495 70 : }
496 :
497 :
498 :
499 : template<typename T>
500 1244422 : void DenseMatrix<T>::get_principal_submatrix (unsigned int sub_m, DenseMatrix<T> & dest) const
501 : {
502 1244422 : get_principal_submatrix(sub_m, sub_m, dest);
503 1244422 : }
504 :
505 :
506 :
507 : template<typename T>
508 35459 : void DenseMatrix<T>::get_transpose (DenseMatrix<T> & dest) const
509 : {
510 35459 : dest.resize(this->n(), this->m());
511 :
512 180745 : for (auto i : make_range(dest.m()))
513 825328 : for (auto j : make_range(dest.n()))
514 680042 : dest(i,j) = (*this)(j,i);
515 35459 : }
516 :
517 :
518 :
519 :
520 : template<typename T>
521 2925480 : void DenseMatrix<T>::lu_solve (const DenseVector<T> & b,
522 : DenseVector<T> & x)
523 : {
524 : // Check to be sure that the matrix is square before attempting
525 : // an LU-solve. In general, one can compute the LU factorization of
526 : // a non-square matrix, but:
527 : //
528 : // Overdetermined systems (m>n) have a solution only if enough of
529 : // the equations are linearly-dependent.
530 : //
531 : // Underdetermined systems (m<n) typically have infinitely many
532 : // solutions.
533 : //
534 : // We don't want to deal with either of these ambiguous cases here...
535 253040 : libmesh_assert_equal_to (this->m(), this->n());
536 :
537 2925480 : switch(this->_decomposition_type)
538 : {
539 807050 : case NONE:
540 : {
541 807050 : if (this->use_blas_lapack)
542 758121 : this->_lu_decompose_lapack();
543 : else
544 48929 : this->_lu_decompose ();
545 60518 : break;
546 : }
547 :
548 2116978 : case LU_BLAS_LAPACK:
549 : {
550 : // Already factored, just need to call back_substitute.
551 2116978 : if (this->use_blas_lapack)
552 192522 : break;
553 : }
554 : libmesh_fallthrough();
555 :
556 : case LU:
557 : {
558 : // Already factored, just need to call back_substitute.
559 1452 : if (!(this->use_blas_lapack))
560 0 : break;
561 : }
562 : libmesh_fallthrough();
563 :
564 : default:
565 0 : libmesh_error_msg("Error! This matrix already has a different decomposition...");
566 : }
567 :
568 2925480 : if (this->use_blas_lapack)
569 2875099 : this->_lu_back_substitute_lapack (b, x);
570 : else
571 50381 : this->_lu_back_substitute (b, x);
572 2925480 : }
573 :
574 :
575 :
576 :
577 :
578 :
579 : template<typename T>
580 50381 : void DenseMatrix<T>::_lu_back_substitute (const DenseVector<T> & b,
581 : DenseVector<T> & x ) const
582 : {
583 : const unsigned int
584 0 : n_cols = this->n();
585 :
586 0 : libmesh_assert_equal_to (this->m(), n_cols);
587 0 : libmesh_assert_equal_to (this->m(), b.size());
588 :
589 50381 : x.resize (n_cols);
590 :
591 : // A convenient reference to *this
592 0 : const DenseMatrix<T> & A = *this;
593 :
594 : // Temporary vector storage. We use this instead of
595 : // modifying the RHS.
596 0 : DenseVector<T> z = b;
597 :
598 : // Lower-triangular "top to bottom" solve step, taking into account pivots
599 460467 : for (unsigned int i=0; i<n_cols; ++i)
600 : {
601 : // Swap
602 410086 : if (_pivots[i] != static_cast<pivot_index_t>(i))
603 173216 : std::swap( z(i), z(_pivots[i]) );
604 :
605 410086 : x(i) = z(i);
606 :
607 2242650 : for (unsigned int j=0; j<i; ++j)
608 0 : x(i) -= A(i,j)*x(j);
609 :
610 0 : x(i) /= A(i,i);
611 : }
612 :
613 : // Upper-triangular "bottom to top" solve step
614 50381 : const unsigned int last_row = n_cols-1;
615 :
616 460467 : for (int i=last_row; i>=0; --i)
617 : {
618 2242650 : for (int j=i+1; j<static_cast<int>(n_cols); ++j)
619 1832564 : x(i) -= A(i,j)*x(j);
620 : }
621 50381 : }
622 :
623 :
624 :
625 :
626 :
627 :
628 :
629 :
630 : template<typename T>
631 48929 : void DenseMatrix<T>::_lu_decompose ()
632 : {
633 : using std::abs;
634 :
635 : // If this function was called, there better not be any
636 : // previous decomposition of the matrix.
637 0 : libmesh_assert_equal_to (this->_decomposition_type, NONE);
638 :
639 : // Get the matrix size and make sure it is square
640 : const unsigned int
641 0 : n_rows = this->m();
642 :
643 : // A convenient reference to *this
644 0 : DenseMatrix<T> & A = *this;
645 :
646 48929 : _pivots.resize(n_rows);
647 :
648 450303 : for (unsigned int i=0; i<n_rows; ++i)
649 : {
650 : // Find the pivot row by searching down the i'th column
651 401374 : _pivots[i] = i;
652 :
653 : // abs(complex) must return a Real!
654 0 : auto the_max = abs( A(i,i) );
655 2212158 : for (unsigned int j=i+1; j<n_rows; ++j)
656 : {
657 0 : auto candidate_max = abs( A(j,i) );
658 1810784 : if (the_max < candidate_max)
659 : {
660 0 : the_max = candidate_max;
661 241267 : _pivots[i] = j;
662 : }
663 : }
664 :
665 : // libMesh::out << "the_max=" << the_max << " found at row " << _pivots[i] << std::endl;
666 :
667 : // If the max was found in a different row, interchange rows.
668 : // Here we interchange the *entire* row, in Gaussian elimination
669 : // you would only interchange the subrows A(i,j) and A(p(i),j), for j>i
670 401374 : if (_pivots[i] != static_cast<pivot_index_t>(i))
671 : {
672 1838455 : for (unsigned int j=0; j<n_rows; ++j)
673 1666589 : std::swap( A(i,j), A(_pivots[i], j) );
674 : }
675 :
676 : // If the max abs entry found is zero, the matrix is singular
677 0 : libmesh_error_msg_if(A(i,i) == libMesh::zero, "Matrix A is singular!");
678 :
679 : // Scale upper triangle entries of row i by the diagonal entry
680 : // Note: don't scale the diagonal entry itself!
681 0 : const T diag_inv = 1. / A(i,i);
682 2212158 : for (unsigned int j=i+1; j<n_rows; ++j)
683 0 : A(i,j) *= diag_inv;
684 :
685 : // Update the remaining sub-matrix A[i+1:m][i+1:m]
686 : // by subtracting off (the diagonal-scaled)
687 : // upper-triangular part of row i, scaled by the
688 : // i'th column entry of each row. In terms of
689 : // row operations, this is:
690 : // for each r > i
691 : // SubRow(r) = SubRow(r) - A(r,i)*SubRow(i)
692 : //
693 : // If we were scaling the i'th column as well, like
694 : // in Gaussian elimination, this would 'zero' the
695 : // entry in the i'th column.
696 2212158 : for (unsigned int row=i+1; row<n_rows; ++row)
697 14771792 : for (unsigned int col=i+1; col<n_rows; ++col)
698 0 : A(row,col) -= A(row,i) * A(i,col);
699 :
700 : } // end i loop
701 :
702 : // Set the flag for LU decomposition
703 48929 : this->_decomposition_type = LU;
704 48929 : }
705 :
706 :
707 :
708 : template<typename T>
709 0 : void DenseMatrix<T>::svd (DenseVector<Real> & sigma)
710 : {
711 : // We use the LAPACK svd implementation
712 0 : _svd_lapack(sigma);
713 0 : }
714 :
715 :
716 : template<typename T>
717 1470 : void DenseMatrix<T>::svd (DenseVector<Real> & sigma,
718 : DenseMatrix<Number> & U,
719 : DenseMatrix<Number> & VT)
720 : {
721 : // We use the LAPACK svd implementation
722 1470 : _svd_lapack(sigma, U, VT);
723 1470 : }
724 :
725 :
726 :
727 : template<typename T>
728 5654 : void DenseMatrix<T>::svd_solve(const DenseVector<T> & rhs,
729 : DenseVector<T> & x,
730 : Real rcond) const
731 : {
732 5654 : _svd_solve_lapack(rhs, x, rcond);
733 5654 : }
734 :
735 :
736 :
737 : template<typename T>
738 0 : void DenseMatrix<T>::evd (DenseVector<T> & lambda_real,
739 : DenseVector<T> & lambda_imag)
740 : {
741 : // We use the LAPACK eigenvalue problem implementation
742 0 : _evd_lapack(lambda_real, lambda_imag);
743 0 : }
744 :
745 :
746 :
747 : template<typename T>
748 0 : void DenseMatrix<T>::evd_left(DenseVector<T> & lambda_real,
749 : DenseVector<T> & lambda_imag,
750 : DenseMatrix<T> & VL)
751 : {
752 : // We use the LAPACK eigenvalue problem implementation
753 0 : _evd_lapack(lambda_real, lambda_imag, &VL, nullptr);
754 0 : }
755 :
756 :
757 :
758 : template<typename T>
759 0 : void DenseMatrix<T>::evd_right(DenseVector<T> & lambda_real,
760 : DenseVector<T> & lambda_imag,
761 : DenseMatrix<T> & VR)
762 : {
763 : // We use the LAPACK eigenvalue problem implementation
764 0 : _evd_lapack(lambda_real, lambda_imag, nullptr, &VR);
765 0 : }
766 :
767 :
768 :
769 : template<typename T>
770 140 : void DenseMatrix<T>::evd_left_and_right(DenseVector<T> & lambda_real,
771 : DenseVector<T> & lambda_imag,
772 : DenseMatrix<T> & VL,
773 : DenseMatrix<T> & VR)
774 : {
775 : // We use the LAPACK eigenvalue problem implementation
776 140 : _evd_lapack(lambda_real, lambda_imag, &VL, &VR);
777 140 : }
778 :
779 :
780 :
781 : template<typename T>
782 0 : T DenseMatrix<T>::det ()
783 : {
784 0 : switch(this->_decomposition_type)
785 : {
786 0 : case NONE:
787 : {
788 : // First LU decompose the matrix.
789 : // Note that the lu_decompose routine will check to see if the
790 : // matrix is square so we don't worry about it.
791 0 : if (this->use_blas_lapack)
792 0 : this->_lu_decompose_lapack();
793 : else
794 0 : this->_lu_decompose ();
795 : }
796 : case LU:
797 : case LU_BLAS_LAPACK:
798 : {
799 : // Already decomposed, don't do anything
800 0 : break;
801 : }
802 0 : default:
803 0 : libmesh_error_msg("Error! Can't compute the determinant under the current decomposition.");
804 : }
805 :
806 : // A variable to keep track of the running product of diagonal terms.
807 0 : T determinant = 1.;
808 :
809 : // Loop over diagonal terms, computing the product. In practice,
810 : // be careful because this value could easily become too large to
811 : // fit in a double or float. To be safe, one should keep track of
812 : // the power (of 10) of the determinant in a separate variable
813 : // and maintain an order 1 value for the determinant itself.
814 0 : unsigned int n_interchanges = 0;
815 0 : for (auto i : make_range(this->m()))
816 : {
817 0 : if (this->_decomposition_type==LU)
818 0 : if (_pivots[i] != static_cast<pivot_index_t>(i))
819 0 : n_interchanges++;
820 :
821 : // Lapack pivots are 1-based!
822 0 : if (this->_decomposition_type==LU_BLAS_LAPACK)
823 0 : if (_pivots[i] != static_cast<pivot_index_t>(i+1))
824 0 : n_interchanges++;
825 :
826 0 : determinant *= (*this)(i,i);
827 : }
828 :
829 : // Compute sign of determinant, depends on number of row interchanges!
830 : // The sign should be (-1)^{n}, where n is the number of interchanges.
831 0 : Real sign = n_interchanges % 2 == 0 ? 1. : -1.;
832 :
833 0 : return sign*determinant;
834 : }
835 :
836 :
837 :
838 : // The cholesky solve function first decomposes the matrix
839 : // with cholesky_decompose and then uses the cholesky_back_substitute
840 : // routine to find the solution x.
841 : template <typename T>
842 : template <typename T2>
843 1030566 : void DenseMatrix<T>::cholesky_solve (const DenseVector<T2> & b,
844 : DenseVector<T2> & x)
845 : {
846 : // Check for a previous decomposition
847 1030566 : switch(this->_decomposition_type)
848 : {
849 771914 : case NONE:
850 : {
851 771914 : this->_cholesky_decompose ();
852 771914 : break;
853 : }
854 :
855 21459 : case CHOLESKY:
856 : {
857 : // Already factored, just need to call back_substitute.
858 21459 : break;
859 : }
860 :
861 0 : default:
862 0 : libmesh_error_msg("Error! This matrix already has a different decomposition...");
863 : }
864 :
865 : // Perform back substitution
866 1030566 : this->_cholesky_back_substitute (b, x);
867 1030566 : }
868 :
869 :
870 :
871 :
872 : // This algorithm is based on the Cholesky decomposition in
873 : // the Numerical Recipes in C book.
874 : template<typename T>
875 771914 : void DenseMatrix<T>::_cholesky_decompose ()
876 : {
877 : using std::sqrt;
878 :
879 : // If we called this function, there better not be any
880 : // previous decomposition of the matrix.
881 62303 : libmesh_assert_equal_to (this->_decomposition_type, NONE);
882 :
883 : // Shorthand notation for number of rows and columns.
884 : const unsigned int
885 124606 : n_rows = this->m(),
886 124606 : n_cols = this->n();
887 :
888 : // Just to be really sure...
889 62303 : libmesh_assert_equal_to (n_rows, n_cols);
890 :
891 : // A convenient reference to *this
892 62303 : DenseMatrix<T> & A = *this;
893 :
894 5675372 : for (unsigned int i=0; i<n_rows; ++i)
895 : {
896 49404091 : for (unsigned int j=i; j<n_cols; ++j)
897 : {
898 541779784 : for (unsigned int k=0; k<i; ++k)
899 579265893 : A(i,j) -= A(i,k) * A(j,k);
900 :
901 44500633 : if (i == j)
902 : {
903 : #ifndef LIBMESH_USE_COMPLEX_NUMBERS
904 4520189 : libmesh_error_msg_if(A(i,j) <= 0.0,
905 : "Error! Can only use Cholesky decomposition with symmetric positive definite matrices.");
906 : #endif
907 :
908 5299235 : A(i,i) = sqrt(A(i,j));
909 : }
910 : else
911 46046291 : A(j,i) = A(i,j) / A(i,i);
912 : }
913 : }
914 :
915 : // Set the flag for CHOLESKY decomposition
916 771914 : this->_decomposition_type = CHOLESKY;
917 771914 : }
918 :
919 :
920 :
921 : template <typename T>
922 : template <typename T2>
923 1030566 : void DenseMatrix<T>::_cholesky_back_substitute (const DenseVector<T2> & b,
924 : DenseVector<T2> & x) const
925 : {
926 : // Shorthand notation for number of rows and columns.
927 : const unsigned int
928 167524 : n_rows = this->m(),
929 167524 : n_cols = this->n();
930 :
931 : // Just to be really sure...
932 83762 : libmesh_assert_equal_to (n_rows, n_cols);
933 :
934 : // A convenient reference to *this
935 83762 : const DenseMatrix<T> & A = *this;
936 :
937 : // Now compute the solution to Ax =b using the factorization.
938 946804 : x.resize(n_rows);
939 :
940 : // Solve for Ly=b
941 10352256 : for (unsigned int i=0; i<n_cols; ++i)
942 : {
943 9321690 : T2 temp = b(i);
944 :
945 147702985 : for (unsigned int k=0; k<i; ++k)
946 135370805 : temp -= A(i,k)*x(k);
947 :
948 10053903 : x(i) = temp / A(i,i);
949 : }
950 :
951 : // Solve for L^T x = y
952 10352256 : for (unsigned int i=0; i<n_cols; ++i)
953 : {
954 9321690 : const unsigned int ib = (n_cols-1)-i;
955 :
956 147702985 : for (unsigned int k=(ib+1); k<n_cols; ++k)
957 146367858 : x(ib) -= A(k,ib) * x(k);
958 :
959 8962884 : x(ib) /= A(ib,ib);
960 : }
961 1030566 : }
962 :
963 :
964 :
965 :
966 :
967 :
968 :
969 :
970 : // This routine is commented out since it is not really a memory
971 : // efficient implementation. Also, you don't *need* the inverse
972 : // for anything, instead just use lu_solve to solve Ax=b.
973 : // template<typename T>
974 : // void DenseMatrix<T>::inverse ()
975 : // {
976 : // // First LU decompose the matrix
977 : // // Note that the lu_decompose routine will check to see if the
978 : // // matrix is square so we don't worry about it.
979 : // if (!this->_lu_decomposed)
980 : // this->_lu_decompose();
981 :
982 : // // A unit vector which will be used as a rhs
983 : // // to pick off a single value each time.
984 : // DenseVector<T> e;
985 : // e.resize(this->m());
986 :
987 : // // An empty vector which will be used to hold the solution
988 : // // to the back substitutions.
989 : // DenseVector<T> x;
990 : // x.resize(this->m());
991 :
992 : // // An empty dense matrix to store the resulting inverse
993 : // // temporarily until we can overwrite A.
994 : // DenseMatrix<T> inv;
995 : // inv.resize(this->m(), this->n());
996 :
997 : // // Resize the passed in matrix to hold the inverse
998 : // inv.resize(this->m(), this->n());
999 :
1000 : // for (unsigned int j=0; j<this->n(); ++j)
1001 : // {
1002 : // e.zero();
1003 : // e(j) = 1.;
1004 : // this->_lu_back_substitute(e, x, false);
1005 : // for (unsigned int i=0; i<this->n(); ++i)
1006 : // inv(i,j) = x(i);
1007 : // }
1008 :
1009 : // // Now overwrite all the entries
1010 : // *this = inv;
1011 : // }
1012 :
1013 :
1014 : } // namespace libMesh
1015 :
1016 : #endif // LIBMESH_DENSE_MATRIX_IMPL_H
|
