dense_matrix.h
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 
20 #ifndef LIBMESH_DENSE_MATRIX_H
21 #define LIBMESH_DENSE_MATRIX_H
22 
23 // Local Includes
24 #include "libmesh/libmesh_common.h"
26 
27 // C++ includes
28 #include <vector>
29 #include <algorithm>
30 #include <cstring> // std::memset
31 
32 namespace libMesh
33 {
34 
35 // Forward Declarations
36 template <typename T> class DenseVector;
37 
38 
39 
48 // ------------------------------------------------------------
49 // Dense Matrix class definition
50 template<typename T>
51 class DenseMatrix : public DenseMatrixBase<T>
52 {
53 public:
54 
58  DenseMatrix(const unsigned int new_m=0,
59  const unsigned int new_n=0);
60 
64  //DenseMatrix (const DenseMatrix<T>& other_matrix);
65 
69  virtual ~DenseMatrix() {}
70 
71 
75  virtual void zero();
76 
80  T operator() (const unsigned int i,
81  const unsigned int j) const;
82 
86  T & operator() (const unsigned int i,
87  const unsigned int j);
88 
92  virtual T el(const unsigned int i,
93  const unsigned int j) const { return (*this)(i,j); }
94 
98  virtual T & el(const unsigned int i,
99  const unsigned int j) { return (*this)(i,j); }
100 
104  virtual void left_multiply (const DenseMatrixBase<T>& M2);
105 
109  virtual void right_multiply (const DenseMatrixBase<T>& M3);
110 
115  void vector_mult(DenseVector<T>& dest, const DenseVector<T>& arg) const;
116 
121  void vector_mult_transpose(DenseVector<T>& dest, const DenseVector<T>& arg) const;
122 
127  void vector_mult_add (DenseVector<T>& dest,
128  const T factor,
129  const DenseVector<T>& arg) const;
130 
134  void get_principal_submatrix (unsigned int sub_m, unsigned int sub_n, DenseMatrix<T>& dest) const;
135 
139  void get_principal_submatrix (unsigned int sub_m, DenseMatrix<T>& dest) const;
140 
144  DenseMatrix<T>& operator = (const DenseMatrix<T>& other_matrix);
145 
153  template <typename T2>
154  DenseMatrix<T>& operator = (const DenseMatrix<T2>& other_matrix);
155 
159  void swap(DenseMatrix<T>& other_matrix);
160 
165  void resize(const unsigned int new_m,
166  const unsigned int new_n);
167 
171  void scale (const T factor);
172 
173 
177  void scale_column (const unsigned int col, const T factor);
178 
182  DenseMatrix<T>& operator *= (const T factor);
183 
187  template<typename T2, typename T3>
188  typename boostcopy::enable_if_c<
189  ScalarTraits<T2>::value, void >::type add (const T2 factor,
190  const DenseMatrix<T3>& mat);
191 
195  bool operator== (const DenseMatrix<T> &mat) const;
196 
200  bool operator!= (const DenseMatrix<T> &mat) const;
201 
206 
211 
217  Real min () const;
218 
224  Real max () const;
225 
236  Real l1_norm () const;
237 
249  Real linfty_norm () const;
250 
255 
256 
261 
265  T transpose (const unsigned int i,
266  const unsigned int j) const;
267 
271  void get_transpose(DenseMatrix<T>& dest) const;
272 
278  std::vector<T>& get_values() { return _val; }
279 
283  const std::vector<T>& get_values() const { return _val; }
284 
291  void condense(const unsigned int i,
292  const unsigned int j,
293  const T val,
294  DenseVector<T>& rhs)
295  { DenseMatrixBase<T>::condense (i, j, val, rhs); }
296 
302  void lu_solve (const DenseVector<T>& b,
303  DenseVector<T>& x);
304 
305 
306 
316  template <typename T2>
317  void cholesky_solve(const DenseVector<T2>& b,
318  DenseVector<T2>& x);
319 
320 
329  void svd(DenseVector<T>& sigma);
330 
331 
343  void svd(DenseVector<T>& sigma, DenseMatrix<T>& U, DenseMatrix<T>& VT);
344 
345 
351  T det();
352 
361  // void inverse();
362 
369 
370 private:
371 
375  std::vector<T> _val;
376 
382  void _lu_decompose ();
383 
389  void _lu_back_substitute (const DenseVector<T>& b,
390  DenseVector<T>& x) const;
391 
398  void _cholesky_decompose();
399 
406  template <typename T2>
408  DenseVector<T2>& x) const;
409 
418 
424 
434  };
435 
443  void _multiply_blas(const DenseMatrixBase<T>& other,
444  _BLAS_Multiply_Flag flag);
445 
455  void _lu_decompose_lapack();
456 
462  void _svd_lapack(DenseVector<T>& sigma);
463 
469  void _svd_lapack(DenseVector<T>& sigma,
470  DenseMatrix<T>& U,
471  DenseMatrix<T>& VT);
472 
477  void _svd_helper (char JOBU,
478  char JOBVT,
479  std::vector<T>& sigma_val,
480  std::vector<T>& U_val,
481  std::vector<T>& VT_val);
482 
488  std::vector<int> _pivots;
489 
500  DenseVector<T>& x);
501 
514  void _matvec_blas(T alpha, T beta,
515  DenseVector<T>& dest,
516  const DenseVector<T>& arg,
517  bool trans=false) const;
518 };
519 
520 
521 
522 
523 
524 // ------------------------------------------------------------
528 namespace DenseMatrices
529 {
530 
536 
547 
548 }
549 
550 
551 
552 using namespace DenseMatrices;
553 
554 
555 
556 
557 
558 // ------------------------------------------------------------
559 // Dense Matrix member functions
560 template<typename T>
561 inline
562 DenseMatrix<T>::DenseMatrix(const unsigned int new_m,
563  const unsigned int new_n)
564  : DenseMatrixBase<T>(new_m,new_n),
565 #if defined(LIBMESH_HAVE_PETSC) && defined(LIBMESH_USE_REAL_NUMBERS) && defined(LIBMESH_DEFAULT_DOUBLE_PRECISION)
566  use_blas_lapack(true),
567 #else
568  use_blas_lapack(false),
569 #endif
570  _val(),
571  _decomposition_type(NONE),
572  _pivots()
573 {
574  this->resize(new_m,new_n);
575 }
576 
577 
578 
579 // FIXME[JWP]: This copy ctor has not been maintained along with
580 // the rest of the class...
581 // Can we just use the compiler-generated copy ctor here?
582 // template<typename T>
583 // inline
584 // DenseMatrix<T>::DenseMatrix (const DenseMatrix<T>& other_matrix)
585 // : DenseMatrixBase<T>(other_matrix._m, other_matrix._n)
586 // {
587 // _val = other_matrix._val;
588 // }
589 
590 
591 
592 template<typename T>
593 inline
595 {
596  std::swap(this->_m, other_matrix._m);
597  std::swap(this->_n, other_matrix._n);
598  _val.swap(other_matrix._val);
599  DecompositionType _temp = _decomposition_type;
600  _decomposition_type = other_matrix._decomposition_type;
601  other_matrix._decomposition_type = _temp;
602 }
603 
604 
605 template <typename T>
606 template <typename T2>
607 inline
610 {
611  unsigned int mat_m = mat.m(), mat_n = mat.n();
612  this->resize(mat_m, mat_n);
613  for (unsigned int i=0; i<mat_m; i++)
614  for (unsigned int j=0; j<mat_n; j++)
615  (*this)(i,j) = mat(i,j);
616 
617  return *this;
618 }
619 
620 
621 
622 template<typename T>
623 inline
624 void DenseMatrix<T>::resize(const unsigned int new_m,
625  const unsigned int new_n)
626 {
627  _val.resize(new_m*new_n);
628 
629  this->_m = new_m;
630  this->_n = new_n;
631 
632  _decomposition_type = NONE;
633  this->zero();
634 }
635 
636 
637 
638 template<typename T>
639 inline
641 {
642  _decomposition_type = NONE;
643 
644  // Just doing this ifdef to be completely safe
645 #ifndef LIBMESH_USE_COMPLEX_NUMBERS
646  if(_val.size())
647  std::memset(&_val[0], 0, sizeof(T) * _val.size());
648 #else
649  std::fill (_val.begin(), _val.end(), 0.);
650 #endif
651 }
652 
653 
654 
655 template<typename T>
656 inline
658 {
659  this->_m = other_matrix._m;
660  this->_n = other_matrix._n;
661 
662  _val = other_matrix._val;
663  _decomposition_type = other_matrix._decomposition_type;
664 
665  return *this;
666 }
667 
668 
669 
670 template<typename T>
671 inline
672 T DenseMatrix<T>::operator () (const unsigned int i,
673  const unsigned int j) const
674 {
675  libmesh_assert_less (i*j, _val.size());
676  libmesh_assert_less (i, this->_m);
677  libmesh_assert_less (j, this->_n);
678 
679 
680  // return _val[(i) + (this->_m)*(j)]; // col-major
681  return _val[(i)*(this->_n) + (j)]; // row-major
682 }
683 
684 
685 
686 template<typename T>
687 inline
688 T & DenseMatrix<T>::operator () (const unsigned int i,
689  const unsigned int j)
690 {
691  libmesh_assert_less (i*j, _val.size());
692  libmesh_assert_less (i, this->_m);
693  libmesh_assert_less (j, this->_n);
694 
695  //return _val[(i) + (this->_m)*(j)]; // col-major
696  return _val[(i)*(this->_n) + (j)]; // row-major
697 }
698 
699 
700 
701 
702 
703 template<typename T>
704 inline
705 void DenseMatrix<T>::scale (const T factor)
706 {
707  for (unsigned int i=0; i<_val.size(); i++)
708  _val[i] *= factor;
709 }
710 
711 
712 template<typename T>
713 inline
714 void DenseMatrix<T>::scale_column (const unsigned int col, const T factor)
715 {
716  for (unsigned int i=0; i<this->m(); i++)
717  (*this)(i, col) *= factor;
718 }
719 
720 
721 
722 template<typename T>
723 inline
725 {
726  this->scale(factor);
727  return *this;
728 }
729 
730 
731 
732 template<typename T>
733 template<typename T2, typename T3>
734 inline
735 typename boostcopy::enable_if_c<
736 ScalarTraits<T2>::value, void >::type
737 DenseMatrix<T>::add (const T2 factor,
738  const DenseMatrix<T3>& mat)
739 {
740  libmesh_assert_equal_to (this->m(), mat.m());
741  libmesh_assert_equal_to (this->n(), mat.n());
742 
743  for (unsigned int i=0; i<this->m(); i++)
744  for (unsigned int j=0; j<this->n(); j++)
745  (*this)(i,j) += factor * mat(i,j);
746 }
747 
748 
749 
750 template<typename T>
751 inline
753 {
754  for (unsigned int i=0; i<_val.size(); i++)
755  if (_val[i] != mat._val[i])
756  return false;
757 
758  return true;
759 }
760 
761 
762 
763 template<typename T>
764 inline
766 {
767  for (unsigned int i=0; i<_val.size(); i++)
768  if (_val[i] != mat._val[i])
769  return true;
770 
771  return false;
772 }
773 
774 
775 
776 template<typename T>
777 inline
779 {
780  for (unsigned int i=0; i<_val.size(); i++)
781  _val[i] += mat._val[i];
782 
783  return *this;
784 }
785 
786 
787 
788 template<typename T>
789 inline
791 {
792  for (unsigned int i=0; i<_val.size(); i++)
793  _val[i] -= mat._val[i];
794 
795  return *this;
796 }
797 
798 
799 
800 template<typename T>
801 inline
803 {
804  libmesh_assert (this->_m);
805  libmesh_assert (this->_n);
806  Real my_min = libmesh_real((*this)(0,0));
807 
808  for (unsigned int i=0; i!=this->_m; i++)
809  {
810  for (unsigned int j=0; j!=this->_n; j++)
811  {
812  Real current = libmesh_real((*this)(i,j));
813  my_min = (my_min < current? my_min : current);
814  }
815  }
816  return my_min;
817 }
818 
819 
820 
821 template<typename T>
822 inline
824 {
825  libmesh_assert (this->_m);
826  libmesh_assert (this->_n);
827  Real my_max = libmesh_real((*this)(0,0));
828 
829  for (unsigned int i=0; i!=this->_m; i++)
830  {
831  for (unsigned int j=0; j!=this->_n; j++)
832  {
833  Real current = libmesh_real((*this)(i,j));
834  my_max = (my_max > current? my_max : current);
835  }
836  }
837  return my_max;
838 }
839 
840 
841 
842 template<typename T>
843 inline
845 {
846  libmesh_assert (this->_m);
847  libmesh_assert (this->_n);
848 
849  Real columnsum = 0.;
850  for (unsigned int i=0; i!=this->_m; i++)
851  {
852  columnsum += std::abs((*this)(i,0));
853  }
854  Real my_max = columnsum;
855  for (unsigned int j=1; j!=this->_n; j++)
856  {
857  columnsum = 0.;
858  for (unsigned int i=0; i!=this->_m; i++)
859  {
860  columnsum += std::abs((*this)(i,j));
861  }
862  my_max = (my_max > columnsum? my_max : columnsum);
863  }
864  return my_max;
865 }
866 
867 
868 
869 template<typename T>
870 inline
872 {
873  libmesh_assert (this->_m);
874  libmesh_assert (this->_n);
875 
876  Real rowsum = 0.;
877  for (unsigned int j=0; j!=this->_n; j++)
878  {
879  rowsum += std::abs((*this)(0,j));
880  }
881  Real my_max = rowsum;
882  for (unsigned int i=1; i!=this->_m; i++)
883  {
884  rowsum = 0.;
885  for (unsigned int j=0; j!=this->_n; j++)
886  {
887  rowsum += std::abs((*this)(i,j));
888  }
889  my_max = (my_max > rowsum? my_max : rowsum);
890  }
891  return my_max;
892 }
893 
894 
895 
896 template<typename T>
897 inline
898 T DenseMatrix<T>::transpose (const unsigned int i,
899  const unsigned int j) const
900 {
901  // Implement in terms of operator()
902  return (*this)(j,i);
903 }
904 
905 
906 
907 
908 
909 // template<typename T>
910 // inline
911 // void DenseMatrix<T>::condense(const unsigned int iv,
912 // const unsigned int jv,
913 // const T val,
914 // DenseVector<T>& rhs)
915 // {
916 // libmesh_assert_equal_to (this->_m, rhs.size());
917 // libmesh_assert_equal_to (iv, jv);
918 
919 
920 // // move the known value into the RHS
921 // // and zero the column
922 // for (unsigned int i=0; i<this->m(); i++)
923 // {
924 // rhs(i) -= ((*this)(i,jv))*val;
925 // (*this)(i,jv) = 0.;
926 // }
927 
928 // // zero the row
929 // for (unsigned int j=0; j<this->n(); j++)
930 // (*this)(iv,j) = 0.;
931 
932 // (*this)(iv,jv) = 1.;
933 // rhs(iv) = val;
934 
935 // }
936 
937 
938 
939 
940 } // namespace libMesh
941 
942 
943 
944 
945 #endif // LIBMESH_DENSE_MATRIX_H

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

Hosted By:
SourceForge.net Logo