5#ifndef DUNE_DENSEMATRIX_HH
6#define DUNE_DENSEMATRIX_HH
28 template<
typename M>
class DenseMatrix;
31 struct FieldTraits< DenseMatrix<M> >
33 typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::field_type field_type;
34 typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::real_type real_type;
37 template<
class K,
int N,
int M>
class FieldMatrix;
38 template<
class K,
int N>
class FieldVector;
58 template<
class DenseMatrix,
class RHS >
65 template<
class DenseMatrix,
class RHS,
class =
void >
69 template<
class DenseMatrix,
class RHS >
73 static void apply (
DenseMatrix &denseMatrix,
const RHS &rhs )
76 std::fill( denseMatrix.
begin(), denseMatrix.
end(),
static_cast< field_type
>( rhs ) );
80 template<
class DenseMatrix,
class RHS >
81 class DenseMatrixAssigner< DenseMatrix, RHS,
std::enable_if_t< !std::is_same< typename RHS::const_iterator, void >::value
82 && std::is_convertible< typename RHS::const_iterator::value_type, typename DenseMatrix::iterator::value_type >::value > >
85 static void apply ( DenseMatrix &denseMatrix,
const RHS &rhs )
90 typename RHS::const_iterator sIt = std::begin(rhs);
91 for(; sIt != std::end(rhs); ++tIt, ++sIt)
92 std::copy(std::begin(*sIt), std::end(*sIt), std::begin(*tIt));
100 template<
class DenseMatrix,
class RHS >
101 struct DenseMatrixAssigner
102 :
public Impl::DenseMatrixAssigner< DenseMatrix, RHS >
109 template<
class DenseMatrix,
class RHS >
112 std::false_type hasDenseMatrixAssigner ( ... );
116 template<
class DenseMatrix,
class RHS >
117 struct HasDenseMatrixAssigner
118 :
public decltype( Impl::hasDenseMatrixAssigner( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) )
138 template<
typename MAT>
141 typedef DenseMatVecTraits<MAT> Traits;
144 constexpr MAT & asImp() {
return static_cast<MAT&
>(*this); }
145 constexpr const MAT & asImp()
const {
return static_cast<const MAT&
>(*this); }
191 return asImp().mat_access(i);
196 return asImp().mat_access(i);
213 typedef typename std::remove_reference<row_reference>::type::Iterator
ColIterator;
248 typedef typename std::remove_reference<const_row_reference>::type::ConstIterator
ConstColIterator;
278 template< class RHS, class = std::enable_if_t< HasDenseMatrixAssigner< MAT, RHS >::value > >
288 template <
class Other>
301 using idx_type =
typename decltype(result)
::size_type;
303 for (idx_type i = 0; i <
rows(); ++i)
304 for (idx_type j = 0; j <
cols(); ++j)
305 result[i][j] = - asImp()[i][j];
311 template <
class Other>
337 template <
class Other>
342 (*
this)[ i ].axpy( a, x[ i ] );
347 template <
class Other>
352 if ((*
this)[i]!=x[i])
357 template <
class Other>
367 template<
class X,
class Y>
368 void mv (
const X& x, Y& y)
const
370 auto&& xx = Impl::asVector(x);
371 auto&& yy = Impl::asVector(y);
376 using y_field_type =
typename FieldTraits<Y>::field_type;
379 yy[i] = y_field_type(0);
381 yy[i] += (*
this)[i][j] * xx[j];
386 template<
class X,
class Y >
387 void mtv (
const X &x, Y &y )
const
389 auto&& xx = Impl::asVector(x);
390 auto&& yy = Impl::asVector(y);
395 using y_field_type =
typename FieldTraits<Y>::field_type;
398 yy[i] = y_field_type(0);
400 yy[i] += (*
this)[j][i] * xx[j];
405 template<
class X,
class Y>
406 void umv (
const X& x, Y& y)
const
408 auto&& xx = Impl::asVector(x);
409 auto&& yy = Impl::asVector(y);
414 yy[i] += (*
this)[i][j] * xx[j];
418 template<
class X,
class Y>
419 void umtv (
const X& x, Y& y)
const
421 auto&& xx = Impl::asVector(x);
422 auto&& yy = Impl::asVector(y);
427 yy[j] += (*
this)[i][j]*xx[i];
431 template<
class X,
class Y>
432 void umhv (
const X& x, Y& y)
const
434 auto&& xx = Impl::asVector(x);
435 auto&& yy = Impl::asVector(y);
444 template<
class X,
class Y>
445 void mmv (
const X& x, Y& y)
const
447 auto&& xx = Impl::asVector(x);
448 auto&& yy = Impl::asVector(y);
453 yy[i] -= (*
this)[i][j] * xx[j];
457 template<
class X,
class Y>
458 void mmtv (
const X& x, Y& y)
const
460 auto&& xx = Impl::asVector(x);
461 auto&& yy = Impl::asVector(y);
466 yy[j] -= (*
this)[i][j]*xx[i];
470 template<
class X,
class Y>
471 void mmhv (
const X& x, Y& y)
const
473 auto&& xx = Impl::asVector(x);
474 auto&& yy = Impl::asVector(y);
483 template<
class X,
class Y>
484 void usmv (
const typename FieldTraits<Y>::field_type & alpha,
485 const X& x, Y& y)
const
487 auto&& xx = Impl::asVector(x);
488 auto&& yy = Impl::asVector(y);
493 yy[i] += alpha * (*
this)[i][j] * xx[j];
497 template<
class X,
class Y>
498 void usmtv (
const typename FieldTraits<Y>::field_type & alpha,
499 const X& x, Y& y)
const
501 auto&& xx = Impl::asVector(x);
502 auto&& yy = Impl::asVector(y);
507 yy[j] += alpha*(*
this)[i][j]*xx[i];
511 template<
class X,
class Y>
512 void usmhv (
const typename FieldTraits<Y>::field_type & alpha,
513 const X& x, Y& y)
const
515 auto&& xx = Impl::asVector(x);
516 auto&& yy = Impl::asVector(y);
530 typename FieldTraits<value_type>::real_type sum=(0.0);
531 for (
size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
532 return fvmeta::sqrt(sum);
538 typename FieldTraits<value_type>::real_type sum=(0.0);
539 for (
size_type i=0; i<
rows(); ++i) sum += (*
this)[i].two_norm2();
545 typename std::enable_if<!HasNaN<vt>::value,
int>::type = 0>
547 using real_type =
typename FieldTraits<vt>::real_type;
551 for (
auto const &x : *
this) {
552 real_type
const a = x.one_norm();
560 typename std::enable_if<!HasNaN<vt>::value,
int>::type = 0>
562 using real_type =
typename FieldTraits<vt>::real_type;
566 for (
auto const &x : *
this) {
567 real_type
const a = x.one_norm_real();
575 typename std::enable_if<HasNaN<vt>::value,
int>::type = 0>
577 using real_type =
typename FieldTraits<vt>::real_type;
582 for (
auto const &x : *
this) {
583 real_type
const a = x.one_norm();
592 typename std::enable_if<HasNaN<vt>::value,
int>::type = 0>
594 using real_type =
typename FieldTraits<vt>::real_type;
599 for (
auto const &x : *
this) {
600 real_type
const a = x.one_norm_real();
613 template <
class V1,
class V2>
614 void solve (V1& x,
const V2& b,
bool doPivoting =
true)
const;
626 template<
typename M2>
637 (*
this)[i][j] +=
M[i][k]*C[k][j];
644 template<
typename M2>
655 (*
this)[i][j] += C[i][k]*
M[k][j];
671 C[i][j] +=
M[i][k]*(*
this)[k][j];
679 FieldMatrix<K,rows,l> rightmultiplyany (
const FieldMatrix<K,cols,l>&
M)
const
681 FieldMatrix<K,rows,l> C;
687 C[i][j] += (*
this)[i][k]*
M[k][j];
711 return asImp().mat_rows();
717 return asImp().mat_cols();
735 ElimPivot(std::vector<simd_index_type> & pivot);
737 void swap(std::size_t i, simd_index_type j);
740 void operator()(
const T&,
int,
int)
743 std::vector<simd_index_type> & pivot_;
751 void swap(std::size_t i, simd_index_type j);
753 void operator()(
const typename V::field_type& factor,
int k,
int i);
763 void swap(std::size_t i, simd_index_type j)
815 template<
class Func,
class Mask>
817 Mask &nonsingularLanes,
bool throwEarly,
bool doPivoting);
821 template<
typename MAT>
825 typedef typename std::vector<size_type>::size_type
size_type;
826 for(
size_type i=0; i < pivot_.size(); ++i) pivot_[i]=i;
829 template<
typename MAT>
830 void DenseMatrix<MAT>::ElimPivot::swap(std::size_t i, simd_index_type j)
833 Simd::cond(Simd::Scalar<simd_index_type>(i) == j, pivot_[i], j);
836 template<
typename MAT>
838 DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
842 template<
typename MAT>
844 void DenseMatrix<MAT>::Elim<V>::swap(std::size_t i, simd_index_type j)
854 template<
typename MAT>
856 void DenseMatrix<MAT>::
857 Elim<V>::operator()(
const typename V::field_type& factor,
int k,
int i)
859 (*rhs_)[k] -= factor*(*rhs_)[i];
862 template<
typename MAT>
863 template<
typename Func,
class Mask>
866 bool throwEarly,
bool doPivoting)
871 typedef typename FieldTraits<value_type>::real_type real_type;
874 for (size_type i=0; i<A.rows(); i++)
876 real_type pivmax = fvmeta::absreal(A[i][i]);
881 simd_index_type imax=i;
882 for (size_type k=i+1; k<A.rows(); k++)
884 auto abs = fvmeta::absreal(A[k][i]);
885 auto mask = abs > pivmax;
890 for (size_type j=0; j<A.rows(); j++)
900 for(std::size_t l = 0; l <
Simd::lanes(A[i][j]); ++l)
908 nonsingularLanes = nonsingularLanes && (pivmax != real_type(0));
911 DUNE_THROW(FMatrixError,
"matrix is singular");
919 for (size_type k=i+1; k<A.rows(); k++)
923 field_type factor = A[k][i]/A[i][i];
925 for (size_type j=i+1; j<A.rows(); j++)
926 A[k][j] -= factor*A[i][j];
932 template<
typename MAT>
933 template <
class V1,
class V2>
936 using real_type =
typename FieldTraits<value_type>::real_type;
939 DUNE_THROW(FMatrixError,
"Can't solve for a " << rows() <<
"x" << cols() <<
" matrix!");
943#ifdef DUNE_FMatrix_WITH_CHECKING
946 DUNE_THROW(FMatrixError,
"matrix is singular");
948 x[0] = b[0]/(*this)[0][0];
951 else if (rows()==2) {
953 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
954#ifdef DUNE_FMatrix_WITH_CHECKING
957 DUNE_THROW(FMatrixError,
"matrix is singular");
959 detinv = real_type(1.0)/detinv;
961 x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
962 x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
965 else if (rows()==3) {
967 field_type d = determinant(doPivoting);
968#ifdef DUNE_FMatrix_WITH_CHECKING
971 DUNE_THROW(FMatrixError,
"matrix is singular");
974 x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
975 - b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
976 + b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
978 x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
979 - (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
980 + (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
982 x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
983 - (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
984 + (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
992 AutonomousValue<MAT> A(asImp());
993 Simd::Mask<typename FieldTraits<value_type>::real_type>
994 nonsingularLanes(
true);
996 AutonomousValue<MAT>::luDecomposition(A, elim, nonsingularLanes,
true, doPivoting);
999 for(
int i=rows()-1; i>=0; i--) {
1000 for (size_type j=i+1; j<rows(); j++)
1001 rhs[i] -= A[i][j]*x[j];
1002 x[i] = rhs[i]/A[i][i];
1007 template<
typename MAT>
1010 using real_type =
typename FieldTraits<MAT>::real_type;
1015 DUNE_THROW(FMatrixError,
"Can't invert a " << rows() <<
"x" << cols() <<
" matrix!");
1019#ifdef DUNE_FMatrix_WITH_CHECKING
1022 DUNE_THROW(FMatrixError,
"matrix is singular");
1024 (*this)[0][0] = real_type( 1 ) / (*this)[0][0];
1027 else if (rows()==2) {
1029 field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
1030#ifdef DUNE_FMatrix_WITH_CHECKING
1033 DUNE_THROW(FMatrixError,
"matrix is singular");
1035 detinv = real_type( 1 ) / detinv;
1037 field_type temp=(*this)[0][0];
1038 (*this)[0][0] = (*this)[1][1]*detinv;
1039 (*this)[0][1] = -(*this)[0][1]*detinv;
1040 (*this)[1][0] = -(*this)[1][0]*detinv;
1041 (*this)[1][1] = temp*detinv;
1046 using K = field_type;
1048 K t4 = (*this)[0][0] * (*this)[1][1];
1049 K t6 = (*this)[0][0] * (*this)[1][2];
1050 K t8 = (*this)[0][1] * (*this)[1][0];
1051 K t10 = (*this)[0][2] * (*this)[1][0];
1052 K t12 = (*this)[0][1] * (*this)[2][0];
1053 K t14 = (*this)[0][2] * (*this)[2][0];
1055 K det = (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
1056 t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
1059 K matrix01 = (*this)[0][1];
1060 K matrix00 = (*this)[0][0];
1061 K matrix10 = (*this)[1][0];
1062 K matrix11 = (*this)[1][1];
1064 (*this)[0][0] = ((*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1])*t17;
1065 (*this)[0][1] = -((*this)[0][1] * (*this)[2][2] - (*this)[0][2] * (*this)[2][1])*t17;
1066 (*this)[0][2] = (matrix01 * (*this)[1][2] - (*this)[0][2] * (*this)[1][1])*t17;
1067 (*this)[1][0] = -((*this)[1][0] * (*this)[2][2] - (*this)[1][2] * (*this)[2][0])*t17;
1068 (*this)[1][1] = (matrix00 * (*this)[2][2] - t14) * t17;
1069 (*this)[1][2] = -(t6-t10) * t17;
1070 (*this)[2][0] = (matrix10 * (*this)[2][1] - matrix11 * (*this)[2][0]) * t17;
1071 (*this)[2][1] = -(matrix00 * (*this)[2][1] - t12) * t17;
1072 (*this)[2][2] = (t4-t8) * t17;
1077 AutonomousValue<MAT> A(asImp());
1078 std::vector<simd_index_type> pivot(rows());
1079 Simd::Mask<typename FieldTraits<value_type>::real_type>
1080 nonsingularLanes(
true);
1081 AutonomousValue<MAT>::luDecomposition(A, ElimPivot(pivot), nonsingularLanes,
true, doPivoting);
1088 for(size_type i=0; i<rows(); ++i)
1092 for (size_type i=0; i<rows(); i++)
1093 for (size_type j=0; j<i; j++)
1094 for (size_type k=0; k<rows(); k++)
1095 (*
this)[i][k] -= L[i][j]*(*this)[j][k];
1098 for (size_type i=rows(); i>0;) {
1100 for (size_type k=0; k<rows(); k++) {
1101 for (size_type j=i+1; j<rows(); j++)
1102 (*
this)[i][k] -= U[i][j]*(*this)[j][k];
1103 (*this)[i][k] /= U[i][i];
1107 for(size_type i=rows(); i>0; ) {
1109 for(std::size_t l = 0; l <
Simd::lanes((*
this)[0][0]); ++l)
1113 for(size_type j=0; j<rows(); ++j)
1122 template<
typename MAT>
1128 DUNE_THROW(FMatrixError,
"There is no determinant for a " << rows() <<
"x" << cols() <<
" matrix!");
1131 return (*
this)[0][0];
1134 return (*
this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
1138 field_type t4 = (*this)[0][0] * (*this)[1][1];
1139 field_type t6 = (*this)[0][0] * (*this)[1][2];
1140 field_type t8 = (*this)[0][1] * (*this)[1][0];
1141 field_type t10 = (*this)[0][2] * (*this)[1][0];
1142 field_type t12 = (*this)[0][1] * (*this)[2][0];
1143 field_type t14 = (*this)[0][2] * (*this)[2][0];
1145 return (t4*(*
this)[2][2]-t6*(*
this)[2][1]-t8*(*
this)[2][2]+
1146 t10*(*
this)[2][1]+t12*(*
this)[1][2]-t14*(*
this)[1][1]);
1150 AutonomousValue<MAT> A(asImp());
1152 Simd::Mask<typename FieldTraits<value_type>::real_type>
1153 nonsingularLanes(
true);
1155 AutonomousValue<MAT>::luDecomposition(A, ElimDet(det), nonsingularLanes,
false, doPivoting);
1156 det =
Simd::cond(nonsingularLanes, det, field_type(0));
1158 for (size_type i = 0; i < rows(); ++i)
1165 namespace DenseMatrixHelp {
1168 template <
typename MAT,
typename V1,
typename V2>
1175 for(size_type i=0; i<matrix.
rows(); ++i)
1178 for(size_type j=0; j<matrix.
cols(); ++j)
1180 ret[i] += matrix[i][j]*x[j];
1187 template <
typename K,
int rows,
int cols>
1190 typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
1192 for(size_type i=0; i<cols(); ++i)
1195 for(size_type j=0; j<rows(); ++j)
1196 ret[i] += matrix[j][i]*x[j];
1201 template <
typename K,
int rows,
int cols>
1202 static inline FieldVector<K,rows> mult(
const FieldMatrix<K,rows,cols> &matrix,
const FieldVector<K,cols> & x)
1204 FieldVector<K,rows> ret;
1210 template <
typename K,
int rows,
int cols>
1211 static inline FieldVector<K,cols> multTransposed(
const FieldMatrix<K,rows,cols> &matrix,
const FieldVector<K,rows> & x)
1213 FieldVector<K,cols> ret;
1214 multAssignTransposed( matrix, x, ret );
1222 template<
typename MAT>
1226 s << a[i] << std::endl;
Macro for wrapping boundary checks.
Generic iterator class for dense vector and matrix implementations.
Definition: densevector.hh:131
A dense n x m matrix.
Definition: densematrix.hh:140
ConstIterator const_iterator
typedef for stl compliant access
Definition: densematrix.hh:244
derived_type operator-() const
Matrix negation.
Definition: densematrix.hh:298
void solve(V1 &x, const V2 &b, bool doPivoting=true) const
Solve system A x = b.
void mv(const X &x, Y &y) const
y = A x
Definition: densematrix.hh:368
Traits::value_type field_type
export the type representing the field
Definition: densematrix.hh:160
derived_type & axpy(const field_type &a, const DenseMatrix< Other > &x)
vector space axpy operation (*this += a x)
Definition: densematrix.hh:338
ConstIterator beforeEnd() const
Definition: densematrix.hh:264
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: densematrix.hh:458
FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: densematrix.hh:561
std::remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: densematrix.hh:248
ConstIterator beforeBegin() const
Definition: densematrix.hh:271
void invert(bool doPivoting=true)
Compute inverse.
static void luDecomposition(DenseMatrix< MAT > &A, Func func, Mask &nonsingularLanes, bool throwEarly, bool doPivoting)
do an LU-Decomposition on matrix A
Traits::value_type block_type
export the type representing the components
Definition: densematrix.hh:163
void mtv(const X &x, Y &y) const
y = A^T x
Definition: densematrix.hh:387
constexpr size_type cols() const
number of columns
Definition: densematrix.hh:715
size_type size() const
size method (number of rows)
Definition: densematrix.hh:200
constexpr size_type M() const
number of columns
Definition: densematrix.hh:703
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:645
Iterator end()
end iterator
Definition: densematrix.hh:222
Iterator beforeBegin()
Definition: densematrix.hh:236
derived_type & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:329
derived_type & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:321
Iterator beforeEnd()
Definition: densematrix.hh:229
Traits::value_type value_type
export the type representing the field
Definition: densematrix.hh:157
void usmv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition: densematrix.hh:484
void usmhv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: densematrix.hh:512
void mmv(const X &x, Y &y) const
y -= A x
Definition: densematrix.hh:445
constexpr size_type rows() const
number of rows
Definition: densematrix.hh:709
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition: densematrix.hh:627
void usmtv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: densematrix.hh:498
derived_type & operator-=(const DenseMatrix< Other > &x)
vector space subtraction
Definition: densematrix.hh:312
bool operator!=(const DenseMatrix< Other > &x) const
Binary matrix incomparison.
Definition: densematrix.hh:358
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: densematrix.hh:471
Traits::derived_type derived_type
type of derived matrix class
Definition: densematrix.hh:154
row_reference operator[](size_type i)
random access
Definition: densematrix.hh:189
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: densematrix.hh:723
Iterator RowIterator
rename the iterators for easier access
Definition: densematrix.hh:211
static constexpr int blocklevel
The number of block levels we contain. This is the leaf, that is, 1.
Definition: densematrix.hh:178
FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: densematrix.hh:528
void umv(const X &x, Y &y) const
y += A x
Definition: densematrix.hh:406
DenseIterator< const DenseMatrix, const row_type, const_row_reference > ConstIterator
Iterator class for sequential access.
Definition: densematrix.hh:242
FieldTraits< vt >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: densematrix.hh:546
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:169
constexpr size_type N() const
number of rows
Definition: densematrix.hh:697
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:166
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:175
FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: densematrix.hh:536
std::remove_reference< row_reference >::type::Iterator ColIterator
rename the iterators for easier access
Definition: densematrix.hh:213
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:172
bool operator==(const DenseMatrix< Other > &x) const
Binary matrix comparison.
Definition: densematrix.hh:348
Iterator iterator
typedef for stl compliant access
Definition: densematrix.hh:209
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: densematrix.hh:246
DenseIterator< DenseMatrix, row_type, row_reference > Iterator
Iterator class for sequential access.
Definition: densematrix.hh:207
void umtv(const X &x, Y &y) const
y += A^T x
Definition: densematrix.hh:419
ConstIterator begin() const
begin iterator
Definition: densematrix.hh:251
field_type determinant(bool doPivoting=true) const
calculates the determinant of this matrix
Iterator begin()
begin iterator
Definition: densematrix.hh:216
void umhv(const X &x, Y &y) const
y += A^H x
Definition: densematrix.hh:432
derived_type & operator+=(const DenseMatrix< Other > &x)
vector space addition
Definition: densematrix.hh:289
ConstIterator end() const
end iterator
Definition: densematrix.hh:257
Interface for a class of dense vectors over a given field.
Definition: densevector.hh:229
size_type size() const
size method
Definition: densevector.hh:336
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:126
static ctype absolute_limit()
return threshold to declare matrix singular
Definition: precision.hh:28
A dense n x m matrix.
Definition: fmatrix.hh:117
vector space out of a tensor product of fields.
Definition: fvector.hh:91
Default exception class for mathematical errors.
Definition: exceptions.hh:241
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:1169
A few common exception classes.
Implements a vector constructed from a given type representing a field and a compile-time given size.
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition: boundschecking.hh:30
typename AutonomousValueType< T >::type AutonomousValue
Type free of internal references that T can be converted to.
Definition: typetraits.hh:588
#define DUNE_THROW(E, m)
Definition: exceptions.hh:218
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
Mask< V > mask(ADLTag< 0, std::is_same< V, Mask< V > >::value >, const V &v)
implements Simd::mask()
Definition: defaults.hh:153
bool anyTrue(const Mask &mask)
Whether any entry is true
Definition: interface.hh:429
V cond(M &&mask, const V &ifTrue, const V &ifFalse)
Like the ?: operator.
Definition: interface.hh:386
bool allTrue(const Mask &mask)
Whether all entries are true
Definition: interface.hh:439
typename Overloads::RebindType< std::decay_t< S >, std::decay_t< V > >::type Rebind
Construct SIMD type with different scalar type.
Definition: interface.hh:253
constexpr std::size_t lanes()
Number of lanes in a SIMD type.
Definition: interface.hh:305
decltype(auto) lane(std::size_t l, V &&v)
Extract an element of a SIMD type.
Definition: interface.hh:324
Rebind< bool, V > Mask
Mask type type of some SIMD type.
Definition: interface.hh:289
Some useful basic math stuff.
bool isNaN(const FieldVector< K, SIZE > &b, PriorityTag< 2 >, ADLTag)
Returns whether any entry is NaN.
Definition: fvector.hh:627
Dune namespace.
Definition: alignedallocator.hh:13
int sign(const T &val)
Return the sign of the value.
Definition: math.hh:180
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:164
Various precision settings for calculations with FieldMatrix and FieldVector.
Implements a scalar vector view wrapper around an existing scalar.
Include file for users of the SIMD abstraction layer.
you have to specialize this structure for any type that should be assignable to a DenseMatrix
Definition: densematrix.hh:59
Traits for type conversions and type information.