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