3 #ifndef DUNE_DENSEMATRIX_HH 4 #define DUNE_DENSEMATRIX_HH 44 static typename V::size_type
size(
const V & v) {
return v.size(); }
47 template<
class K,
int N>
51 static typename V::size_type
size(
const V & v) {
return N; }
73 template<
class DenseMatrix,
class RHS >
80 template<
class DenseMatrix,
class RHS,
class =
void >
84 template<
class DenseMatrix,
class RHS >
91 std::fill( denseMatrix.
begin(), denseMatrix.
end(),
static_cast< field_type
>( rhs ) );
95 template<
class DenseMatrix,
class RHS >
104 typename RHS::const_iterator sIt = std::begin(rhs);
105 for(; sIt != std::end(rhs); ++tIt, ++sIt)
106 std::copy(std::begin(*sIt), std::end(*sIt), std::begin(*tIt));
114 template<
class DenseMatrix,
class RHS >
116 :
public Impl::DenseMatrixAssigner< DenseMatrix, RHS >
123 template<
class DenseMatrix,
class RHS >
126 std::false_type hasDenseMatrixAssigner ( ... );
130 template<
class DenseMatrix,
class RHS >
131 struct HasDenseMatrixAssigner
132 :
public decltype( Impl::hasDenseMatrixAssigner( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) )
135 #endif // #ifndef DOXYGEN 152 template<
typename MAT>
158 MAT & asImp() {
return static_cast<MAT&
>(*this); }
159 const MAT & asImp()
const {
return static_cast<const MAT&
>(*this); }
197 row_reference operator[] ( size_type i )
199 return asImp().mat_access(i);
202 const_row_reference operator[] ( size_type i )
const 204 return asImp().mat_access(i);
221 typedef typename std::remove_reference<row_reference>::type::Iterator
ColIterator;
226 return Iterator(*
this,0);
232 return Iterator(*
this,rows());
239 return Iterator(*
this,rows()-1);
246 return Iterator(*
this,-1);
256 typedef typename std::remove_reference<const_row_reference>::type::ConstIterator
ConstColIterator;
261 return ConstIterator(*
this,0);
265 ConstIterator
end ()
const 267 return ConstIterator(*
this,rows());
274 return ConstIterator(*
this,rows()-1);
281 return ConstIterator(*
this,-1);
286 template< class RHS, class = std::enable_if_t< HasDenseMatrixAssigner< MAT, RHS >::value > >
296 template <
class Other>
300 for (size_type i=0; i<rows(); i++)
306 template <
class Other>
310 for (size_type i=0; i<rows(); i++)
318 for (size_type i=0; i<rows(); i++)
326 for (size_type i=0; i<rows(); i++)
332 template <
class Other>
336 for( size_type i = 0; i < rows(); ++i )
337 (*
this)[ i ].axpy( k, y[ i ] );
342 template <
class Other>
346 for (size_type i=0; i<rows(); i++)
347 if ((*
this)[i]!=y[i])
352 template <
class Other>
362 template<
class X,
class Y>
363 void mv (
const X& x, Y& y)
const 370 for (size_type i=0; i<rows(); ++i)
373 for (size_type j=0; j<cols(); j++)
374 y[i] += (*
this)[i][j] * x[j];
379 template<
class X,
class Y >
380 void mtv (
const X &x, Y &y )
const 387 for(size_type i = 0; i < cols(); ++i)
390 for(size_type j = 0; j < rows(); ++j)
391 y[i] += (*
this)[j][i] * x[j];
396 template<
class X,
class Y>
397 void umv (
const X& x, Y& y)
const 401 for (size_type i=0; i<rows(); i++)
402 for (size_type j=0; j<cols(); j++)
403 y[i] += (*
this)[i][j] * x[j];
407 template<
class X,
class Y>
408 void umtv (
const X& x, Y& y)
const 412 for (size_type i=0; i<rows(); i++)
413 for (size_type j=0; j<cols(); j++)
414 y[j] += (*
this)[i][j]*x[i];
418 template<
class X,
class Y>
419 void umhv (
const X& x, Y& y)
const 423 for (size_type i=0; i<rows(); i++)
424 for (size_type j=0; j<cols(); j++)
429 template<
class X,
class Y>
430 void mmv (
const X& x, Y& y)
const 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 mmtv (
const X& x, Y& y)
const 445 for (size_type i=0; i<rows(); i++)
446 for (size_type j=0; j<cols(); j++)
447 y[j] -= (*
this)[i][j]*x[i];
451 template<
class X,
class Y>
452 void mmhv (
const X& x, Y& y)
const 456 for (size_type i=0; i<rows(); i++)
457 for (size_type j=0; j<cols(); j++)
462 template<
class X,
class Y>
464 const X& x, Y& y)
const 468 for (size_type i=0; i<rows(); i++)
469 for (size_type j=0; j<cols(); j++)
470 y[i] += alpha * (*
this)[i][j] * x[j];
474 template<
class X,
class Y>
476 const X& x, Y& y)
const 480 for (size_type i=0; i<rows(); i++)
481 for (size_type j=0; j<cols(); j++)
482 y[j] += alpha*(*
this)[i][j]*x[i];
486 template<
class X,
class Y>
488 const X& x, Y& y)
const 492 for (size_type i=0; i<rows(); i++)
493 for (size_type j=0; j<cols(); j++)
503 for (size_type i=0; i<rows(); ++i) sum += (*
this)[i].two_norm2();
504 return fvmeta::sqrt(sum);
511 for (size_type i=0; i<rows(); ++i) sum += (*
this)[i].two_norm2();
516 template <
typename vt = value_type,
517 typename std::enable_if<!has_nan<vt>::value,
int>::type = 0>
523 for (
auto const &x : *
this) {
531 template <
typename vt = value_type,
532 typename std::enable_if<!has_nan<vt>::value,
int>::type = 0>
538 for (
auto const &x : *
this) {
546 template <
typename vt = value_type,
547 typename std::enable_if<has_nan<vt>::value,
int>::type = 0>
554 for (
auto const &x : *
this) {
564 template <
typename vt = value_type,
565 typename std::enable_if<has_nan<vt>::value,
int>::type = 0>
572 for (
auto const &x : *
this) {
588 void solve (V& x,
const V& b)
const;
597 field_type determinant ()
const;
600 template<
typename M2>
607 for (size_type i=0; i<rows(); i++)
608 for (size_type j=0; j<cols(); j++) {
610 for (size_type k=0; k<rows(); k++)
611 (*
this)[i][j] += M[i][k]*C[k][j];
618 template<
typename M2>
625 for (size_type i=0; i<rows(); i++)
626 for (size_type j=0; j<cols(); j++) {
628 for (size_type k=0; k<cols(); k++)
629 (*
this)[i][j] += C[i][k]*M[k][j];
641 for (size_type i=0; i<l; i++) {
642 for (size_type j=0; j<cols(); j++) {
644 for (size_type k=0; k<rows(); k++)
645 C[i][j] += M[i][k]*(*
this)[k][j];
657 for (size_type i=0; i<rows(); i++) {
658 for (size_type j=0; j<l; j++) {
660 for (size_type k=0; k<cols(); k++)
661 C[i][j] += (*
this)[i][k]*M[k][j];
685 return asImp().mat_rows();
691 return asImp().mat_cols();
697 bool exists (size_type i, size_type j)
const 711 ElimPivot(std::vector<size_type> & pivot);
713 void swap(
int i,
int j);
716 void operator()(
const T&,
int,
int)
719 std::vector<size_type> & pivot_;
727 void swap(
int i,
int j);
729 void operator()(
const typename V::field_type& factor,
int k,
int i);
736 ElimDet(field_type&
sign) : sign_(sign)
742 void operator()(
const field_type&,
int,
int)
754 template<
typename MAT>
758 typedef typename std::vector<size_type>::size_type
size_type;
759 for(size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
762 template<
typename MAT>
768 template<
typename MAT>
774 template<
typename MAT>
778 std::swap((*rhs_)[i], (*rhs_)[j]);
781 template<
typename MAT>
786 (*rhs_)[k] -= factor*(*rhs_)[i];
788 template<
typename MAT>
789 template<
typename Func>
810 if ((abs=fvmeta::absreal(A[k][i]))>pivmax)
812 pivmax = abs; imax = k;
817 std::swap(A[i][j],A[imax][j]);
823 if (pivmax<singthres)
832 A[k][j] -= factor*A[i][j];
838 template<
typename MAT>
848 #ifdef DUNE_FMatrix_WITH_CHECKING 852 x[0] = b[0]/(*this)[0][0];
855 else if (rows()==2) {
857 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
858 #ifdef DUNE_FMatrix_WITH_CHECKING 864 x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
865 x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
868 else if (rows()==3) {
871 #ifdef DUNE_FMatrix_WITH_CHECKING 876 x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
877 - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
878 + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
880 x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
881 - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
882 + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
884 x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
885 - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
886 + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
896 luDecomposition(A, elim);
899 for(
int i=rows()-1; i>=0; i--) {
901 rhs[i] -= A[i][j]*x[j];
902 x[i] = rhs[i]/A[i][i];
907 template<
typename MAT>
916 #ifdef DUNE_FMatrix_WITH_CHECKING 920 (*this)[0][0] =
field_type( 1 ) / (*this)[0][0];
923 else if (rows()==2) {
925 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
926 #ifdef DUNE_FMatrix_WITH_CHECKING 933 (*this)[0][0] = (*this)[1][1]*detinv;
934 (*this)[0][1] = -(*this)[0][1]*detinv;
935 (*this)[1][0] = -(*this)[1][0]*detinv;
936 (*this)[1][1] = temp*detinv;
943 K t4 = (*this)[0][0] * (*this)[1][1];
944 K t6 = (*this)[0][0] * (*this)[1][2];
945 K t8 = (*this)[0][1] * (*this)[1][0];
946 K t10 = (*this)[0][2] * (*this)[1][0];
947 K t12 = (*this)[0][1] * (*this)[2][0];
948 K t14 = (*this)[0][2] * (*this)[2][0];
950 K det = (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
951 t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
954 K matrix01 = (*this)[0][1];
955 K matrix00 = (*this)[0][0];
956 K matrix10 = (*this)[1][0];
957 K matrix11 = (*this)[1][1];
959 (*this)[0][0] = ((*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1])*t17;
960 (*this)[0][1] = -((*this)[0][1] * (*this)[2][2] - (*this)[0][2] * (*this)[2][1])*t17;
961 (*this)[0][2] = (matrix01 * (*this)[1][2] - (*this)[0][2] * (*this)[1][1])*t17;
962 (*this)[1][0] = -((*this)[1][0] * (*this)[2][2] - (*this)[1][2] * (*this)[2][0])*t17;
963 (*this)[1][1] = (matrix00 * (*this)[2][2] - t14) * t17;
964 (*this)[1][2] = -(t6-t10) * t17;
965 (*this)[2][0] = (matrix10 * (*this)[2][1] - matrix11 * (*this)[2][0]) * t17;
966 (*this)[2][1] = -(matrix00 * (*this)[2][1] - t12) * t17;
967 (*this)[2][2] = (t4-t8) * t17;
972 std::vector<size_type> pivot(rows());
973 luDecomposition(A, ElimPivot(pivot));
987 (*
this)[i][k] -= L[i][j]*(*this)[j][k];
994 (*
this)[i][k] -= U[i][j]*(*this)[j][k];
995 (*this)[i][k] /= U[i][i];
1003 std::swap((*
this)[j][pivot[i]], (*this)[j][i]);
1009 template<
typename MAT>
1018 return (*
this)[0][0];
1021 return (*
this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
1025 field_type t4 = (*this)[0][0] * (*this)[1][1];
1026 field_type t6 = (*this)[0][0] * (*this)[1][2];
1027 field_type t8 = (*this)[0][1] * (*this)[1][0];
1028 field_type t10 = (*this)[0][2] * (*this)[1][0];
1029 field_type t12 = (*this)[0][1] * (*this)[2][0];
1030 field_type t14 = (*this)[0][2] * (*this)[2][0];
1032 return (t4*(*
this)[2][2]-t6*(*
this)[2][1]-t8*(*
this)[2][2]+
1033 t10*(*
this)[2][1]+t12*(*
this)[1][2]-t14*(*
this)[1][1]);
1041 luDecomposition(A, ElimDet(det));
1054 namespace DenseMatrixHelp {
1057 template <
typename MAT,
typename V1,
typename V2>
1064 for(size_type i=0; i<matrix.
rows(); ++i)
1067 for(size_type j=0; j<matrix.
cols(); ++j)
1069 ret[i] += matrix[i][j]*x[j];
1075 template <
typename K,
int rows,
int cols>
1081 for(size_type i=0; i<cols(); ++i)
1084 for(size_type j=0; j<rows(); ++j)
1085 ret[i] += matrix[j][i]*x[j];
1090 template <
typename K,
int rows,
int cols>
1099 template <
typename K,
int rows,
int cols>
1111 template<
typename MAT>
1112 std::ostream& operator<< (std::ostream& s, const DenseMatrix<MAT>& a)
1115 s << a[i] << std::endl;
ConstIterator end() const
end iterator
Definition: densematrix.hh:265
EnableIfInterOperable< T1, T2, bool >::type operator==(const ForwardIteratorFacade< T1, V1, R1, D > &lhs, const ForwardIteratorFacade< T2, V2, R2, D > &rhs)
Checks for equality.
Definition: iteratorfacades.hh:233
Implements a vector constructed from a given type representing a field and a compile-time given size...
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:126
Various precision settings for calculations with FieldMatrix and FieldVector.
std::remove_reference< row_reference >::type::Iterator ColIterator
rename the iterators for easier access
Definition: densematrix.hh:221
Traits::value_type field_type
export the type representing the field
Definition: densematrix.hh:171
DenseMatrix & axpy(const field_type &k, const DenseMatrix< Other > &y)
vector space axpy operation (*this += k y)
Definition: densematrix.hh:333
void umv(const X &x, Y &y) const
y += A x
Definition: densematrix.hh:397
ConstIterator const_iterator
typedef for stl compliant access
Definition: densematrix.hh:252
size_type size() const
size method
Definition: densevector.hh:297
Traits::value_type value_type
export the type representing the field
Definition: densematrix.hh:168
Default exception class for mathematical errors.
Definition: exceptions.hh:239
DenseIterator< DenseMatrix, row_type, row_reference > Iterator
Iterator class for sequential access.
Definition: densematrix.hh:215
container_type::size_type size_type
Definition: dynmatrix.hh:43
void solve(V &x, const V &b) const
Solve system A x = b.
void umhv(const X &x, Y &y) const
y += A^H x
Definition: densematrix.hh:419
static FieldVector< K, cols > multTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x)
calculates ret = matrix^T * x
Definition: fmatrix.hh:505
std::remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: densematrix.hh:256
field_type determinant() const
calculates the determinant of this matrix
ConstIterator beforeBegin() const
Definition: densematrix.hh:279
size_type size() const
size method (number of rows)
Definition: densematrix.hh:208
decltype(auto) apply(F &&f, ArgTuple &&args)
Apply function with arguments given as tuple.
Definition: apply.hh:54
ConstIterator beforeEnd() const
Definition: densematrix.hh:272
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:619
void usmhv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: densematrix.hh:487
const FieldTraits< typename DenseMatVecTraits< M >::value_type >::field_type field_type
Definition: densematrix.hh:30
A dense n x m matrix.
Definition: densematrix.hh:38
Construct a matrix with a dynamic size.
Definition: dynmatrix.hh:29
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: densematrix.hh:441
K value_type
Definition: dynmatrix.hh:42
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:482
Dune namespace.
Definition: alignment.hh:10
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:140
Traits for type conversions and type information.
Some useful basic math stuff.
DenseIterator< const DenseMatrix, const row_type, const_row_reference > ConstIterator
Iterator class for sequential access.
Definition: densematrix.hh:250
void umtv(const X &x, Y &y) const
y += A^T x
Definition: densematrix.hh:408
void usmtv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: densematrix.hh:475
A few common exception classes.
Macro for wrapping boundary checks.
void mv(const X &x, Y &y) const
y = A x
Definition: densematrix.hh:363
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:183
A free function to provide the demangled class name of a given object or type as a string...
#define DUNE_THROW(E, m)
Definition: exceptions.hh:216
Definition of the DUNE_UNUSED macro for the case that config.h is not available.
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:103
size_type cols() const
number of columns
Definition: densematrix.hh:689
ConstIterator begin() const
begin iterator
Definition: densematrix.hh:259
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: densematrix.hh:500
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:177
Iterator RowIterator
rename the iterators for easier access
Definition: densematrix.hh:219
T real_type
export the type representing the real type of the field
Definition: ftraits.hh:28
void invert()
Compute inverse.
Definition: matvectraits.hh:29
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densematrix.hh:533
const FieldTraits< typename DenseMatVecTraits< M >::value_type >::real_type real_type
Definition: densematrix.hh:31
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: densematrix.hh:508
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition: densematrix.hh:1058
static FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition: fmatrix.hh:496
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: densematrix.hh:452
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: densematrix.hh:254
Precisions for calculations with FieldMatrix and FieldVector.
Definition: precision.hh:23
Traits::derived_type derived_type
type of derived matrix class
Definition: densematrix.hh:165
EnableIfInterOperable< T1, T2, bool >::type operator!=(const ForwardIteratorFacade< T1, V1, R1, D > &lhs, const ForwardIteratorFacade< T2, V2, R2, D > &rhs)
Checks for inequality.
Definition: iteratorfacades.hh:255
size_type M() const
number of columns
Definition: densematrix.hh:677
you have to specialize this structure for any type that should be assignable to a DenseMatrix ...
Definition: densematrix.hh:74
Iterator iterator
typedef for stl compliant access
Definition: densematrix.hh:217
vector space out of a tensor product of fields.
Definition: densematrix.hh:39
size_type rows() const
number of rows
Definition: densematrix.hh:683
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:119
A dense n x m matrix.
Definition: densematrix.hh:25
void usmv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition: densematrix.hh:463
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:18
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:180
Interface for a class of dense vectors over a given field.
Definition: densevector.hh:19
Definition: ftraits.hh:23
constexpr auto size(const Dune::FieldVector< T, i > *, const PriorityTag< 5 > &) -> decltype(std::integral_constant< std::size_t, i >())
Definition: hybridutilities.hh:22
size_type N() const
number of rows
Definition: densematrix.hh:671
Iterator begin()
begin iterator
Definition: densematrix.hh:224
Iterator beforeEnd()
Definition: densematrix.hh:237
Iterator end()
end iterator
Definition: densematrix.hh:230
Iterator beforeBegin()
Definition: densematrix.hh:244
T field_type
export the type representing the field
Definition: ftraits.hh:26
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: densematrix.hh:518
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition: densematrix.hh:601
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:186
void mmv(const X &x, Y &y) const
y -= A x
Definition: densematrix.hh:430
void mtv(const X &x, Y &y) const
y = A^T x
Definition: densematrix.hh:380
Traits::value_type block_type
export the type representing the components
Definition: densematrix.hh:174
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition: boundschecking.hh:28
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: densematrix.hh:697
Base::size_type size_type
Definition: fmatrix.hh:81