4 #ifndef DUNE_DENSEMATRIX_HH
5 #define DUNE_DENSEMATRIX_HH
23 template<
typename M>
class DenseMatrix;
42 static typename V::size_type
size(
const V & v) {
return v.size(); }
45 template<
class K,
int N>
48 typedef FieldVector<K,N> V;
49 static typename V::size_type
size(
const V & v) {
return N; }
68 template<
typename M,
typename T>
71 DUNE_THROW(
NotImplemented,
"You need to specialise the method istl_assign_to_fmatrix(DenseMatrix<M>& f, const T& t) "
72 <<
"(with M being " << className<M>() <<
") "
73 <<
"for T == " << className<T>() <<
"!");
79 struct DenseMatrixAssigner
81 template<
typename M,
typename T>
82 static void assign(DenseMatrix<M>& fm,
const T&
t)
91 struct DenseMatrixAssigner<true>
93 template<
typename M,
typename T>
94 static void assign(DenseMatrix<M>& fm,
const T&
t)
96 fm =
static_cast<const typename DenseMatVecTraits<M>::value_type
>(
t);
114 template<
typename MAT>
120 MAT & asImp() {
return static_cast<MAT&
>(*this); }
121 const MAT & asImp()
const {
return static_cast<const MAT&
>(*this); }
161 return asImp().mat_access(i);
166 return asImp().mat_access(i);
257 DenseMatrixAssigner<Conversion<T,field_type>::exists>::assign(*
this, t);
263 template <
class Other>
272 template <
class Other>
297 template <
class Other>
301 (*
this)[ i ].
axpy( k, y[ i ] );
306 template <
class Other>
310 if ((*
this)[i]!=y[i])
315 template <
class Other>
325 template<
class X,
class Y>
326 void mv (
const X& x, Y& y)
const
328 #ifdef DUNE_FMatrix_WITH_CHECKING
336 y[i] += (*
this)[i][j] * x[j];
341 template<
class X,
class Y >
342 void mtv (
const X &x, Y &y )
const
344 #ifdef DUNE_FMatrix_WITH_CHECKING
358 y[ i ] += (*
this)[ j ][ i ] * x[ j ];
363 template<
class X,
class Y>
364 void umv (
const X& x, Y& y)
const
366 #ifdef DUNE_FMatrix_WITH_CHECKING
368 DUNE_THROW(
FMatrixError,
"y += A x -- index out of range (sizes: x: " << x.N() <<
", y: " << y.N() <<
", A: " << this->
N() <<
" x " << this->
M() <<
")" << std::endl);
370 DUNE_THROW(
FMatrixError,
"y += A x -- index out of range (sizes: x: " << x.N() <<
", y: " << y.N() <<
", A: " << this->
N() <<
" x " << this->
M() <<
")" << std::endl);
374 y[i] += (*
this)[i][j] * x[j];
378 template<
class X,
class Y>
379 void umtv (
const X& x, Y& y)
const
381 #ifdef DUNE_FMatrix_WITH_CHECKING
388 y[j] += (*
this)[i][j]*x[i];
392 template<
class X,
class Y>
393 void umhv (
const X& x, Y& y)
const
395 #ifdef DUNE_FMatrix_WITH_CHECKING
406 template<
class X,
class Y>
407 void mmv (
const X& x, Y& y)
const
409 #ifdef DUNE_FMatrix_WITH_CHECKING
415 y[i] -= (*
this)[i][j] * x[j];
419 template<
class X,
class Y>
420 void mmtv (
const X& x, Y& y)
const
422 #ifdef DUNE_FMatrix_WITH_CHECKING
429 y[j] -= (*
this)[i][j]*x[i];
433 template<
class X,
class Y>
434 void mmhv (
const X& x, Y& y)
const
436 #ifdef DUNE_FMatrix_WITH_CHECKING
447 template<
class X,
class Y>
450 #ifdef DUNE_FMatrix_WITH_CHECKING
456 y[i] += alpha * (*
this)[i][j] * x[j];
460 template<
class X,
class Y>
463 #ifdef DUNE_FMatrix_WITH_CHECKING
470 y[j] += alpha*(*
this)[i][j]*x[i];
474 template<
class X,
class Y>
477 #ifdef DUNE_FMatrix_WITH_CHECKING
493 for (
size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
494 return fvmeta::sqrt(sum);
501 for (
size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
513 for (it = it + 1; it !=
end(); ++it)
514 max = std::max(max, it->one_norm());
527 for (it = it + 1; it !=
end(); ++it)
528 max = std::max(max, it->one_norm_real());
540 void solve (V& x,
const V& b)
const;
552 template<
typename M2>
562 (*
this)[i][j] += M[i][k]*C[k][j];
569 template<
typename M2>
579 (*
this)[i][j] += C[i][k]*M[k][j];
595 C[i][j] += M[i][k]*(*
this)[k][j];
603 FieldMatrix<K,rows,l> rightmultiplyany (
const FieldMatrix<K,cols,l>& M)
const
605 FieldMatrix<K,rows,l> C;
611 C[i][j] += (*
this)[i][k]*M[k][j];
635 return asImp().mat_rows();
641 return asImp().mat_cols();
649 #ifdef DUNE_FMatrix_WITH_CHECKING
661 ElimPivot(std::vector<size_type> & pivot);
663 void swap(
int i,
int j);
666 void operator()(
const T&,
int k,
int i)
669 std::vector<size_type> & pivot_;
677 void swap(
int i,
int j);
679 void operator()(
const typename V::field_type& factor,
int k,
int i);
689 void swap(
int i,
int j)
692 void operator()(
const field_type&,
int k,
int i)
700 void luDecomposition(DenseMatrix<MAT>& A, Func func)
const;
704 template<
typename MAT>
705 DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<size_type> & pivot)
708 typedef typename std::vector<size_type>::size_type size_type;
709 for(size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
712 template<
typename MAT>
713 void DenseMatrix<MAT>::ElimPivot::swap(
int i,
int j)
718 template<
typename MAT>
720 DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
724 template<
typename MAT>
726 void DenseMatrix<MAT>::Elim<V>::swap(
int i,
int j)
728 std::swap((*rhs_)[i], (*rhs_)[j]);
731 template<
typename MAT>
733 void DenseMatrix<MAT>::
734 Elim<V>::operator()(
const typename V::field_type& factor,
int k,
int i)
736 (*rhs_)[k] -= factor*(*rhs_)[i];
738 template<
typename MAT>
739 template<
typename Func>
740 inline void DenseMatrix<MAT>::luDecomposition(DenseMatrix<MAT>& A, Func func)
const
742 typedef typename FieldTraits<value_type>::real_type
744 real_type norm = A.infinity_norm_real();
745 real_type pivthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::pivoting_limit() );
746 real_type singthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::singular_limit() );
749 for (size_type i=0; i<rows(); i++)
751 typename FieldTraits<value_type>::real_type pivmax=fvmeta::absreal(A[i][i]);
758 typename FieldTraits<value_type>::real_type abs(0.0);
759 for (size_type k=i+1; k<rows(); k++)
760 if ((abs=fvmeta::absreal(A[k][i]))>pivmax)
762 pivmax = abs; imax = k;
766 for (size_type j=0; j<rows(); j++)
767 std::swap(A[i][j],A[imax][j]);
773 if (pivmax<singthres)
774 DUNE_THROW(FMatrixError,
"matrix is singular");
777 for (size_type k=i+1; k<rows(); k++)
779 field_type factor = A[k][i]/A[i][i];
781 for (size_type j=i+1; j<rows(); j++)
782 A[k][j] -= factor*A[i][j];
788 template<
typename MAT>
790 inline void DenseMatrix<MAT>::solve(V& x,
const V& b)
const
794 DUNE_THROW(FMatrixError,
"Can't solve for a " << rows() <<
"x" << cols() <<
" matrix!");
798 #ifdef DUNE_FMatrix_WITH_CHECKING
799 if (fvmeta::absreal((*
this)[0][0])<FMatrixPrecision<>::absolute_limit())
800 DUNE_THROW(FMatrixError,
"matrix is singular");
802 x[0] = b[0]/(*this)[0][0];
805 else if (rows()==2) {
807 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
808 #ifdef DUNE_FMatrix_WITH_CHECKING
809 if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
810 DUNE_THROW(FMatrixError,
"matrix is singular");
814 x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
815 x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
818 else if (rows()==3) {
820 field_type d = determinant();
821 #ifdef DUNE_FMatrix_WITH_CHECKING
822 if (fvmeta::absreal(d)<FMatrixPrecision<>::absolute_limit())
823 DUNE_THROW(FMatrixError,
"matrix is singular");
826 x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
827 - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
828 + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
830 x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
831 - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
832 + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
834 x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
835 - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
836 + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
846 luDecomposition(A, elim);
849 for(
int i=rows()-1; i>=0; i--){
850 for (size_type j=i+1; j<rows(); j++)
851 rhs[i] -= A[i][j]*x[j];
852 x[i] = rhs[i]/A[i][i];
857 template<
typename MAT>
858 inline void DenseMatrix<MAT>::invert()
862 DUNE_THROW(FMatrixError,
"Can't invert a " << rows() <<
"x" << cols() <<
" matrix!");
866 #ifdef DUNE_FMatrix_WITH_CHECKING
867 if (fvmeta::absreal((*
this)[0][0])<FMatrixPrecision<>::absolute_limit())
868 DUNE_THROW(FMatrixError,
"matrix is singular");
870 (*this)[0][0] = 1.0/(*this)[0][0];
873 else if (rows()==2) {
875 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
876 #ifdef DUNE_FMatrix_WITH_CHECKING
877 if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
878 DUNE_THROW(FMatrixError,
"matrix is singular");
882 field_type temp=(*this)[0][0];
883 (*this)[0][0] = (*this)[1][1]*detinv;
884 (*this)[0][1] = -(*this)[0][1]*detinv;
885 (*this)[1][0] = -(*this)[1][0]*detinv;
886 (*this)[1][1] = temp*detinv;
892 std::vector<size_type> pivot(rows());
893 luDecomposition(A, ElimPivot(pivot));
894 DenseMatrix<MAT>& L=A;
895 DenseMatrix<MAT>& U=A;
900 for(size_type i=0; i<rows(); ++i)
904 for (size_type i=0; i<rows(); i++)
905 for (size_type j=0; j<i; j++)
906 for (size_type k=0; k<rows(); k++)
907 (*
this)[i][k] -= L[i][j]*(*this)[j][k];
910 for (size_type i=rows(); i>0;){
912 for (size_type k=0; k<rows(); k++){
913 for (size_type j=i+1; j<rows(); j++)
914 (*
this)[i][k] -= U[i][j]*(*this)[j][k];
915 (*this)[i][k] /= U[i][i];
919 for(size_type i=rows(); i>0; ){
922 for(size_type j=0; j<rows(); ++j)
923 std::swap((*
this)[j][pivot[i]], (*this)[j][i]);
929 template<
typename MAT>
930 inline typename DenseMatrix<MAT>::field_type
931 DenseMatrix<MAT>::determinant()
const
935 DUNE_THROW(FMatrixError,
"There is no determinant for a " << rows() <<
"x" << cols() <<
" matrix!");
938 return (*
this)[0][0];
941 return (*
this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
945 field_type t4 = (*this)[0][0] * (*this)[1][1];
946 field_type t6 = (*this)[0][0] * (*this)[1][2];
947 field_type t8 = (*this)[0][1] * (*this)[1][0];
948 field_type t10 = (*this)[0][2] * (*this)[1][0];
949 field_type t12 = (*this)[0][1] * (*this)[2][0];
950 field_type t14 = (*this)[0][2] * (*this)[2][0];
952 return (t4*(*
this)[2][2]-t6*(*
this)[2][1]-t8*(*
this)[2][2]+
953 t10*(*
this)[2][1]+t12*(*
this)[1][2]-t14*(*
this)[1][1]);
961 luDecomposition(A, ElimDet(det));
963 catch (FMatrixError&)
967 for (size_type i = 0; i < rows(); ++i)
974 namespace DenseMatrixHelp {
977 template <
typename K>
980 inverse[0][0] = 1.0/matrix[0][0];
985 template <
typename K>
993 template <
typename K>
997 field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
998 field_type det_1 = 1.0/det;
999 inverse[0][0] = matrix[1][1] * det_1;
1000 inverse[0][1] = - matrix[0][1] * det_1;
1001 inverse[1][0] = - matrix[1][0] * det_1;
1002 inverse[1][1] = matrix[0][0] * det_1;
1008 template <
typename K>
1012 field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
1013 field_type det_1 = 1.0/det;
1014 inverse[0][0] = matrix[1][1] * det_1;
1015 inverse[1][0] = - matrix[0][1] * det_1;
1016 inverse[0][1] = - matrix[1][0] * det_1;
1017 inverse[1][1] = matrix[0][0] * det_1;
1022 template <
typename K>
1026 field_type t4 = matrix[0][0] * matrix[1][1];
1027 field_type t6 = matrix[0][0] * matrix[1][2];
1028 field_type t8 = matrix[0][1] * matrix[1][0];
1029 field_type t10 = matrix[0][2] * matrix[1][0];
1030 field_type t12 = matrix[0][1] * matrix[2][0];
1031 field_type t14 = matrix[0][2] * matrix[2][0];
1033 field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
1034 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
1035 field_type t17 = 1.0/det;
1037 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
1038 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
1039 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
1040 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
1041 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
1042 inverse[1][2] = -(t6-t10) * t17;
1043 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
1044 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
1045 inverse[2][2] = (t4-t8) * t17;
1051 template <
typename K>
1055 field_type t4 = matrix[0][0] * matrix[1][1];
1056 field_type t6 = matrix[0][0] * matrix[1][2];
1057 field_type t8 = matrix[0][1] * matrix[1][0];
1058 field_type t10 = matrix[0][2] * matrix[1][0];
1059 field_type t12 = matrix[0][1] * matrix[2][0];
1060 field_type t14 = matrix[0][2] * matrix[2][0];
1062 field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
1063 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
1064 field_type t17 = 1.0/det;
1066 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
1067 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
1068 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
1069 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
1070 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
1071 inverse[2][1] = -(t6-t10) * t17;
1072 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
1073 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
1074 inverse[2][2] = (t4-t8) * t17;
1080 template<
class K,
int m,
int n,
int p >
1087 for( size_type i = 0; i < m; ++i )
1089 for( size_type j = 0; j < p; ++j )
1091 ret[ i ][ j ] = K( 0 );
1092 for( size_type k = 0; k < n; ++k )
1093 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
1099 template <
typename K,
int rows,
int cols>
1104 for(size_type i=0; i<cols(); i++)
1105 for(size_type j=0; j<cols(); j++)
1108 for(size_type k=0; k<rows(); k++)
1109 ret[i][j]+=matrix[k][i]*matrix[k][j];
1115 template <
typename MAT,
typename V1,
typename V2>
1119 assert(ret.
size() == matrix.
rows());
1122 for(size_type i=0; i<matrix.
rows(); ++i)
1125 for(size_type j=0; j<matrix.
cols(); ++j)
1127 ret[i] += matrix[i][j]*x[j];
1134 template <
typename K,
int rows,
int cols>
1139 for(size_type i=0; i<cols(); ++i)
1142 for(size_type j=0; j<rows(); ++j)
1143 ret[i] += matrix[j][i]*x[j];
1148 template <
typename K,
int rows,
int cols>
1149 static inline FieldVector<K,rows>
mult(
const FieldMatrix<K,rows,cols> &matrix,
const FieldVector<K,cols> & x)
1151 FieldVector<K,rows> ret;
1157 template <
typename K,
int rows,
int cols>
1158 static inline FieldVector<K,cols>
multTransposed(
const FieldMatrix<K,rows,cols> &matrix,
const FieldVector<K,rows> & x)
1160 FieldVector<K,cols> ret;
1169 template<
typename MAT>
1170 std::ostream& operator<< (std::ostream& s, const DenseMatrix<MAT>& a)
1173 s << a[i] << std::endl;