3#ifndef DUNE_DENSEMATRIX_HH
4#define DUNE_DENSEMATRIX_HH
22 template<
typename M>
class DenseMatrix;
25 struct FieldTraits< DenseMatrix<M> >
27 typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::field_type field_type;
28 typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::real_type real_type;
35 template<
class K,
int N,
int M>
class FieldMatrix;
36 template<
class K,
int N>
class FieldVector;
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 >
76 void istl_assign_to_fmatrix (
DenseMatrix &denseMatrix,
const K (&values)[ M ][ N ] )
78 for(
int i = 0; i < N; ++i )
79 for(
int j = 0; j < M; ++j )
80 denseMatrix[ i ][ j ] = values[ i ][ j ];
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 ]; };
116 static yes test(
decltype( istl_assign_to_fmatrix( get_ref< M >(), get_ref< C >() ) ) * );
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 )
133 istl_assign_to_fmatrix( m, t );
138 template<
class M,
class T >
139 struct DefaultImplementation< M, T, false >
141 static void apply ( M &m,
const T &t )
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 );
183 template<
typename MAT>
186 typedef DenseMatVecTraits<MAT> Traits;
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;
317 template<
class RHS >
327 template <
class Other>
336 template <
class Other>
361 template <
class Other>
365 (*
this)[ i ].
axpy( k, y[ i ] );
370 template <
class Other>
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) );
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) );
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);
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
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
468 template<
class X,
class Y>
469 void mmv (
const X& x, Y& y)
const
471#ifdef DUNE_FMatrix_WITH_CHECKING
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
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
509 template<
class X,
class Y>
512#ifdef DUNE_FMatrix_WITH_CHECKING
518 y[i] += alpha * (*
this)[i][j] * x[j];
522 template<
class X,
class Y>
525#ifdef DUNE_FMatrix_WITH_CHECKING
532 y[j] += alpha*(*
this)[i][j]*x[i];
536 template<
class X,
class Y>
539#ifdef DUNE_FMatrix_WITH_CHECKING
554 typename FieldTraits<value_type>::real_type sum=(0.0);
555 for (
size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
556 return fvmeta::sqrt(sum);
562 typename FieldTraits<value_type>::real_type sum=(0.0);
563 for (
size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
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());
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>
617 assert(
M.rows() ==
M.cols() &&
M.rows() ==
rows());
624 (*
this)[i][j] +=
M[i][k]*C[k][j];
631 template<
typename M2>
634 assert(
M.rows() ==
M.cols() &&
M.cols() ==
cols());
641 (*
this)[i][j] += C[i][k]*
M[k][j];
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;
673 C[i][j] += (*
this)[i][k]*
M[k][j];
697 return asImp().mat_rows();
703 return asImp().mat_cols();
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);
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
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>
859 DUNE_THROW(FMatrixError,
"Can't solve for a " << rows() <<
"x" << cols() <<
" matrix!");
863#ifdef DUNE_FMatrix_WITH_CHECKING
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
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
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>
927 DUNE_THROW(FMatrixError,
"Can't invert a " << rows() <<
"x" << cols() <<
" matrix!");
931#ifdef DUNE_FMatrix_WITH_CHECKING
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
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>
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 {
1042 template <
typename K>
1043 static inline K invertMatrix (
const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
1045 inverse[0][0] = 1.0/matrix[0][0];
1046 return matrix[0][0];
1050 template <
typename K>
1051 static inline K invertMatrix_retTransposed (
const FieldMatrix<K,1,1> &matrix, FieldMatrix<K,1,1> &inverse)
1053 return invertMatrix(matrix,inverse);
1058 template <
typename K>
1059 static inline K invertMatrix (
const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
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>
1074 static inline K invertMatrix_retTransposed (
const FieldMatrix<K,2,2> &matrix, FieldMatrix<K,2,2> &inverse)
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>
1088 static inline K invertMatrix (
const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
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>
1117 static inline K invertMatrix_retTransposed (
const FieldMatrix<K,3,3> &matrix, FieldMatrix<K,3,3> &inverse)
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 >
1146 static inline void multMatrix (
const FieldMatrix< K, m, n > &A,
1147 const FieldMatrix< K, n, p > &B,
1148 FieldMatrix< K, m, p > &ret )
1150 typedef typename FieldMatrix< K, m, p > :: size_type size_type;
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>
1165 static inline void multTransposedMatrix(
const FieldMatrix<K,rows,cols> &matrix, FieldMatrix<K,cols,cols>& ret)
1167 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
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];
1199 template <
typename K,
int rows,
int cols>
1202 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
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;
1226 multAssignTransposed( matrix, x, ret );
1234 template<
typename MAT>
1238 s << a[i] << std::endl;
Checks wether a type is convertible to another.
Definition: typetraits.hh:177
@ exists
True if the conversion exists.
Definition: typetraits.hh:187
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:127
A dense n x m matrix.
Definition: densematrix.hh:185
DenseMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:353
ConstIterator const_iterator
typedef for stl compliant access
Definition: densematrix.hh:283
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: densematrix.hh:523
bool operator!=(const DenseMatrix< Other > &y) const
Binary matrix incomparison.
Definition: densematrix.hh:380
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
Traits::value_type field_type
export the type representing the field
Definition: densematrix.hh:202
void solve(V &x, const V &b) const
Solve system A x = b.
size_type cols() const
number of columns
Definition: densematrix.hh:701
ConstIterator beforeEnd() const
Definition: densematrix.hh:303
bool operator==(const DenseMatrix< Other > &y) const
Binary matrix comparison.
Definition: densematrix.hh:371
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: densematrix.hh:482
DenseMatrix & operator+=(const DenseMatrix< Other > &y)
vector space addition
Definition: densematrix.hh:328
ConstIterator beforeBegin() const
Definition: densematrix.hh:310
Traits::value_type block_type
export the type representing the components
Definition: densematrix.hh:205
void mtv(const X &x, Y &y) const
y = A^T x
Definition: densematrix.hh:407
DenseMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:345
remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: densematrix.hh:287
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: densematrix.hh:537
DenseMatrix & operator-=(const DenseMatrix< Other > &y)
vector space subtraction
Definition: densematrix.hh:337
FieldTraits< value_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densematrix.hh:582
size_type size() const
size method (number of rows)
Definition: densematrix.hh:239
size_type M() const
number of columns
Definition: densematrix.hh:689
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:632
remove_reference< row_reference >::type::Iterator ColIterator
rename the iterators for easier access
Definition: densematrix.hh:252
void usmv(const field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition: densematrix.hh:510
Iterator end()
end iterator
Definition: densematrix.hh:261
Iterator beforeBegin()
Definition: densematrix.hh:275
Iterator beforeEnd()
Definition: densematrix.hh:268
Traits::value_type value_type
export the type representing the field
Definition: densematrix.hh:199
@ blocklevel
The number of block levels we contain. This is 1.
Definition: densematrix.hh:222
size_type rows() const
number of rows
Definition: densematrix.hh:695
void mmv(const X &x, Y &y) const
y -= A x
Definition: densematrix.hh:469
void invert()
Compute inverse.
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition: densematrix.hh:615
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: densematrix.hh:496
Traits::derived_type derived_type
type of derived matrix class
Definition: densematrix.hh:196
row_reference operator[](size_type i)
random access
Definition: densematrix.hh:228
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: densematrix.hh:709
DenseMatrix & axpy(const field_type &k, const DenseMatrix< Other > &y)
vector space axpy operation (*this += k y)
Definition: densematrix.hh:362
Iterator RowIterator
rename the iterators for easier access
Definition: densematrix.hh:250
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: densematrix.hh:552
void umv(const X &x, Y &y) const
y += A x
Definition: densematrix.hh:426
DenseIterator< const DenseMatrix, const row_type, const_row_reference > ConstIterator
Iterator class for sequential access.
Definition: densematrix.hh:281
size_type N() const
number of rows
Definition: densematrix.hh:683
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:211
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:208
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
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: densematrix.hh:560
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:214
Iterator iterator
typedef for stl compliant access
Definition: densematrix.hh:248
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: densematrix.hh:285
field_type determinant() const
calculates the determinant of this matrix
DenseIterator< DenseMatrix, row_type, row_reference > Iterator
Iterator class for sequential access.
Definition: densematrix.hh:246
void umtv(const X &x, Y &y) const
y += A^T x
Definition: densematrix.hh:441
ConstIterator begin() const
begin iterator
Definition: densematrix.hh:290
Iterator begin()
begin iterator
Definition: densematrix.hh:255
void umhv(const X &x, Y &y) const
y += A^H x
Definition: densematrix.hh:455
ConstIterator end() const
end iterator
Definition: densematrix.hh:296
Interface for a class of dense vectors over a given field.
Definition: densevector.hh:234
size_type size() const
size method
Definition: densevector.hh:296
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition: densevector.hh:643
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:171
static ctype pivoting_limit()
return threshold to do pivoting
Definition: precision.hh:26
static ctype absolute_limit()
return threshold to declare matrix singular
Definition: precision.hh:50
static ctype singular_limit()
return threshold to declare matrix singular
Definition: precision.hh:38
A dense n x m matrix.
Definition: fmatrix.hh:67
vector space out of a tensor product of fields.
Definition: fvector.hh:94
FieldVector()
Constructor making default-initialized vector.
Definition: fvector.hh:108
Default exception class for mathematical errors.
Definition: exceptions.hh:266
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: matrix.hh:518
A free function to provide the demangled class name of a given object or type as a string.
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition: densematrix.hh:1181
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
std::ostream & operator<<(std::ostream &s, const array< T, N > &e)
Output operator for array.
Definition: array.hh:26
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
Some useful basic math stuff.
Dune namespace.
Definition: alignment.hh:10
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:119
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:103
Various precision settings for calculations with FieldMatrix and FieldVector.
you have to specialize this structure for any type that should be assignable to a DenseMatrix
Definition: densematrix.hh:71
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentional unused function parameters with.
Definition: unused.hh:18