3#ifndef DUNE_DENSEMATRIX_HH
4#define DUNE_DENSEMATRIX_HH
24 template<
typename M>
class DenseMatrix;
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 >
79 bool primitive = std::is_convertible< RHS, typename DenseMatrix::field_type >::value >
80 class DenseMatrixAssignerImplementation;
82 template<
class DenseMatrix,
class RHS >
83 class DenseMatrixAssignerImplementation<
DenseMatrix, RHS, true >
86 static void apply (
DenseMatrix &denseMatrix,
const RHS &rhs )
89 std::fill( denseMatrix.
begin(), denseMatrix.
end(),
static_cast< field_type
>( rhs ) );
93 template<
class DenseMatrix,
class RHS >
94 class DenseMatrixAssignerImplementation< DenseMatrix, RHS, false >
97 static void apply ( DenseMatrix &denseMatrix,
const RHS &rhs )
99 static_assert( (std::is_convertible< const RHS, const DenseMatrix >::value),
100 "No template specialization of DenseMatrixAssigner found" );
101 denseMatrix =
static_cast< const DenseMatrix &
>( rhs );
108 template<
class DenseMatrix,
class RHS >
109 struct DenseMatrixAssigner
111 static void apply ( DenseMatrix &denseMatrix,
const RHS &rhs )
113 DenseMatrixAssignerImplementation< DenseMatrix, RHS >::apply( denseMatrix, rhs );
133 template<
typename MAT>
139 MAT & asImp() {
return static_cast<MAT&
>(*this); }
140 const MAT & asImp()
const {
return static_cast<const MAT&
>(*this); }
180 return asImp().mat_access(i);
185 return asImp().mat_access(i);
202 typedef typename std::remove_reference<row_reference>::type::Iterator
ColIterator;
237 typedef typename std::remove_reference<const_row_reference>::type::ConstIterator
ConstColIterator;
267 template<
class RHS >
277 template <
class Other>
287 template <
class Other>
313 template <
class Other>
318 (*
this)[ i ].axpy( k, y[ i ] );
323 template <
class Other>
328 if ((*
this)[i]!=y[i])
333 template <
class Other>
343 template<
class X,
class Y>
344 void mv (
const X& x, Y& y)
const
355 y[i] += (*
this)[i][j] * x[j];
360 template<
class X,
class Y >
361 void mtv (
const X &x, Y &y )
const
372 y[i] += (*
this)[j][i] * x[j];
377 template<
class X,
class Y>
378 void umv (
const X& x, Y& y)
const
384 y[i] += (*
this)[i][j] * x[j];
388 template<
class X,
class Y>
389 void umtv (
const X& x, Y& y)
const
395 y[j] += (*
this)[i][j]*x[i];
399 template<
class X,
class Y>
400 void umhv (
const X& x, Y& y)
const
410 template<
class X,
class Y>
411 void mmv (
const X& x, Y& y)
const
417 y[i] -= (*
this)[i][j] * x[j];
421 template<
class X,
class Y>
422 void mmtv (
const X& x, Y& y)
const
428 y[j] -= (*
this)[i][j]*x[i];
432 template<
class X,
class Y>
433 void mmhv (
const X& x, Y& y)
const
443 template<
class X,
class Y>
445 const X& x, Y& y)
const
451 y[i] += alpha * (*
this)[i][j] * x[j];
455 template<
class X,
class Y>
457 const X& x, Y& y)
const
463 y[j] += alpha*(*
this)[i][j]*x[i];
467 template<
class X,
class Y>
469 const X& x, Y& y)
const
484 for (
size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
485 return fvmeta::sqrt(sum);
492 for (
size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
498 typename std::enable_if<!has_nan<vt>::value,
int>::type = 0>
504 for (
auto const &x : *
this) {
505 real_type
const a = x.one_norm();
513 typename std::enable_if<!has_nan<vt>::value,
int>::type = 0>
519 for (
auto const &x : *
this) {
520 real_type
const a = x.one_norm_real();
528 typename std::enable_if<has_nan<vt>::value,
int>::type = 0>
535 for (
auto const &x : *
this) {
536 real_type
const a = x.one_norm();
546 typename std::enable_if<has_nan<vt>::value,
int>::type = 0>
553 for (
auto const &x : *
this) {
554 real_type
const a = x.one_norm_real();
569 void solve (V& x,
const V& b)
const;
581 template<
typename M2>
592 (*
this)[i][j] +=
M[i][k]*C[k][j];
599 template<
typename M2>
610 (*
this)[i][j] += C[i][k]*
M[k][j];
626 C[i][j] +=
M[i][k]*(*
this)[k][j];
634 FieldMatrix<K,rows,l> rightmultiplyany (
const FieldMatrix<K,cols,l>&
M)
const
636 FieldMatrix<K,rows,l> C;
642 C[i][j] += (*
this)[i][k]*
M[k][j];
666 return asImp().mat_rows();
672 return asImp().mat_cols();
692 ElimPivot(std::vector<size_type> & pivot);
694 void swap(
int i,
int j);
697 void operator()(
const T&,
int,
int)
700 std::vector<size_type> & pivot_;
708 void swap(
int i,
int j);
710 void operator()(
const typename V::field_type& factor,
int k,
int i);
717 ElimDet(field_type& sign) : sign_(
sign)
723 void operator()(
const field_type&,
int,
int)
731 void luDecomposition(DenseMatrix<MAT>& A, Func func)
const;
735 template<
typename MAT>
736 DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<size_type> & pivot)
739 typedef typename std::vector<size_type>::size_type
size_type;
740 for(
size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
743 template<
typename MAT>
744 void DenseMatrix<MAT>::ElimPivot::swap(
int i,
int j)
749 template<
typename MAT>
751 DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
755 template<
typename MAT>
757 void DenseMatrix<MAT>::Elim<V>::swap(
int i,
int j)
759 std::swap((*rhs_)[i], (*rhs_)[j]);
762 template<
typename MAT>
764 void DenseMatrix<MAT>::
765 Elim<V>::operator()(
const typename V::field_type& factor,
int k,
int i)
767 (*rhs_)[k] -= factor*(*rhs_)[i];
769 template<
typename MAT>
770 template<
typename Func>
771 inline void DenseMatrix<MAT>::luDecomposition(DenseMatrix<MAT>& A, Func func)
const
775 real_type norm = A.infinity_norm_real();
780 for (size_type i=0; i<rows(); i++)
790 for (size_type k=i+1; k<rows(); k++)
791 if ((abs=fvmeta::absreal(A[k][i]))>pivmax)
793 pivmax = abs; imax = k;
797 for (size_type j=0; j<rows(); j++)
798 std::swap(A[i][j],A[imax][j]);
804 if (pivmax<singthres)
805 DUNE_THROW(FMatrixError,
"matrix is singular");
808 for (size_type k=i+1; k<rows(); k++)
810 field_type factor = A[k][i]/A[i][i];
812 for (size_type j=i+1; j<rows(); j++)
813 A[k][j] -= factor*A[i][j];
819 template<
typename MAT>
825 DUNE_THROW(FMatrixError,
"Can't solve for a " << rows() <<
"x" << cols() <<
" matrix!");
829#ifdef DUNE_FMatrix_WITH_CHECKING
831 DUNE_THROW(FMatrixError,
"matrix is singular");
833 x[0] = b[0]/(*this)[0][0];
836 else if (rows()==2) {
838 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
839#ifdef DUNE_FMatrix_WITH_CHECKING
841 DUNE_THROW(FMatrixError,
"matrix is singular");
845 x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
846 x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
849 else if (rows()==3) {
851 field_type d = determinant();
852#ifdef DUNE_FMatrix_WITH_CHECKING
854 DUNE_THROW(FMatrixError,
"matrix is singular");
857 x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
858 - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
859 + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
861 x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
862 - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
863 + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
865 x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
866 - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
867 + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
877 luDecomposition(A, elim);
880 for(
int i=rows()-1; i>=0; i--) {
881 for (size_type j=i+1; j<rows(); j++)
882 rhs[i] -= A[i][j]*x[j];
883 x[i] = rhs[i]/A[i][i];
888 template<
typename MAT>
893 DUNE_THROW(FMatrixError,
"Can't invert a " << rows() <<
"x" << cols() <<
" matrix!");
897#ifdef DUNE_FMatrix_WITH_CHECKING
899 DUNE_THROW(FMatrixError,
"matrix is singular");
901 (*this)[0][0] = field_type( 1 ) / (*this)[0][0];
904 else if (rows()==2) {
906 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
907#ifdef DUNE_FMatrix_WITH_CHECKING
909 DUNE_THROW(FMatrixError,
"matrix is singular");
911 detinv = field_type( 1 ) / detinv;
913 field_type temp=(*this)[0][0];
914 (*this)[0][0] = (*this)[1][1]*detinv;
915 (*this)[0][1] = -(*this)[0][1]*detinv;
916 (*this)[1][0] = -(*this)[1][0]*detinv;
917 (*this)[1][1] = temp*detinv;
923 std::vector<size_type> pivot(rows());
924 luDecomposition(A, ElimPivot(pivot));
925 DenseMatrix<MAT>& L=A;
926 DenseMatrix<MAT>& U=A;
931 for(size_type i=0; i<rows(); ++i)
935 for (size_type i=0; i<rows(); i++)
936 for (size_type j=0; j<i; j++)
937 for (size_type k=0; k<rows(); k++)
938 (*
this)[i][k] -= L[i][j]*(*this)[j][k];
941 for (size_type i=rows(); i>0;) {
943 for (size_type k=0; k<rows(); k++) {
944 for (size_type j=i+1; j<rows(); j++)
945 (*
this)[i][k] -= U[i][j]*(*this)[j][k];
946 (*this)[i][k] /= U[i][i];
950 for(size_type i=rows(); i>0; ) {
953 for(size_type j=0; j<rows(); ++j)
954 std::swap((*
this)[j][pivot[i]], (*this)[j][i]);
960 template<
typename MAT>
966 DUNE_THROW(FMatrixError,
"There is no determinant for a " << rows() <<
"x" << cols() <<
" matrix!");
969 return (*
this)[0][0];
972 return (*
this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
976 field_type t4 = (*this)[0][0] * (*this)[1][1];
977 field_type t6 = (*this)[0][0] * (*this)[1][2];
978 field_type t8 = (*this)[0][1] * (*this)[1][0];
979 field_type t10 = (*this)[0][2] * (*this)[1][0];
980 field_type t12 = (*this)[0][1] * (*this)[2][0];
981 field_type t14 = (*this)[0][2] * (*this)[2][0];
983 return (t4*(*
this)[2][2]-t6*(*
this)[2][1]-t8*(*
this)[2][2]+
984 t10*(*
this)[2][1]+t12*(*
this)[1][2]-t14*(*
this)[1][1]);
992 luDecomposition(A, ElimDet(det));
994 catch (FMatrixError&)
998 for (size_type i = 0; i < rows(); ++i)
1005 namespace DenseMatrixHelp {
1008 template <
typename K>
1011 inverse[0][0] = 1.0/matrix[0][0];
1012 return matrix[0][0];
1016 template <
typename K>
1019 return invertMatrix(matrix,inverse);
1024 template <
typename K>
1028 field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
1029 field_type det_1 = 1.0/det;
1030 inverse[0][0] = matrix[1][1] * det_1;
1031 inverse[0][1] = - matrix[0][1] * det_1;
1032 inverse[1][0] = - matrix[1][0] * det_1;
1033 inverse[1][1] = matrix[0][0] * det_1;
1039 template <
typename K>
1043 field_type det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
1044 field_type det_1 = 1.0/det;
1045 inverse[0][0] = matrix[1][1] * det_1;
1046 inverse[1][0] = - matrix[0][1] * det_1;
1047 inverse[0][1] = - matrix[1][0] * det_1;
1048 inverse[1][1] = matrix[0][0] * det_1;
1053 template <
typename K>
1057 field_type t4 = matrix[0][0] * matrix[1][1];
1058 field_type t6 = matrix[0][0] * matrix[1][2];
1059 field_type t8 = matrix[0][1] * matrix[1][0];
1060 field_type t10 = matrix[0][2] * matrix[1][0];
1061 field_type t12 = matrix[0][1] * matrix[2][0];
1062 field_type t14 = matrix[0][2] * matrix[2][0];
1064 field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
1065 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
1066 field_type t17 = 1.0/det;
1068 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
1069 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
1070 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
1071 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
1072 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
1073 inverse[1][2] = -(t6-t10) * t17;
1074 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
1075 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
1076 inverse[2][2] = (t4-t8) * t17;
1082 template <
typename K>
1086 field_type t4 = matrix[0][0] * matrix[1][1];
1087 field_type t6 = matrix[0][0] * matrix[1][2];
1088 field_type t8 = matrix[0][1] * matrix[1][0];
1089 field_type t10 = matrix[0][2] * matrix[1][0];
1090 field_type t12 = matrix[0][1] * matrix[2][0];
1091 field_type t14 = matrix[0][2] * matrix[2][0];
1093 field_type det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
1094 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
1095 field_type t17 = 1.0/det;
1097 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
1098 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
1099 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
1100 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
1101 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
1102 inverse[2][1] = -(t6-t10) * t17;
1103 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
1104 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
1105 inverse[2][2] = (t4-t8) * t17;
1111 template<
class K,
int m,
int n,
int p >
1118 for( size_type i = 0; i < m; ++i )
1120 for( size_type j = 0; j < p; ++j )
1122 ret[ i ][ j ] = K( 0 );
1123 for( size_type k = 0; k < n; ++k )
1124 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
1130 template <
typename K,
int rows,
int cols>
1135 for(size_type i=0; i<cols(); i++)
1136 for(size_type j=0; j<cols(); j++)
1139 for(size_type k=0; k<rows(); k++)
1140 ret[i][j]+=matrix[k][i]*matrix[k][j];
1146 template <
typename MAT,
typename V1,
typename V2>
1153 for(size_type i=0; i<matrix.
rows(); ++i)
1156 for(size_type j=0; j<matrix.
cols(); ++j)
1158 ret[i] += matrix[i][j]*x[j];
1165 template <
typename K,
int rows,
int cols>
1170 for(size_type i=0; i<cols(); ++i)
1173 for(size_type j=0; j<rows(); ++j)
1174 ret[i] += matrix[j][i]*x[j];
1179 template <
typename K,
int rows,
int cols>
1180 static inline FieldVector<K,rows> mult(
const FieldMatrix<K,rows,cols> &matrix,
const FieldVector<K,cols> & x)
1182 FieldVector<K,rows> ret;
1188 template <
typename K,
int rows,
int cols>
1189 static inline FieldVector<K,cols> multTransposed(
const FieldMatrix<K,rows,cols> &matrix,
const FieldVector<K,rows> & x)
1191 FieldVector<K,cols> ret;
1192 multAssignTransposed( matrix, x, ret );
1200 template<
typename MAT>
1204 s << a[i] << std::endl;
#define DUNE_ASSERT_BOUNDS(cond)
Definition boundschecking.hh:20
Implements a vector constructed from a given type representing a field and a compile-time given size.
Various precision settings for calculations with FieldMatrix and FieldVector.
Some useful basic math stuff.
A few common exception classes.
A free function to provide the demangled class name of a given object or type as a string.
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition unused.hh:18
std::ostream & operator<<(std::ostream &s, const std::array< T, N > &e)
Output operator for array.
Definition array.hh:28
#define DUNE_THROW(E, m)
Definition exceptions.hh:216
Dune namespace.
Definition alignment.hh:11
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
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition densematrix.hh:1147
constexpr auto size(const Dune::FieldVector< T, i > *, const PriorityTag< 5 > &) -> decltype(std::integral_constant< std::size_t, i >())
Definition hybridutilities.hh:22
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition apply.hh:54
A dense n x m matrix.
Definition densematrix.hh:135
DenseMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition densematrix.hh:305
ConstIterator const_iterator
typedef for stl compliant access
Definition densematrix.hh:233
bool operator!=(const DenseMatrix< Other > &y) const
Binary matrix incomparison.
Definition densematrix.hh:334
void mv(const X &x, Y &y) const
y = A x
Definition densematrix.hh:344
Traits::value_type field_type
export the type representing the field
Definition densematrix.hh:152
void solve(V &x, const V &b) const
Solve system A x = b.
size_type cols() const
number of columns
Definition densematrix.hh:670
ConstIterator beforeEnd() const
Definition densematrix.hh:253
bool operator==(const DenseMatrix< Other > &y) const
Binary matrix comparison.
Definition densematrix.hh:324
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition densematrix.hh:422
DenseMatrix & operator=(const RHS &rhs)
Definition densematrix.hh:268
DenseMatrix & operator+=(const DenseMatrix< Other > &y)
vector space addition
Definition densematrix.hh:278
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition densematrix.hh:514
std::remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition densematrix.hh:237
ConstIterator beforeBegin() const
Definition densematrix.hh:260
Traits::value_type block_type
export the type representing the components
Definition densematrix.hh:155
void mtv(const X &x, Y &y) const
y = A^T x
Definition densematrix.hh:361
DenseMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition densematrix.hh:297
DenseMatrix & operator-=(const DenseMatrix< Other > &y)
vector space subtraction
Definition densematrix.hh:288
size_type size() const
size method (number of rows)
Definition densematrix.hh:189
size_type M() const
number of columns
Definition densematrix.hh:658
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition densematrix.hh:600
Iterator end()
end iterator
Definition densematrix.hh:211
Iterator beforeBegin()
Definition densematrix.hh:225
Iterator beforeEnd()
Definition densematrix.hh:218
Traits::value_type value_type
export the type representing the field
Definition densematrix.hh:149
void usmv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition densematrix.hh:444
void usmhv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition densematrix.hh:468
size_type rows() const
number of rows
Definition densematrix.hh:664
void mmv(const X &x, Y &y) const
y -= A x
Definition densematrix.hh:411
void invert()
Compute inverse.
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition densematrix.hh:582
void usmtv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition densematrix.hh:456
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition densematrix.hh:433
Traits::derived_type derived_type
type of derived matrix class
Definition densematrix.hh:146
row_reference operator[](size_type i)
random access
Definition densematrix.hh:178
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition densematrix.hh:678
DenseMatrix & axpy(const field_type &k, const DenseMatrix< Other > &y)
vector space axpy operation (*this += k y)
Definition densematrix.hh:314
Iterator RowIterator
rename the iterators for easier access
Definition densematrix.hh:200
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition densematrix.hh:481
void umv(const X &x, Y &y) const
y += A x
Definition densematrix.hh:378
DenseIterator< const DenseMatrix, const row_type, const_row_reference > ConstIterator
Iterator class for sequential access.
Definition densematrix.hh:231
size_type N() const
number of rows
Definition densematrix.hh:652
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition densematrix.hh:499
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition densematrix.hh:161
Traits::size_type size_type
The type used for the index access and size operation.
Definition densematrix.hh:158
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:167
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition densematrix.hh:489
std::remove_reference< row_reference >::type::Iterator ColIterator
rename the iterators for easier access
Definition densematrix.hh:202
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition densematrix.hh:164
Iterator iterator
typedef for stl compliant access
Definition densematrix.hh:198
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition densematrix.hh:235
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:196
void umtv(const X &x, Y &y) const
y += A^T x
Definition densematrix.hh:389
ConstIterator begin() const
begin iterator
Definition densematrix.hh:240
@ blocklevel
The number of block levels we contain. This is 1.
Definition densematrix.hh:172
Iterator begin()
begin iterator
Definition densematrix.hh:205
void umhv(const X &x, Y &y) const
y += A^H x
Definition densematrix.hh:400
ConstIterator end() const
end iterator
Definition densematrix.hh:246
const FieldTraits< typenameDenseMatVecTraits< M >::value_type >::real_type real_type
Definition densematrix.hh:30
const FieldTraits< typenameDenseMatVecTraits< M >::value_type >::field_type field_type
Definition densematrix.hh:29
A dense n x m matrix.
Definition fmatrix.hh:67
Base::size_type size_type
Definition fmatrix.hh:80
vector space out of a tensor product of fields.
Definition fvector.hh:93
constexpr FieldVector()
Constructor making default-initialized vector.
Definition fvector.hh:113
you have to specialize this structure for any type that should be assignable to a DenseMatrix
Definition densematrix.hh:73
Error thrown if operations of a FieldMatrix fail.
Definition densematrix.hh:121
Interface for a class of dense vectors over a given field.
Definition densevector.hh:235
size_type size() const
size method
Definition densevector.hh:297
size_type N() const
number of blocks in the vector (are of size 1 here)
Definition densevector.hh:678
Generic iterator class for dense vector and matrix implementations.
Definition densevector.hh:128
Default exception class for mathematical errors.
Definition exceptions.hh:239
T field_type
export the type representing the field
Definition ftraits.hh:26
T real_type
export the type representing the real type of the field
Definition ftraits.hh:28
Definition matvectraits.hh:29
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