3 #ifndef DUNE_DENSEMATRIX_HH
4 #define DUNE_DENSEMATRIX_HH
41 static typename V::size_type
size(
const V & v) {
return v.size(); }
44 template<
class K,
int N>
47 typedef FieldVector<K,N> V;
48 static typename V::size_type
size(
const V & v) {
return N; }
70 template<
class DenseMatrix,
class RHS >
75 template<
class DenseMatrix,
class K,
int N,
int M >
78 for(
int i = 0; i < N; ++i )
79 for(
int j = 0; j < M; ++j )
80 denseMatrix[ i ][ j ] = values[ i ][ j ];
88 template<
class DenseMatrix,
class RHS,
89 bool primitive = Conversion< RHS, typename DenseMatrix::field_type >::exists >
90 class DenseMatrixAssignerImplementation;
92 template<
class DenseMatrix,
class RHS >
93 class DenseMatrixAssignerImplementation< DenseMatrix, RHS, true >
96 static void apply ( DenseMatrix &denseMatrix,
const RHS &rhs )
99 std::fill( denseMatrix.begin(), denseMatrix.end(),
static_cast< field_type
>( rhs ) );
103 template<
class DenseMatrix,
class RHS >
104 class DenseMatrixAssignerImplementation< DenseMatrix, RHS, false >
106 template<
class M,
class T>
107 struct have_istl_assign_to_fmatrix
109 struct yes {
char dummy[ 1 ]; };
110 struct no {
char dummy[ 2 ]; };
121 static const bool v =
sizeof( test< const T >( 0 ) ) ==
sizeof( yes );
124 template< class M, class T, bool = have_istl_assign_to_fmatrix< M, T >::v >
125 struct DefaultImplementation;
128 template<
class M,
class T >
129 struct DefaultImplementation< M, T, true >
131 static void apply ( M &m,
const T &
t )
138 template<
class M,
class T >
139 struct DefaultImplementation< M, T, false >
141 static void apply ( M &m,
const T &
t )
143 static_assert( (Conversion< const T, const M >::exists),
"No template specialization of DenseMatrixAssigner found" );
144 m =
static_cast< const M &
>(
t );
149 static void apply ( DenseMatrix &denseMatrix,
const RHS &rhs )
151 DefaultImplementation< DenseMatrix, RHS >::apply( denseMatrix, rhs );
158 template<
class DenseMatrix,
class RHS >
159 struct DenseMatrixAssigner
161 static void apply ( DenseMatrix &denseMatrix,
const RHS &rhs )
163 DenseMatrixAssignerImplementation< DenseMatrix, RHS >::apply( denseMatrix, rhs );
166 #endif // #ifndef DOXYGEN
183 template<
typename MAT>
189 MAT & asImp() {
return static_cast<MAT&
>(*this); }
190 const MAT & asImp()
const {
return static_cast<const MAT&
>(*this); }
230 return asImp().mat_access(i);
235 return asImp().mat_access(i);
252 typedef typename remove_reference<row_reference>::type::Iterator
ColIterator;
287 typedef typename remove_reference<const_row_reference>::type::ConstIterator
ConstColIterator;
296 ConstIterator
end ()
const
317 template<
class RHS >
327 template <
class Other>
330 for (size_type i=0; i<
rows(); i++)
336 template <
class Other>
339 for (size_type i=0; i<
rows(); i++)
347 for (size_type i=0; i<
rows(); i++)
355 for (size_type i=0; i<
rows(); i++)
361 template <
class Other>
364 for( size_type i = 0; i <
rows(); ++i )
365 (*
this)[ i ].
axpy( k, y[ i ] );
370 template <
class Other>
373 for (size_type i=0; i<
rows(); i++)
374 if ((*
this)[i]!=y[i])
379 template <
class Other>
389 template<
class X,
class Y>
390 void mv (
const X& x, Y& y)
const
392 #ifdef DUNE_FMatrix_WITH_CHECKING
393 assert( (
void*)(&x) != (
void*)(&y) );
397 for (size_type i=0; i<
rows(); ++i)
400 for (size_type j=0; j<
cols(); j++)
401 y[i] += (*
this)[i][j] * x[j];
406 template<
class X,
class Y >
407 void mtv (
const X &x, Y &y )
const
409 #ifdef DUNE_FMatrix_WITH_CHECKING
410 assert( (
void*)(&x) != (
void*)(&y) );
416 for( size_type i = 0; i <
cols(); ++i )
419 for( size_type j = 0; j <
rows(); ++j )
420 y[ i ] += (*
this)[ j ][ i ] * x[ j ];
425 template<
class X,
class Y>
426 void umv (
const X& x, Y& y)
const
428 #ifdef DUNE_FMatrix_WITH_CHECKING
430 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);
432 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);
434 for (size_type i=0; i<
rows(); i++)
435 for (size_type j=0; j<
cols(); j++)
436 y[i] += (*
this)[i][j] * x[j];
440 template<
class X,
class Y>
441 void umtv (
const X& x, Y& y)
const
443 #ifdef DUNE_FMatrix_WITH_CHECKING
448 for (size_type i=0; i<
rows(); i++)
449 for (size_type j=0; j<
cols(); j++)
450 y[j] += (*
this)[i][j]*x[i];
454 template<
class X,
class Y>
455 void umhv (
const X& x, Y& y)
const
457 #ifdef DUNE_FMatrix_WITH_CHECKING
462 for (size_type i=0; i<
rows(); i++)
463 for (size_type j=0; j<
cols(); j++)
468 template<
class X,
class Y>
469 void mmv (
const X& x, Y& y)
const
471 #ifdef DUNE_FMatrix_WITH_CHECKING
475 for (size_type i=0; i<
rows(); i++)
476 for (size_type j=0; j<
cols(); j++)
477 y[i] -= (*
this)[i][j] * x[j];
481 template<
class X,
class Y>
482 void mmtv (
const X& x, Y& y)
const
484 #ifdef DUNE_FMatrix_WITH_CHECKING
489 for (size_type i=0; i<
rows(); i++)
490 for (size_type j=0; j<
cols(); j++)
491 y[j] -= (*
this)[i][j]*x[i];
495 template<
class X,
class Y>
496 void mmhv (
const X& x, Y& y)
const
498 #ifdef DUNE_FMatrix_WITH_CHECKING
503 for (size_type i=0; i<
rows(); i++)
504 for (size_type j=0; j<
cols(); j++)
509 template<
class X,
class Y>
510 void usmv (
const field_type& alpha,
const X& x, Y& y)
const
512 #ifdef DUNE_FMatrix_WITH_CHECKING
516 for (size_type i=0; i<
rows(); i++)
517 for (size_type j=0; j<
cols(); j++)
518 y[i] += alpha * (*
this)[i][j] * x[j];
522 template<
class X,
class Y>
523 void usmtv (
const field_type& alpha,
const X& x, Y& y)
const
525 #ifdef DUNE_FMatrix_WITH_CHECKING
530 for (size_type i=0; i<
rows(); i++)
531 for (size_type j=0; j<
cols(); j++)
532 y[j] += alpha*(*
this)[i][j]*x[i];
536 template<
class X,
class Y>
537 void usmhv (
const field_type& alpha,
const X& x, Y& y)
const
539 #ifdef DUNE_FMatrix_WITH_CHECKING
544 for (size_type i=0; i<
rows(); i++)
545 for (size_type j=0; j<
cols(); j++)
555 for (size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
556 return fvmeta::sqrt(sum);
563 for (size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
573 ConstIterator it =
begin();
574 typename remove_const< typename FieldTraits<value_type>::real_type >::type max = (*it).one_norm();
575 for (it = it + 1; it !=
end(); ++it)
576 max = std::max(max, (*it).one_norm());
587 ConstIterator it =
begin();
588 typename remove_const< typename FieldTraits<value_type>::real_type >::type max = (*it).one_norm_real();
589 for (it = it + 1; it !=
end(); ++it)
590 max = std::max(max, (*it).one_norm_real());
602 void solve (V& x,
const V& b)
const;
614 template<
typename M2>
620 for (size_type i=0; i<
rows(); i++)
621 for (size_type j=0; j<
cols(); j++) {
623 for (size_type k=0; k<
rows(); k++)
624 (*
this)[i][j] += M[i][k]*C[k][j];
631 template<
typename M2>
637 for (size_type i=0; i<
rows(); i++)
638 for (size_type j=0; j<
cols(); j++) {
640 for (size_type k=0; k<
cols(); k++)
641 (*
this)[i][j] += C[i][k]*M[k][j];
653 for (size_type i=0; i<l; i++) {
654 for (size_type j=0; j<
cols(); j++) {
656 for (size_type k=0; k<
rows(); k++)
657 C[i][j] += M[i][k]*(*
this)[k][j];
665 FieldMatrix<K,rows,l> rightmultiplyany (
const FieldMatrix<K,cols,l>& M)
const
667 FieldMatrix<K,rows,l> C;
669 for (size_type i=0; i<
rows(); i++) {
670 for (size_type j=0; j<l; j++) {
672 for (size_type k=0; k<
cols(); k++)
673 C[i][j] += (*
this)[i][k]*M[k][j];
697 return asImp().mat_rows();
703 return asImp().mat_cols();
709 bool exists (size_type i, size_type j)
const
711 #ifdef DUNE_FMatrix_WITH_CHECKING
726 ElimPivot(std::vector<size_type> & pivot);
728 void swap(
int i,
int j);
731 void operator()(
const T&,
int,
int)
734 std::vector<size_type> & pivot_;
742 void swap(
int i,
int j);
744 void operator()(
const typename V::field_type& factor,
int k,
int i);
751 ElimDet(field_type&
sign) : sign_(sign)
757 void operator()(
const field_type&,
int,
int)
765 void luDecomposition(DenseMatrix<MAT>& A, Func func)
const;
769 template<
typename MAT>
770 DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<size_type> & pivot)
773 typedef typename std::vector<size_type>::size_type size_type;
774 for(size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
777 template<
typename MAT>
778 void DenseMatrix<MAT>::ElimPivot::swap(
int i,
int j)
783 template<
typename MAT>
785 DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
789 template<
typename MAT>
791 void DenseMatrix<MAT>::Elim<V>::swap(
int i,
int j)
793 std::swap((*rhs_)[i], (*rhs_)[j]);
796 template<
typename MAT>
798 void DenseMatrix<MAT>::
799 Elim<V>::operator()(
const typename V::field_type& factor,
int k,
int i)
801 (*rhs_)[k] -= factor*(*rhs_)[i];
803 template<
typename MAT>
804 template<
typename Func>
805 inline void DenseMatrix<MAT>::luDecomposition(DenseMatrix<MAT>& A, Func func)
const
807 typedef typename FieldTraits<value_type>::real_type
809 real_type norm = A.infinity_norm_real();
810 real_type pivthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::pivoting_limit() );
811 real_type singthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::singular_limit() );
814 for (size_type i=0; i<rows(); i++)
816 typename FieldTraits<value_type>::real_type pivmax=fvmeta::absreal(A[i][i]);
823 typename FieldTraits<value_type>::real_type abs(0.0);
824 for (size_type k=i+1; k<rows(); k++)
825 if ((abs=fvmeta::absreal(A[k][i]))>pivmax)
827 pivmax = abs; imax = k;
831 for (size_type j=0; j<rows(); j++)
832 std::swap(A[i][j],A[imax][j]);
838 if (pivmax<singthres)
839 DUNE_THROW(FMatrixError,
"matrix is singular");
842 for (size_type k=i+1; k<rows(); k++)
844 field_type factor = A[k][i]/A[i][i];
846 for (size_type j=i+1; j<rows(); j++)
847 A[k][j] -= factor*A[i][j];
853 template<
typename MAT>
855 inline void DenseMatrix<MAT>::solve(V& x,
const V& b)
const
859 DUNE_THROW(FMatrixError,
"Can't solve for a " << rows() <<
"x" << cols() <<
" matrix!");
863 #ifdef DUNE_FMatrix_WITH_CHECKING
864 if (fvmeta::absreal((*
this)[0][0])<FMatrixPrecision<>::absolute_limit())
865 DUNE_THROW(FMatrixError,
"matrix is singular");
867 x[0] = b[0]/(*this)[0][0];
870 else if (rows()==2) {
872 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
873 #ifdef DUNE_FMatrix_WITH_CHECKING
874 if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
875 DUNE_THROW(FMatrixError,
"matrix is singular");
879 x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
880 x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
883 else if (rows()==3) {
885 field_type d = determinant();
886 #ifdef DUNE_FMatrix_WITH_CHECKING
887 if (fvmeta::absreal(d)<FMatrixPrecision<>::absolute_limit())
888 DUNE_THROW(FMatrixError,
"matrix is singular");
891 x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
892 - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
893 + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
895 x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
896 - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
897 + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
899 x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
900 - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
901 + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
911 luDecomposition(A, elim);
914 for(
int i=rows()-1; i>=0; i--) {
915 for (size_type j=i+1; j<rows(); j++)
916 rhs[i] -= A[i][j]*x[j];
917 x[i] = rhs[i]/A[i][i];
922 template<
typename MAT>
923 inline void DenseMatrix<MAT>::invert()
927 DUNE_THROW(FMatrixError,
"Can't invert a " << rows() <<
"x" << cols() <<
" matrix!");
931 #ifdef DUNE_FMatrix_WITH_CHECKING
932 if (fvmeta::absreal((*
this)[0][0])<FMatrixPrecision<>::absolute_limit())
933 DUNE_THROW(FMatrixError,
"matrix is singular");
935 (*this)[0][0] = field_type( 1 ) / (*this)[0][0];
938 else if (rows()==2) {
940 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
941 #ifdef DUNE_FMatrix_WITH_CHECKING
942 if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
943 DUNE_THROW(FMatrixError,
"matrix is singular");
945 detinv = field_type( 1 ) / detinv;
947 field_type temp=(*this)[0][0];
948 (*this)[0][0] = (*this)[1][1]*detinv;
949 (*this)[0][1] = -(*this)[0][1]*detinv;
950 (*this)[1][0] = -(*this)[1][0]*detinv;
951 (*this)[1][1] = temp*detinv;
957 std::vector<size_type> pivot(rows());
958 luDecomposition(A, ElimPivot(pivot));
959 DenseMatrix<MAT>& L=A;
960 DenseMatrix<MAT>& U=A;
965 for(size_type i=0; i<rows(); ++i)
969 for (size_type i=0; i<rows(); i++)
970 for (size_type j=0; j<i; j++)
971 for (size_type k=0; k<rows(); k++)
972 (*
this)[i][k] -= L[i][j]*(*this)[j][k];
975 for (size_type i=rows(); i>0;) {
977 for (size_type k=0; k<rows(); k++) {
978 for (size_type j=i+1; j<rows(); j++)
979 (*
this)[i][k] -= U[i][j]*(*this)[j][k];
980 (*this)[i][k] /= U[i][i];
984 for(size_type i=rows(); i>0; ) {
987 for(size_type j=0; j<rows(); ++j)
988 std::swap((*
this)[j][pivot[i]], (*this)[j][i]);
994 template<
typename MAT>
995 inline typename DenseMatrix<MAT>::field_type
996 DenseMatrix<MAT>::determinant()
const
1000 DUNE_THROW(FMatrixError,
"There is no determinant for a " << rows() <<
"x" << cols() <<
" matrix!");
1003 return (*
this)[0][0];
1006 return (*
this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
1010 field_type t4 = (*this)[0][0] * (*this)[1][1];
1011 field_type t6 = (*this)[0][0] * (*this)[1][2];
1012 field_type t8 = (*this)[0][1] * (*this)[1][0];
1013 field_type t10 = (*this)[0][2] * (*this)[1][0];
1014 field_type t12 = (*this)[0][1] * (*this)[2][0];
1015 field_type t14 = (*this)[0][2] * (*this)[2][0];
1017 return (t4*(*
this)[2][2]-t6*(*
this)[2][1]-t8*(*
this)[2][2]+
1018 t10*(*
this)[2][1]+t12*(*
this)[1][2]-t14*(*
this)[1][1]);
1026 luDecomposition(A, ElimDet(det));
1028 catch (FMatrixError&)
1032 for (size_type i = 0; i < rows(); ++i)
1039 namespace DenseMatrixHelp {
1041 template <
typename K>
1045 inverse[0][0] = 1.0/matrix[0][0];
1046 return matrix[0][0];
1050 template <
typename K>
1058 template <
typename K>
1062 field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
1063 field_type det_1 = 1.0/det;
1064 inverse[0][0] = matrix[1][1] * det_1;
1065 inverse[0][1] = - matrix[0][1] * det_1;
1066 inverse[1][0] = - matrix[1][0] * det_1;
1067 inverse[1][1] = matrix[0][0] * det_1;
1073 template <
typename K>
1077 field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
1078 field_type det_1 = 1.0/det;
1079 inverse[0][0] = matrix[1][1] * det_1;
1080 inverse[1][0] = - matrix[0][1] * det_1;
1081 inverse[0][1] = - matrix[1][0] * det_1;
1082 inverse[1][1] = matrix[0][0] * det_1;
1087 template <
typename K>
1091 field_type t4 = matrix[0][0] * matrix[1][1];
1092 field_type t6 = matrix[0][0] * matrix[1][2];
1093 field_type t8 = matrix[0][1] * matrix[1][0];
1094 field_type t10 = matrix[0][2] * matrix[1][0];
1095 field_type t12 = matrix[0][1] * matrix[2][0];
1096 field_type t14 = matrix[0][2] * matrix[2][0];
1098 field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
1099 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
1100 field_type t17 = 1.0/det;
1102 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
1103 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
1104 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
1105 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
1106 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
1107 inverse[1][2] = -(t6-t10) * t17;
1108 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
1109 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
1110 inverse[2][2] = (t4-t8) * t17;
1116 template <
typename K>
1120 field_type t4 = matrix[0][0] * matrix[1][1];
1121 field_type t6 = matrix[0][0] * matrix[1][2];
1122 field_type t8 = matrix[0][1] * matrix[1][0];
1123 field_type t10 = matrix[0][2] * matrix[1][0];
1124 field_type t12 = matrix[0][1] * matrix[2][0];
1125 field_type t14 = matrix[0][2] * matrix[2][0];
1127 field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
1128 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
1129 field_type t17 = 1.0/det;
1131 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
1132 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
1133 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
1134 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
1135 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
1136 inverse[2][1] = -(t6-t10) * t17;
1137 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
1138 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
1139 inverse[2][2] = (t4-t8) * t17;
1145 template<
class K,
int m,
int n,
int p >
1152 for( size_type i = 0; i < m; ++i )
1154 for( size_type j = 0; j < p; ++j )
1156 ret[ i ][ j ] = K( 0 );
1157 for( size_type k = 0; k < n; ++k )
1158 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
1164 template <
typename K,
int rows,
int cols>
1169 for(size_type i=0; i<cols(); i++)
1170 for(size_type j=0; j<cols(); j++)
1173 for(size_type k=0; k<rows(); k++)
1174 ret[i][j]+=matrix[k][i]*matrix[k][j];
1180 template <
typename MAT,
typename V1,
typename V2>
1184 assert(ret.
size() == matrix.
rows());
1187 for(size_type i=0; i<matrix.
rows(); ++i)
1190 for(size_type j=0; j<matrix.
cols(); ++j)
1192 ret[i] += matrix[i][j]*x[j];
1198 template <
typename K,
int rows,
int cols>
1204 for(size_type i=0; i<cols(); ++i)
1207 for(size_type j=0; j<rows(); ++j)
1208 ret[i] += matrix[j][i]*x[j];
1213 template <
typename K,
int rows,
int cols>
1214 static inline FieldVector<K,rows>
mult(
const FieldMatrix<K,rows,cols> &matrix,
const FieldVector<K,cols> & x)
1216 FieldVector<K,rows> ret;
1222 template <
typename K,
int rows,
int cols>
1223 static inline FieldVector<K,cols>
multTransposed(
const FieldMatrix<K,rows,cols> &matrix,
const FieldVector<K,rows> & x)
1225 FieldVector<K,cols> ret;
1234 template<
typename MAT>
1235 std::ostream& operator<< (std::ostream& s, const DenseMatrix<MAT>& a)
1238 s << a[i] << std::endl;
DenseMatrix & axpy(const field_type &k, const DenseMatrix< Other > &y)
vector space axpy operation (*this += k y)
Definition: densematrix.hh:362
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: densematrix.hh:552
DUNE_CONSTEXPR size_type size() const
Definition: fvector.hh:160
ConstIterator end() const
end iterator
Definition: densematrix.hh:296
ConstIterator beforeBegin() const
Definition: densematrix.hh:310
Traits::value_type value_type
export the type representing the field
Definition: densematrix.hh:199
ConstIterator begin() const
begin iterator
Definition: densematrix.hh:290
Iterator RowIterator
rename the iterators for easier access
Definition: densematrix.hh:250
A dense n x m matrix.
Definition: densematrix.hh:22
size_type size() const
size method
Definition: densevector.hh:296
const FieldTraits< typename DenseMatVecTraits< M >::value_type >::real_type real_type
Definition: densematrix.hh:28
Iterator begin()
begin iterator
Definition: densematrix.hh:255
void umhv(const X &x, Y &y) const
y += A^H x
Definition: densematrix.hh:455
Traits::const_row_reference const_row_reference
The type used to represent a reference to a constant row (usually const row_type &) ...
Definition: densematrix.hh:217
Definition: ftraits.hh:22
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:632
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:103
void istl_assign_to_fmatrix(DenseMatrix &denseMatrix, const K(&values)[M][N])
Definition: densematrix.hh:76
DenseIterator< const DenseMatrix, const row_type, const_row_reference > ConstIterator
Iterator class for sequential access.
Definition: densematrix.hh:281
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:119
void umtv(const X &x, Y &y) const
y += A^T x
Definition: densematrix.hh:441
static void multAssignTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x, FieldVector< K, cols > &ret)
calculates ret = matrix^T * x
Definition: fmatrix.hh:480
Iterator end()
end iterator
Definition: densematrix.hh:261
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:214
DenseIterator< DenseMatrix, row_type, row_reference > Iterator
Iterator class for sequential access.
Definition: densematrix.hh:246
A few common exception classes.
DenseMatrix & operator-=(const DenseMatrix< Other > &y)
vector space subtraction
Definition: densematrix.hh:337
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition: densevector.hh:643
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:125
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition: densematrix.hh:615
T real_type
export the type representing the real type of the field
Definition: ftraits.hh:27
The number of block levels we contain. This is 1.
Definition: densematrix.hh:222
size_type N() const
number of rows
Definition: densematrix.hh:683
FieldTraits< value_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: densematrix.hh:568
void mv(const X &x, Y &y) const
y = A x
Definition: densematrix.hh:390
Definition: matvectraits.hh:29
DenseMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:353
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: densematrix.hh:560
Interface for a class of dense vectors over a given field.
Definition: densevector.hh:18
DenseMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:345
ConstIterator beforeEnd() const
Definition: densematrix.hh:303
bool operator==(const DenseMatrix< Other > &y) const
Binary matrix comparison.
Definition: densematrix.hh:371
size_type M() const
number of columns
Definition: densematrix.hh:689
static FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition: fmatrix.hh:494
Traits::value_type block_type
export the type representing the components
Definition: densematrix.hh:205
T t
Definition: alignment.hh:34
field_type determinant() const
calculates the determinant of this matrix
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: densematrix.hh:523
Iterator beforeEnd()
Definition: densematrix.hh:268
remove_reference< row_reference >::type::Iterator ColIterator
rename the iterators for easier access
Definition: densematrix.hh:252
void invert()
Compute inverse.
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:211
static void multMatrix(const FieldMatrix< K, m, n > &A, const FieldMatrix< K, n, p > &B, FieldMatrix< K, m, p > &ret)
calculates ret = A * B
Definition: fmatrix.hh:444
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
Dune namespace.
Definition: alignment.hh:9
Base::size_type size_type
Definition: fmatrix.hh:80
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: densematrix.hh:709
Default exception class for mathematical errors.
Definition: exceptions.hh:266
bool operator!=(const DenseMatrix< Other > &y) const
Binary matrix incomparison.
Definition: densematrix.hh:380
Implements a vector constructed from a given type representing a field and a compile-time given size...
A free function to provide the demangled class name of a given object or type as a string...
size_type cols() const
number of columns
Definition: densematrix.hh:701
void umv(const X &x, Y &y) const
y += A x
Definition: densematrix.hh:426
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: densematrix.hh:496
static FieldVector< K, cols > multTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x)
calculates ret = matrix^T * x
Definition: fmatrix.hh:503
Some useful basic math stuff.
static K invertMatrix_retTransposed(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition: fmatrix.hh:349
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: densematrix.hh:537
ConstIterator const_iterator
typedef for stl compliant access
Definition: densematrix.hh:283
vector space out of a tensor product of fields.
Definition: densematrix.hh:36
static void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition: fmatrix.hh:463
DenseMatrix & operator=(const RHS &rhs)
Definition: densematrix.hh:318
row_reference operator[](size_type i)
random access
Definition: densematrix.hh:228
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18
you have to specialize this structure for any type that should be assignable to a DenseMatrix ...
Definition: densematrix.hh:71
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: densematrix.hh:482
void solve(V &x, const V &b) const
Solve system A x = b.
DenseMatrix & operator+=(const DenseMatrix< Other > &y)
vector space addition
Definition: densematrix.hh:328
Iterator iterator
typedef for stl compliant access
Definition: densematrix.hh:248
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition: densematrix.hh:1181
void mmv(const X &x, Y &y) const
y -= A x
Definition: densematrix.hh:469
Various precision settings for calculations with FieldMatrix and FieldVector.
void mtv(const X &x, Y &y) const
y = A^T x
Definition: densematrix.hh:407
FieldTraits< value_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densematrix.hh:582
Traits::value_type field_type
export the type representing the field
Definition: densematrix.hh:202
void usmv(const field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition: densematrix.hh:510
const FieldTraits< typename DenseMatVecTraits< M >::value_type >::field_type field_type
Definition: densematrix.hh:27
static K invertMatrix(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition: fmatrix.hh:341
Iterator beforeBegin()
Definition: densematrix.hh:275
remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: densematrix.hh:287
size_type size() const
size method (number of rows)
Definition: densematrix.hh:239
A dense n x m matrix.
Definition: densematrix.hh:35
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: densematrix.hh:285
size_type rows() const
number of rows
Definition: densematrix.hh:695
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:208
Traits::derived_type derived_type
type of derived matrix class
Definition: densematrix.hh:196
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:171
FieldVector()
Constructor making default-initialized vector.
Definition: fvector.hh:108