dense_matrix.C
Go to the documentation of this file.
1 // The libMesh Finite Element Library.
2 // Copyright (C) 2002-2014 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 
63 template<typename T>
65 {
66  if (this->use_blas_lapack)
67  this->_multiply_blas(A, LEFT_MULTIPLY_TRANSPOSE);
68  else
69  {
70  //Check to see if we are doing (A^T)*A
71  if (this == &A)
72  {
73  //libmesh_here();
74  DenseMatrix<T> B(*this);
75 
76  // Simple but inefficient way
77  // return this->left_multiply_transpose(B);
78 
79  // More efficient, but more code way
80  // If A is mxn, the result will be a square matrix of Size n x n.
81  const unsigned int n_rows = A.m();
82  const unsigned int n_cols = A.n();
83 
84  // resize() *this and also zero out all entries.
85  this->resize(n_cols,n_cols);
86 
87  // Compute the lower-triangular part
88  for (unsigned int i=0; i<n_cols; ++i)
89  for (unsigned int j=0; j<=i; ++j)
90  for (unsigned int k=0; k<n_rows; ++k) // inner products are over n_rows
91  (*this)(i,j) += B(k,i)*B(k,j);
92 
93  // Copy lower-triangular part into upper-triangular part
94  for (unsigned int i=0; i<n_cols; ++i)
95  for (unsigned int j=i+1; j<n_cols; ++j)
96  (*this)(i,j) = (*this)(j,i);
97  }
98 
99  else
100  {
101  DenseMatrix<T> B(*this);
102 
103  this->resize (A.n(), B.n());
104 
105  libmesh_assert_equal_to (A.m(), B.m());
106  libmesh_assert_equal_to (this->m(), A.n());
107  libmesh_assert_equal_to (this->n(), B.n());
108 
109  const unsigned int m_s = A.n();
110  const unsigned int p_s = A.m();
111  const unsigned int n_s = this->n();
112 
113  // Do it this way because there is a
114  // decent chance (at least for constraint matrices)
115  // that A.transpose(i,k) = 0.
116  for (unsigned int i=0; i<m_s; i++)
117  for (unsigned int k=0; k<p_s; k++)
118  if (A.transpose(i,k) != 0.)
119  for (unsigned int j=0; j<n_s; j++)
120  (*this)(i,j) += A.transpose(i,k)*B(k,j);
121  }
122  }
123 
124 }
125 
126 
127 
128 
129 
130 
131 template<typename T>
133 {
134  if (this->use_blas_lapack)
135  this->_multiply_blas(M3, RIGHT_MULTIPLY);
136  else
137  {
138  // (*this) <- M3 * (*this)
139  // Where:
140  // (*this) = (m x n),
141  // M2 = (m x p),
142  // M3 = (p x n)
143 
144  // M2 is a copy of *this before it gets resize()d
145  DenseMatrix<T> M2(*this);
146 
147  // Resize *this so that the result can fit
148  this->resize (M2.m(), M3.n());
149 
150  this->multiply(*this, M2, M3);
151  }
152 }
153 
154 
155 
156 
157 template<typename T>
159 {
160  if (this->use_blas_lapack)
161  this->_multiply_blas(B, RIGHT_MULTIPLY_TRANSPOSE);
162  else
163  {
164  //Check to see if we are doing B*(B^T)
165  if (this == &B)
166  {
167  //libmesh_here();
168  DenseMatrix<T> A(*this);
169 
170  // Simple but inefficient way
171  // return this->right_multiply_transpose(A);
172 
173  // More efficient, more code
174  // If B is mxn, the result will be a square matrix of Size m x m.
175  const unsigned int n_rows = B.m();
176  const unsigned int n_cols = B.n();
177 
178  // resize() *this and also zero out all entries.
179  this->resize(n_rows,n_rows);
180 
181  // Compute the lower-triangular part
182  for (unsigned int i=0; i<n_rows; ++i)
183  for (unsigned int j=0; j<=i; ++j)
184  for (unsigned int k=0; k<n_cols; ++k) // inner products are over n_cols
185  (*this)(i,j) += A(i,k)*A(j,k);
186 
187  // Copy lower-triangular part into upper-triangular part
188  for (unsigned int i=0; i<n_rows; ++i)
189  for (unsigned int j=i+1; j<n_rows; ++j)
190  (*this)(i,j) = (*this)(j,i);
191  }
192 
193  else
194  {
195  DenseMatrix<T> A(*this);
196 
197  this->resize (A.m(), B.m());
198 
199  libmesh_assert_equal_to (A.n(), B.n());
200  libmesh_assert_equal_to (this->m(), A.m());
201  libmesh_assert_equal_to (this->n(), B.m());
202 
203  const unsigned int m_s = A.m();
204  const unsigned int p_s = A.n();
205  const unsigned int n_s = this->n();
206 
207  // Do it this way because there is a
208  // decent chance (at least for constraint matrices)
209  // that B.transpose(k,j) = 0.
210  for (unsigned int j=0; j<n_s; j++)
211  for (unsigned int k=0; k<p_s; k++)
212  if (B.transpose(k,j) != 0.)
213  for (unsigned int i=0; i<m_s; i++)
214  (*this)(i,j) += A(i,k)*B.transpose(k,j);
215  }
216  }
217 }
218 
219 
220 
221 template<typename T>
223  const DenseVector<T>& arg) const
224 {
225  // Make sure the input sizes are compatible
226  libmesh_assert_equal_to (this->n(), arg.size());
227 
228  // Resize and clear dest.
229  // Note: DenseVector::resize() also zeros the vector.
230  dest.resize(this->m());
231 
232  // Short-circuit if the matrix is empty
233  if(this->m() == 0 || this->n() == 0)
234  return;
235 
236  if (this->use_blas_lapack)
237  this->_matvec_blas(1., 0., dest, arg);
238  else
239  {
240  const unsigned int n_rows = this->m();
241  const unsigned int n_cols = this->n();
242 
243  for(unsigned int i=0; i<n_rows; i++)
244  for(unsigned int j=0; j<n_cols; j++)
245  dest(i) += (*this)(i,j)*arg(j);
246  }
247 }
248 
249 
250 
251 
252 template<typename T>
254  const DenseVector<T>& arg) const
255 {
256  // Make sure the input sizes are compatible
257  libmesh_assert_equal_to (this->m(), arg.size());
258 
259  // Resize and clear dest.
260  // Note: DenseVector::resize() also zeros the vector.
261  dest.resize(this->n());
262 
263  // Short-circuit if the matrix is empty
264  if(this->m() == 0)
265  return;
266 
267  if (this->use_blas_lapack)
268  {
269  this->_matvec_blas(1., 0., dest, arg, /*trans=*/true);
270  }
271  else
272  {
273  const unsigned int n_rows = this->m();
274  const unsigned int n_cols = this->n();
275 
276  // WORKS
277  // for(unsigned int j=0; j<n_cols; j++)
278  // for(unsigned int i=0; i<n_rows; i++)
279  // dest(j) += (*this)(i,j)*arg(i);
280 
281  // ALSO WORKS, (i,j) just swapped
282  for(unsigned int i=0; i<n_cols; i++)
283  for(unsigned int j=0; j<n_rows; j++)
284  dest(i) += (*this)(j,i)*arg(j);
285  }
286 }
287 
288 
289 
290 template<typename T>
292  const T factor,
293  const DenseVector<T>& arg) const
294 {
295  // Short-circuit if the matrix is empty
296  if(this->m() == 0)
297  {
298  dest.resize(0);
299  return;
300  }
301 
302  if (this->use_blas_lapack)
303  this->_matvec_blas(factor, 1., dest, arg);
304  else
305  {
306  DenseVector<T> temp(arg.size());
307  this->vector_mult(temp, arg);
308  dest.add(factor, temp);
309  }
310 }
311 
312 
313 
314 template<typename T>
316  unsigned int sub_n,
317  DenseMatrix<T>& dest) const
318 {
319  libmesh_assert( (sub_m <= this->m()) && (sub_n <= this->n()) );
320 
321  dest.resize(sub_m, sub_n);
322  for(unsigned int i=0; i<sub_m; i++)
323  for(unsigned int j=0; j<sub_n; j++)
324  dest(i,j) = (*this)(i,j);
325 }
326 
327 
328 
329 template<typename T>
330 void DenseMatrix<T>::get_principal_submatrix (unsigned int sub_m, DenseMatrix<T>& dest) const
331 {
332  get_principal_submatrix(sub_m, sub_m, dest);
333 }
334 
335 
336 
337 template<typename T>
339 {
340  dest.resize(this->n(), this->m());
341 
342  for (unsigned int i=0; i<dest.m(); i++)
343  for (unsigned int j=0; j<dest.n(); j++)
344  dest(i,j) = (*this)(j,i);
345 }
346 
347 
348 
349 
350 template<typename T>
352  DenseVector<T>& x)
353 {
354  // Check to be sure that the matrix is square before attempting
355  // an LU-solve. In general, one can compute the LU factorization of
356  // a non-square matrix, but:
357  //
358  // Overdetermined systems (m>n) have a solution only if enough of
359  // the equations are linearly-dependent.
360  //
361  // Underdetermined systems (m<n) typically have infinitely many
362  // solutions.
363  //
364  // We don't want to deal with either of these ambiguous cases here...
365  libmesh_assert_equal_to (this->m(), this->n());
366 
367  switch(this->_decomposition_type)
368  {
369  case NONE:
370  {
371  if (this->use_blas_lapack)
372  this->_lu_decompose_lapack();
373  else
374  this->_lu_decompose ();
375  break;
376  }
377 
378  case LU_BLAS_LAPACK:
379  {
380  // Already factored, just need to call back_substitute.
381  if (this->use_blas_lapack)
382  break;
383  }
384 
385  case LU:
386  {
387  // Already factored, just need to call back_substitute.
388  if ( !(this->use_blas_lapack) )
389  break;
390  }
391 
392  default:
393  {
394  libMesh::err << "Error! This matrix already has a "
395  << "different decomposition..."
396  << std::endl;
397  libmesh_error();
398  }
399  }
400 
401  if (this->use_blas_lapack)
402  this->_lu_back_substitute_lapack (b, x);
403  else
404  this->_lu_back_substitute (b, x);
405 }
406 
407 
408 
409 
410 
411 
412 template<typename T>
414  DenseVector<T>& x ) const
415 {
416  const unsigned int
417  n_cols = this->n();
418 
419  libmesh_assert_equal_to (this->m(), n_cols);
420  libmesh_assert_equal_to (this->m(), b.size());
421 
422  x.resize (n_cols);
423 
424  // A convenient reference to *this
425  const DenseMatrix<T>& A = *this;
426 
427  // Temporary vector storage. We use this instead of
428  // modifying the RHS.
429  DenseVector<T> z = b;
430 
431  // Lower-triangular "top to bottom" solve step, taking into account pivots
432  for (unsigned int i=0; i<n_cols; ++i)
433  {
434  // Swap
435  if (_pivots[i] != static_cast<int>(i))
436  std::swap( z(i), z(_pivots[i]) );
437 
438  x(i) = z(i);
439 
440  for (unsigned int j=0; j<i; ++j)
441  x(i) -= A(i,j)*x(j);
442 
443  x(i) /= A(i,i);
444  }
445 
446  // Upper-triangular "bottom to top" solve step
447  const unsigned int last_row = n_cols-1;
448 
449  for (int i=last_row; i>=0; --i)
450  {
451  for (int j=i+1; j<static_cast<int>(n_cols); ++j)
452  x(i) -= A(i,j)*x(j);
453  }
454 }
455 
456 
457 
458 
459 
460 
461 
462 
463 template<typename T>
465 {
466  // If this function was called, there better not be any
467  // previous decomposition of the matrix.
468  libmesh_assert_equal_to (this->_decomposition_type, NONE);
469 
470  // Get the matrix size and make sure it is square
471  const unsigned int
472  n_rows = this->m();
473 
474  // A convenient reference to *this
475  DenseMatrix<T>& A = *this;
476 
477  _pivots.resize(n_rows);
478 
479  for (unsigned int i=0; i<n_rows; ++i)
480  {
481  // Find the pivot row by searching down the i'th column
482  _pivots[i] = i;
483 
484  // std::abs(complex) must return a Real!
485  Real the_max = std::abs( A(i,i) );
486  for (unsigned int j=i+1; j<n_rows; ++j)
487  {
488  Real candidate_max = std::abs( A(j,i) );
489  if (the_max < candidate_max)
490  {
491  the_max = candidate_max;
492  _pivots[i] = j;
493  }
494  }
495 
496  // libMesh::out << "the_max=" << the_max << " found at row " << _pivots[i] << std::endl;
497 
498  // If the max was found in a different row, interchange rows.
499  // Here we interchange the *entire* row, in Gaussian elimination
500  // you would only interchange the subrows A(i,j) and A(p(i),j), for j>i
501  if (_pivots[i] != static_cast<int>(i))
502  {
503  for (unsigned int j=0; j<n_rows; ++j)
504  std::swap( A(i,j), A(_pivots[i], j) );
505  }
506 
507 
508  // If the max abs entry found is zero, the matrix is singular
509  if (A(i,i) == libMesh::zero)
510  {
511  libMesh::out << "Matrix A is singular!" << std::endl;
512  libmesh_error();
513  }
514 
515  // Scale upper triangle entries of row i by the diagonal entry
516  // Note: don't scale the diagonal entry itself!
517  const T diag_inv = 1. / A(i,i);
518  for (unsigned int j=i+1; j<n_rows; ++j)
519  A(i,j) *= diag_inv;
520 
521  // Update the remaining sub-matrix A[i+1:m][i+1:m]
522  // by subtracting off (the diagonal-scaled)
523  // upper-triangular part of row i, scaled by the
524  // i'th column entry of each row. In terms of
525  // row operations, this is:
526  // for each r > i
527  // SubRow(r) = SubRow(r) - A(r,i)*SubRow(i)
528  //
529  // If we were scaling the i'th column as well, like
530  // in Gaussian elimination, this would 'zero' the
531  // entry in the i'th column.
532  for (unsigned int row=i+1; row<n_rows; ++row)
533  for (unsigned int col=i+1; col<n_rows; ++col)
534  A(row,col) -= A(row,i) * A(i,col);
535 
536  } // end i loop
537 
538  // Set the flag for LU decomposition
539  this->_decomposition_type = LU;
540 }
541 
542 
543 
544 template<typename T>
546 {
547  // We use the LAPACK svd implementation
548  _svd_lapack(sigma);
549 }
550 
551 
552 template<typename T>
554 {
555  // We use the LAPACK svd implementation
556  _svd_lapack(sigma, U, VT);
557 }
558 
559 
560 
561 template<typename T>
563 {
564  switch(this->_decomposition_type)
565  {
566  case NONE:
567  {
568  // First LU decompose the matrix.
569  // Note that the lu_decompose routine will check to see if the
570  // matrix is square so we don't worry about it.
571  if (this->use_blas_lapack)
572  this->_lu_decompose_lapack();
573  else
574  this->_lu_decompose ();
575  }
576  case LU:
577  case LU_BLAS_LAPACK:
578  {
579  // Already decomposed, don't do anything
580  break;
581  }
582  default:
583  {
584  libMesh::err << "Error! Can't compute the determinant under "
585  << "the current decomposition."
586  << std::endl;
587  libmesh_error();
588  }
589  }
590 
591  // A variable to keep track of the running product of diagonal terms.
592  T determinant = 1.;
593 
594  // Loop over diagonal terms, computing the product. In practice,
595  // be careful because this value could easily become too large to
596  // fit in a double or float. To be safe, one should keep track of
597  // the power (of 10) of the determinant in a separate variable
598  // and maintain an order 1 value for the determinant itself.
599  unsigned int n_interchanges = 0;
600  for (unsigned int i=0; i<this->m(); i++)
601  {
602  if (this->_decomposition_type==LU)
603  if (_pivots[i] != static_cast<int>(i))
604  n_interchanges++;
605 
606  // Lapack pivots are 1-based!
607  if (this->_decomposition_type==LU_BLAS_LAPACK)
608  if (_pivots[i] != static_cast<int>(i+1))
609  n_interchanges++;
610 
611  determinant *= (*this)(i,i);
612  }
613 
614  // Compute sign of determinant, depends on number of row interchanges!
615  // The sign should be (-1)^{n}, where n is the number of interchanges.
616  Real sign = n_interchanges % 2 == 0 ? 1. : -1.;
617 
618  return sign*determinant;
619 }
620 
621 
622 
623 // The cholesky solve function first decomposes the matrix
624 // with cholesky_decompose and then uses the cholesky_back_substitute
625 // routine to find the solution x.
626 template <typename T>
627 template <typename T2>
629  DenseVector<T2>& x)
630 {
631  // Check for a previous decomposition
632  switch(this->_decomposition_type)
633  {
634  case NONE:
635  {
636  this->_cholesky_decompose ();
637  break;
638  }
639 
640  case CHOLESKY:
641  {
642  // Already factored, just need to call back_substitute.
643  break;
644  }
645 
646  default:
647  {
648  libMesh::err << "Error! This matrix already has a "
649  << "different decomposition..."
650  << std::endl;
651  libmesh_error();
652  }
653  }
654 
655  // Perform back substitution
656  this->_cholesky_back_substitute (b, x);
657 }
658 
659 
660 
661 
662 // This algorithm is based on the Cholesky decomposition in
663 // the Numerical Recipes in C book.
664 template<typename T>
666 {
667  // If we called this function, there better not be any
668  // previous decomposition of the matrix.
669  libmesh_assert_equal_to (this->_decomposition_type, NONE);
670 
671  // Shorthand notation for number of rows and columns.
672  const unsigned int
673  n_rows = this->m(),
674  n_cols = this->n();
675 
676  // Just to be really sure...
677  libmesh_assert_equal_to (n_rows, n_cols);
678 
679  // A convenient reference to *this
680  DenseMatrix<T>& A = *this;
681 
682  for (unsigned int i=0; i<n_rows; ++i)
683  {
684  for (unsigned int j=i; j<n_cols; ++j)
685  {
686  for (unsigned int k=0; k<i; ++k)
687  A(i,j) -= A(i,k) * A(j,k);
688 
689  if (i == j)
690  {
691 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
692  if (A(i,j) <= 0.0)
693  {
694  libMesh::err << "Error! Can only use Cholesky decomposition "
695  << "with symmetric positive definite matrices."
696  << std::endl;
697  libmesh_error();
698  }
699 #endif
700 
701  A(i,i) = std::sqrt(A(i,j));
702  }
703  else
704  A(j,i) = A(i,j) / A(i,i);
705  }
706  }
707 
708  // Set the flag for CHOLESKY decomposition
709  this->_decomposition_type = CHOLESKY;
710 }
711 
712 
713 
714 template <typename T>
715 template <typename T2>
717  DenseVector<T2>& x) const
718 {
719  // Shorthand notation for number of rows and columns.
720  const unsigned int
721  n_rows = this->m(),
722  n_cols = this->n();
723 
724  // Just to be really sure...
725  libmesh_assert_equal_to (n_rows, n_cols);
726 
727  // A convenient reference to *this
728  const DenseMatrix<T>& A = *this;
729 
730  // Now compute the solution to Ax =b using the factorization.
731  x.resize(n_rows);
732 
733  // Solve for Ly=b
734  for (unsigned int i=0; i<n_cols; ++i)
735  {
736  T2 temp = b(i);
737 
738  for (unsigned int k=0; k<i; ++k)
739  temp -= A(i,k)*x(k);
740 
741  x(i) = temp / A(i,i);
742  }
743 
744  // Solve for L^T x = y
745  for (unsigned int i=0; i<n_cols; ++i)
746  {
747  const unsigned int ib = (n_cols-1)-i;
748 
749  for (unsigned int k=(ib+1); k<n_cols; ++k)
750  x(ib) -= A(k,ib) * x(k);
751 
752  x(ib) /= A(ib,ib);
753  }
754 }
755 
756 
757 
758 
759 
760 
761 
762 
763 // This routine is commented out since it is not really a memory
764 // efficient implementation. Also, you don't *need* the inverse
765 // for anything, instead just use lu_solve to solve Ax=b.
766 // template<typename T>
767 // void DenseMatrix<T>::inverse ()
768 // {
769 // // First LU decompose the matrix
770 // // Note that the lu_decompose routine will check to see if the
771 // // matrix is square so we don't worry about it.
772 // if (!this->_lu_decomposed)
773 // this->_lu_decompose();
774 
775 // // A unit vector which will be used as a rhs
776 // // to pick off a single value each time.
777 // DenseVector<T> e;
778 // e.resize(this->m());
779 
780 // // An empty vector which will be used to hold the solution
781 // // to the back substitutions.
782 // DenseVector<T> x;
783 // x.resize(this->m());
784 
785 // // An empty dense matrix to store the resulting inverse
786 // // temporarily until we can overwrite A.
787 // DenseMatrix<T> inv;
788 // inv.resize(this->m(), this->n());
789 
790 // // Resize the passed in matrix to hold the inverse
791 // inv.resize(this->m(), this->n());
792 
793 // for (unsigned int j=0; j<this->n(); ++j)
794 // {
795 // e.zero();
796 // e(j) = 1.;
797 // this->_lu_back_substitute(e, x, false);
798 // for (unsigned int i=0; i<this->n(); ++i)
799 // inv(i,j) = x(i);
800 // }
801 
802 // // Now overwrite all the entries
803 // *this = inv;
804 // }
805 
806 
807 //--------------------------------------------------------------
808 // Explicit instantiations
809 template class DenseMatrix<Real>;
814 
815 #ifdef LIBMESH_USE_COMPLEX_NUMBERS
816 template class DenseMatrix<Complex>;
819 #endif
820 
821 } // namespace libMesh

Site Created By: libMesh Developers
Last modified: February 07 2014 16:57:04 UTC

Hosted By:
SourceForge.net Logo