4 #ifndef DUNE_DENSEMATRIX_HH
5 #define DUNE_DENSEMATRIX_HH
43 static typename V::size_type
size(
const V & v) {
return v.size(); }
46 template<
class K,
int N>
49 typedef FieldVector<K,N> V;
50 static typename V::size_type
size(
const V & v) {
return N; }
72 template<
class DenseMatrix,
class RHS >
77 template<
class DenseMatrix,
class K,
int N,
int M >
80 for(
int i = 0; i < N; ++i )
81 for(
int j = 0; j < M; ++j )
82 denseMatrix[ i ][ j ] = values[ i ][ j ];
90 template<
class DenseMatrix,
class RHS,
91 bool primitive = Conversion< RHS, typename DenseMatrix::field_type >::exists >
92 class DenseMatrixAssignerImplementation;
94 template<
class DenseMatrix,
class RHS >
95 class DenseMatrixAssignerImplementation< DenseMatrix, RHS, true >
98 static void apply ( DenseMatrix &denseMatrix,
const RHS &rhs )
102 const field_type value =
static_cast< field_type
>( rhs );
104 const std::size_t size = denseMatrix.size();
105 for( std::size_t i = 0; i < size; ++i )
106 denseMatrix[ i ] = value;
110 template<
class DenseMatrix,
class RHS >
111 class DenseMatrixAssignerImplementation< DenseMatrix, RHS, false >
113 template<
class M,
class T>
114 struct have_istl_assign_to_fmatrix
116 struct yes {
char dummy[ 1 ]; };
117 struct no {
char dummy[ 2 ]; };
128 static const bool v =
sizeof( test< const T >( 0 ) ) ==
sizeof( yes );
131 template< class M, class T, bool = have_istl_assign_to_fmatrix< M, T >::v >
132 struct DefaultImplementation;
135 template<
class M,
class T >
136 struct DefaultImplementation< M, T, true >
138 static void apply ( M &m,
const T &
t )
145 template<
class M,
class T >
146 struct DefaultImplementation< M, T, false >
148 static void apply ( M &m,
const T &
t )
150 dune_static_assert( (Conversion< const T, const M >::exists),
"No template specialization of DenseMatrixAssigner found" );
151 m =
static_cast< const M &
>(
t );
156 static void apply ( DenseMatrix &denseMatrix,
const RHS &rhs )
158 DefaultImplementation< DenseMatrix, RHS >::apply( denseMatrix, rhs );
165 template<
class DenseMatrix,
class RHS >
166 struct DenseMatrixAssigner
168 static void apply ( DenseMatrix &denseMatrix,
const RHS &rhs )
170 DenseMatrixAssignerImplementation< DenseMatrix, RHS >::apply( denseMatrix, rhs );
173 #endif // #ifndef DOXYGEN
190 template<
typename MAT>
196 MAT & asImp() {
return static_cast<MAT&
>(*this); }
197 const MAT & asImp()
const {
return static_cast<const MAT&
>(*this); }
237 return asImp().mat_access(i);
242 return asImp().mat_access(i);
303 ConstIterator
end ()
const
324 template<
class RHS >
334 template <
class Other>
337 for (size_type i=0; i<
rows(); i++)
343 template <
class Other>
346 for (size_type i=0; i<
rows(); i++)
354 for (size_type i=0; i<
rows(); i++)
362 for (size_type i=0; i<
rows(); i++)
368 template <
class Other>
371 for( size_type i = 0; i <
rows(); ++i )
372 (*
this)[ i ].
axpy( k, y[ i ] );
377 template <
class Other>
380 for (size_type i=0; i<
rows(); i++)
381 if ((*
this)[i]!=y[i])
386 template <
class Other>
396 template<
class X,
class Y>
397 void mv (
const X& x, Y& y)
const
399 #ifdef DUNE_FMatrix_WITH_CHECKING
403 for (size_type i=0; i<
rows(); ++i)
406 for (size_type j=0; j<
cols(); j++)
407 y[i] += (*
this)[i][j] * x[j];
412 template<
class X,
class Y >
413 void mtv (
const X &x, Y &y )
const
415 #ifdef DUNE_FMatrix_WITH_CHECKING
425 for( size_type i = 0; i <
cols(); ++i )
428 for( size_type j = 0; j <
rows(); ++j )
429 y[ i ] += (*
this)[ j ][ i ] * x[ j ];
434 template<
class X,
class Y>
435 void umv (
const X& x, Y& y)
const
437 #ifdef DUNE_FMatrix_WITH_CHECKING
439 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);
441 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);
443 for (size_type i=0; i<
rows(); i++)
444 for (size_type j=0; j<
cols(); j++)
445 y[i] += (*
this)[i][j] * x[j];
449 template<
class X,
class Y>
450 void umtv (
const X& x, Y& y)
const
452 #ifdef DUNE_FMatrix_WITH_CHECKING
457 for (size_type i=0; i<
rows(); i++)
458 for (size_type j=0; j<
cols(); j++)
459 y[j] += (*
this)[i][j]*x[i];
463 template<
class X,
class Y>
464 void umhv (
const X& x, Y& y)
const
466 #ifdef DUNE_FMatrix_WITH_CHECKING
471 for (size_type i=0; i<
rows(); i++)
472 for (size_type j=0; j<
cols(); j++)
477 template<
class X,
class Y>
478 void mmv (
const X& x, Y& y)
const
480 #ifdef DUNE_FMatrix_WITH_CHECKING
484 for (size_type i=0; i<
rows(); i++)
485 for (size_type j=0; j<
cols(); j++)
486 y[i] -= (*
this)[i][j] * x[j];
490 template<
class X,
class Y>
491 void mmtv (
const X& x, Y& y)
const
493 #ifdef DUNE_FMatrix_WITH_CHECKING
498 for (size_type i=0; i<
rows(); i++)
499 for (size_type j=0; j<
cols(); j++)
500 y[j] -= (*
this)[i][j]*x[i];
504 template<
class X,
class Y>
505 void mmhv (
const X& x, Y& y)
const
507 #ifdef DUNE_FMatrix_WITH_CHECKING
512 for (size_type i=0; i<
rows(); i++)
513 for (size_type j=0; j<
cols(); j++)
518 template<
class X,
class Y>
519 void usmv (
const field_type& alpha,
const X& x, Y& y)
const
521 #ifdef DUNE_FMatrix_WITH_CHECKING
525 for (size_type i=0; i<
rows(); i++)
526 for (size_type j=0; j<
cols(); j++)
527 y[i] += alpha * (*
this)[i][j] * x[j];
531 template<
class X,
class Y>
532 void usmtv (
const field_type& alpha,
const X& x, Y& y)
const
534 #ifdef DUNE_FMatrix_WITH_CHECKING
539 for (size_type i=0; i<
rows(); i++)
540 for (size_type j=0; j<
cols(); j++)
541 y[j] += alpha*(*
this)[i][j]*x[i];
545 template<
class X,
class Y>
546 void usmhv (
const field_type& alpha,
const X& x, Y& y)
const
548 #ifdef DUNE_FMatrix_WITH_CHECKING
553 for (size_type i=0; i<
rows(); i++)
554 for (size_type j=0; j<
cols(); j++)
564 for (size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
565 return fvmeta::sqrt(sum);
572 for (size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
582 ConstIterator it =
begin();
584 for (it = it + 1; it !=
end(); ++it)
585 max = std::max(max, it->one_norm());
596 ConstIterator it =
begin();
598 for (it = it + 1; it !=
end(); ++it)
599 max = std::max(max, it->one_norm_real());
611 void solve (V& x,
const V& b)
const;
623 template<
typename M2>
629 for (size_type i=0; i<
rows(); i++)
630 for (size_type j=0; j<
cols(); j++) {
632 for (size_type k=0; k<
rows(); k++)
633 (*
this)[i][j] += M[i][k]*C[k][j];
640 template<
typename M2>
646 for (size_type i=0; i<
rows(); i++)
647 for (size_type j=0; j<
cols(); j++) {
649 for (size_type k=0; k<
cols(); k++)
650 (*
this)[i][j] += C[i][k]*M[k][j];
662 for (size_type i=0; i<l; i++) {
663 for (size_type j=0; j<
cols(); j++) {
665 for (size_type k=0; k<
rows(); k++)
666 C[i][j] += M[i][k]*(*
this)[k][j];
674 FieldMatrix<K,rows,l> rightmultiplyany (
const FieldMatrix<K,cols,l>& M)
const
676 FieldMatrix<K,rows,l> C;
678 for (size_type i=0; i<
rows(); i++) {
679 for (size_type j=0; j<l; j++) {
681 for (size_type k=0; k<
cols(); k++)
682 C[i][j] += (*
this)[i][k]*M[k][j];
706 return asImp().mat_rows();
712 return asImp().mat_cols();
718 bool exists (size_type i, size_type j)
const
720 #ifdef DUNE_FMatrix_WITH_CHECKING
735 ElimPivot(std::vector<size_type> & pivot);
737 void swap(
int i,
int j);
740 void operator()(
const T&,
int,
int)
743 std::vector<size_type> & pivot_;
751 void swap(
int i,
int j);
753 void operator()(
const typename V::field_type& factor,
int k,
int i);
760 ElimDet(field_type&
sign) : sign_(sign)
766 void operator()(
const field_type&,
int,
int)
774 void luDecomposition(DenseMatrix<MAT>& A, Func func)
const;
778 template<
typename MAT>
779 DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<size_type> & pivot)
782 typedef typename std::vector<size_type>::size_type size_type;
783 for(size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
786 template<
typename MAT>
787 void DenseMatrix<MAT>::ElimPivot::swap(
int i,
int j)
792 template<
typename MAT>
794 DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
798 template<
typename MAT>
800 void DenseMatrix<MAT>::Elim<V>::swap(
int i,
int j)
802 std::swap((*rhs_)[i], (*rhs_)[j]);
805 template<
typename MAT>
807 void DenseMatrix<MAT>::
808 Elim<V>::operator()(
const typename V::field_type& factor,
int k,
int i)
810 (*rhs_)[k] -= factor*(*rhs_)[i];
812 template<
typename MAT>
813 template<
typename Func>
814 inline void DenseMatrix<MAT>::luDecomposition(DenseMatrix<MAT>& A, Func func)
const
816 typedef typename FieldTraits<value_type>::real_type
818 real_type norm = A.infinity_norm_real();
819 real_type pivthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::pivoting_limit() );
820 real_type singthres = std::max( FMatrixPrecision< real_type >::absolute_limit(), norm * FMatrixPrecision< real_type >::singular_limit() );
823 for (size_type i=0; i<rows(); i++)
825 typename FieldTraits<value_type>::real_type pivmax=fvmeta::absreal(A[i][i]);
832 typename FieldTraits<value_type>::real_type abs(0.0);
833 for (size_type k=i+1; k<rows(); k++)
834 if ((abs=fvmeta::absreal(A[k][i]))>pivmax)
836 pivmax = abs; imax = k;
840 for (size_type j=0; j<rows(); j++)
841 std::swap(A[i][j],A[imax][j]);
847 if (pivmax<singthres)
848 DUNE_THROW(FMatrixError,
"matrix is singular");
851 for (size_type k=i+1; k<rows(); k++)
853 field_type factor = A[k][i]/A[i][i];
855 for (size_type j=i+1; j<rows(); j++)
856 A[k][j] -= factor*A[i][j];
862 template<
typename MAT>
864 inline void DenseMatrix<MAT>::solve(V& x,
const V& b)
const
868 DUNE_THROW(FMatrixError,
"Can't solve for a " << rows() <<
"x" << cols() <<
" matrix!");
872 #ifdef DUNE_FMatrix_WITH_CHECKING
873 if (fvmeta::absreal((*
this)[0][0])<FMatrixPrecision<>::absolute_limit())
874 DUNE_THROW(FMatrixError,
"matrix is singular");
876 x[0] = b[0]/(*this)[0][0];
879 else if (rows()==2) {
881 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
882 #ifdef DUNE_FMatrix_WITH_CHECKING
883 if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
884 DUNE_THROW(FMatrixError,
"matrix is singular");
888 x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
889 x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
892 else if (rows()==3) {
894 field_type d = determinant();
895 #ifdef DUNE_FMatrix_WITH_CHECKING
896 if (fvmeta::absreal(d)<FMatrixPrecision<>::absolute_limit())
897 DUNE_THROW(FMatrixError,
"matrix is singular");
900 x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
901 - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
902 + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
904 x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
905 - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
906 + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
908 x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
909 - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
910 + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
920 luDecomposition(A, elim);
923 for(
int i=rows()-1; i>=0; i--) {
924 for (size_type j=i+1; j<rows(); j++)
925 rhs[i] -= A[i][j]*x[j];
926 x[i] = rhs[i]/A[i][i];
931 template<
typename MAT>
932 inline void DenseMatrix<MAT>::invert()
936 DUNE_THROW(FMatrixError,
"Can't invert a " << rows() <<
"x" << cols() <<
" matrix!");
940 #ifdef DUNE_FMatrix_WITH_CHECKING
941 if (fvmeta::absreal((*
this)[0][0])<FMatrixPrecision<>::absolute_limit())
942 DUNE_THROW(FMatrixError,
"matrix is singular");
944 (*this)[0][0] = field_type( 1 ) / (*this)[0][0];
947 else if (rows()==2) {
949 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
950 #ifdef DUNE_FMatrix_WITH_CHECKING
951 if (fvmeta::absreal(detinv)<FMatrixPrecision<>::absolute_limit())
952 DUNE_THROW(FMatrixError,
"matrix is singular");
954 detinv = field_type( 1 ) / detinv;
956 field_type temp=(*this)[0][0];
957 (*this)[0][0] = (*this)[1][1]*detinv;
958 (*this)[0][1] = -(*this)[0][1]*detinv;
959 (*this)[1][0] = -(*this)[1][0]*detinv;
960 (*this)[1][1] = temp*detinv;
966 std::vector<size_type> pivot(rows());
967 luDecomposition(A, ElimPivot(pivot));
968 DenseMatrix<MAT>& L=A;
969 DenseMatrix<MAT>& U=A;
974 for(size_type i=0; i<rows(); ++i)
978 for (size_type i=0; i<rows(); i++)
979 for (size_type j=0; j<i; j++)
980 for (size_type k=0; k<rows(); k++)
981 (*
this)[i][k] -= L[i][j]*(*this)[j][k];
984 for (size_type i=rows(); i>0;) {
986 for (size_type k=0; k<rows(); k++) {
987 for (size_type j=i+1; j<rows(); j++)
988 (*
this)[i][k] -= U[i][j]*(*this)[j][k];
989 (*this)[i][k] /= U[i][i];
993 for(size_type i=rows(); i>0; ) {
996 for(size_type j=0; j<rows(); ++j)
997 std::swap((*
this)[j][pivot[i]], (*this)[j][i]);
1003 template<
typename MAT>
1004 inline typename DenseMatrix<MAT>::field_type
1005 DenseMatrix<MAT>::determinant()
const
1009 DUNE_THROW(FMatrixError,
"There is no determinant for a " << rows() <<
"x" << cols() <<
" matrix!");
1012 return (*
this)[0][0];
1015 return (*
this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
1019 field_type t4 = (*this)[0][0] * (*this)[1][1];
1020 field_type t6 = (*this)[0][0] * (*this)[1][2];
1021 field_type t8 = (*this)[0][1] * (*this)[1][0];
1022 field_type t10 = (*this)[0][2] * (*this)[1][0];
1023 field_type t12 = (*this)[0][1] * (*this)[2][0];
1024 field_type t14 = (*this)[0][2] * (*this)[2][0];
1026 return (t4*(*
this)[2][2]-t6*(*
this)[2][1]-t8*(*
this)[2][2]+
1027 t10*(*
this)[2][1]+t12*(*
this)[1][2]-t14*(*
this)[1][1]);
1035 luDecomposition(A, ElimDet(det));
1037 catch (FMatrixError&)
1041 for (size_type i = 0; i < rows(); ++i)
1048 namespace DenseMatrixHelp {
1050 template <
typename K>
1054 inverse[0][0] = 1.0/matrix[0][0];
1055 return matrix[0][0];
1059 template <
typename K>
1067 template <
typename K>
1071 field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
1072 field_type det_1 = 1.0/det;
1073 inverse[0][0] = matrix[1][1] * det_1;
1074 inverse[0][1] = - matrix[0][1] * det_1;
1075 inverse[1][0] = - matrix[1][0] * det_1;
1076 inverse[1][1] = matrix[0][0] * det_1;
1082 template <
typename K>
1086 field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
1087 field_type det_1 = 1.0/det;
1088 inverse[0][0] = matrix[1][1] * det_1;
1089 inverse[1][0] = - matrix[0][1] * det_1;
1090 inverse[0][1] = - matrix[1][0] * det_1;
1091 inverse[1][1] = matrix[0][0] * det_1;
1096 template <
typename K>
1100 field_type t4 = matrix[0][0] * matrix[1][1];
1101 field_type t6 = matrix[0][0] * matrix[1][2];
1102 field_type t8 = matrix[0][1] * matrix[1][0];
1103 field_type t10 = matrix[0][2] * matrix[1][0];
1104 field_type t12 = matrix[0][1] * matrix[2][0];
1105 field_type t14 = matrix[0][2] * matrix[2][0];
1107 field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
1108 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
1109 field_type t17 = 1.0/det;
1111 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
1112 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
1113 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
1114 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
1115 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
1116 inverse[1][2] = -(t6-t10) * t17;
1117 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
1118 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
1119 inverse[2][2] = (t4-t8) * t17;
1125 template <
typename K>
1129 field_type t4 = matrix[0][0] * matrix[1][1];
1130 field_type t6 = matrix[0][0] * matrix[1][2];
1131 field_type t8 = matrix[0][1] * matrix[1][0];
1132 field_type t10 = matrix[0][2] * matrix[1][0];
1133 field_type t12 = matrix[0][1] * matrix[2][0];
1134 field_type t14 = matrix[0][2] * matrix[2][0];
1136 field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
1137 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
1138 field_type t17 = 1.0/det;
1140 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
1141 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
1142 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
1143 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
1144 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
1145 inverse[2][1] = -(t6-t10) * t17;
1146 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
1147 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
1148 inverse[2][2] = (t4-t8) * t17;
1154 template<
class K,
int m,
int n,
int p >
1161 for( size_type i = 0; i < m; ++i )
1163 for( size_type j = 0; j < p; ++j )
1165 ret[ i ][ j ] = K( 0 );
1166 for( size_type k = 0; k < n; ++k )
1167 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
1173 template <
typename K,
int rows,
int cols>
1178 for(size_type i=0; i<cols(); i++)
1179 for(size_type j=0; j<cols(); j++)
1182 for(size_type k=0; k<rows(); k++)
1183 ret[i][j]+=matrix[k][i]*matrix[k][j];
1189 template <
typename MAT,
typename V1,
typename V2>
1193 assert(ret.
size() == matrix.
rows());
1196 for(size_type i=0; i<matrix.
rows(); ++i)
1199 for(size_type j=0; j<matrix.
cols(); ++j)
1201 ret[i] += matrix[i][j]*x[j];
1207 template <
typename K,
int rows,
int cols>
1213 for(size_type i=0; i<cols(); ++i)
1216 for(size_type j=0; j<rows(); ++j)
1217 ret[i] += matrix[j][i]*x[j];
1222 template <
typename K,
int rows,
int cols>
1223 static inline FieldVector<K,rows>
mult(
const FieldMatrix<K,rows,cols> &matrix,
const FieldVector<K,cols> & x)
1225 FieldVector<K,rows> ret;
1231 template <
typename K,
int rows,
int cols>
1232 static inline FieldVector<K,cols>
multTransposed(
const FieldMatrix<K,rows,cols> &matrix,
const FieldVector<K,rows> & x)
1234 FieldVector<K,cols> ret;
1243 template<
typename MAT>
1244 std::ostream& operator<< (std::ostream& s, const DenseMatrix<MAT>& a)
1247 s << a[i] << std::endl;
ConstIterator begin() const
begin iterator
Definition: densematrix.hh:297
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:178
DenseMatrix & axpy(const field_type &k, const DenseMatrix< Other > &y)
vector space axpy operation (*this += k y)
Definition: densematrix.hh:369
Traits::value_type field_type
export the type representing the field
Definition: densematrix.hh:209
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:120
DenseIterator< const DenseMatrix, const row_type > ConstIterator
Iterator class for sequential access.
Definition: densematrix.hh:288
T t
Definition: alignment.hh:38
static FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition: fmatrix.hh:504
Iterator RowIterator
rename the iterators for easier access
Definition: densematrix.hh:257
Dune namespace.
Definition: alignment.hh:13
A few common exception classes.
Fallback implementation of the C++0x static_assert feature.
const FieldTraits< typename DenseMatVecTraits< M >::value_type >::real_type real_type
Definition: densematrix.hh:30
void umhv(const X &x, Y &y) const
y += A^H x
Definition: densematrix.hh:464
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: densematrix.hh:569
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:341
Interface for a class of dense vectors over a given field.
Definition: densevector.hh:17
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:215
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:641
Traits::derived_type derived_type
type of derived matrix class
Definition: densematrix.hh:203
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:103
DenseIterator< DenseMatrix, row_type > Iterator
Iterator class for sequential access.
Definition: densematrix.hh:253
FieldVector()
Constructor making default-initialized vector.
Definition: fvector.hh:106
A free function to provide the demangled class name of a given object or type as a string...
Iterator beforeEnd()
Definition: densematrix.hh:275
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:218
size_type N() const
number of rows
Definition: densematrix.hh:692
Traits::value_type value_type
export the type representing the field
Definition: densematrix.hh:206
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: densematrix.hh:718
vector space out of a tensor product of fields.
Definition: densematrix.hh:38
you have to specialize this structure for any type that should be assignable to a DenseMatrix ...
Definition: densematrix.hh:73
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition: densematrix.hh:624
ConstIterator beforeEnd() const
Definition: densematrix.hh:310
static K invertMatrix(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition: fmatrix.hh:333
FieldTraits< value_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: densematrix.hh:577
row_type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: densematrix.hh:294
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition: densematrix.hh:1190
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:436
DenseMatrix & operator-=(const DenseMatrix< Other > &y)
vector space subtraction
Definition: densematrix.hh:344
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition: densevector.hh:569
#define dune_static_assert(COND, MSG)
Helper template so that compilation fails if condition is not true.
Definition: static_assert.hh:79
static FieldVector< K, cols > multTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x)
calculates ret = matrix^T * x
Definition: fmatrix.hh:513
bool operator!=(const DenseMatrix< Other > &y) const
Binary matrix incomparison.
Definition: densematrix.hh:387
void invert()
Compute inverse.
void mmv(const X &x, Y &y) const
y -= A x
Definition: densematrix.hh:478
void mv(const X &x, Y &y) const
y = A x
Definition: densematrix.hh:397
DenseMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:360
#define DUNE_THROW(E, m)
Definition: exceptions.hh:244
Removes a const qualifier while preserving others.
Definition: typetraits.hh:175
DenseMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:352
ConstIterator const_iterator
typedef for stl compliant access
Definition: densematrix.hh:290
size_type cols() const
number of columns
Definition: densematrix.hh:710
bool operator==(const DenseMatrix< Other > &y) const
Binary matrix comparison.
Definition: densematrix.hh:378
size_type M() const
number of columns
Definition: densematrix.hh:698
row_reference operator[](size_type i)
random access
Definition: densematrix.hh:235
size_type size() const
size method
Definition: densevector.hh:285
Traits::value_type block_type
export the type representing the components
Definition: densematrix.hh:212
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:119
void solve(V &x, const V &b) const
Solve system A x = b.
Implements a vector constructed from a given type representing a field and a compile-time given size...
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:490
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:532
Iterator iterator
typedef for stl compliant access
Definition: densematrix.hh:255
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
T real_type
export the type representing the real type of the field
Definition: ftraits.hh:28
A dense n x m matrix.
Definition: densematrix.hh:24
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: densematrix.hh:292
void mtv(const X &x, Y &y) const
y = A^T x
Definition: densematrix.hh:413
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18
static void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition: fmatrix.hh:455
FieldTraits< value_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densematrix.hh:591
Definition: ftraits.hh:23
DenseMatrix & operator+=(const DenseMatrix< Other > &y)
vector space addition
Definition: densematrix.hh:335
void usmv(const field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition: densematrix.hh:519
DUNE_CONSTEXPR size_type size() const
Definition: fvector.hh:163
ConstIterator end() const
end iterator
Definition: densematrix.hh:303
The number of block levels we contain. This is 1.
Definition: densematrix.hh:229
ConstIterator beforeBegin() const
Definition: densematrix.hh:317
void umv(const X &x, Y &y) const
y += A x
Definition: densematrix.hh:435
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: densematrix.hh:505
Base::size_type size_type
Definition: fmatrix.hh:80
Iterator begin()
begin iterator
Definition: densematrix.hh:262
Various precision settings for calculations with FieldMatrix and FieldVector.
size_type size() const
size method (number of rows)
Definition: densematrix.hh:246
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: densematrix.hh:546
size_type rows() const
number of rows
Definition: densematrix.hh:704
const FieldTraits< typename DenseMatVecTraits< M >::value_type >::field_type field_type
Definition: densematrix.hh:29
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:224
DenseMatrix & operator=(const RHS &rhs)
Definition: densematrix.hh:325
row_type::Iterator ColIterator
rename the iterators for easier access
Definition: densematrix.hh:259
Base class for Dune-Exceptions.
Definition: exceptions.hh:92
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: densematrix.hh:491
void umtv(const X &x, Y &y) const
y += A^T x
Definition: densematrix.hh:450
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: densematrix.hh:561
Iterator end()
end iterator
Definition: densematrix.hh:268
void istl_assign_to_fmatrix(DenseMatrix &denseMatrix, const K(&values)[M][N])
Definition: densematrix.hh:78
Iterator beforeBegin()
Definition: densematrix.hh:282
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:221
Definition: matvectraits.hh:30
Some useful basic math stuff.
A dense n x m matrix.
Definition: densematrix.hh:37