12#include <initializer_list>
30 class ColumnVectorView
34 using value_type =
typename M::value_type;
35 using size_type =
typename M::size_type;
37 constexpr ColumnVectorView(M& matrix, size_type col) :
42 constexpr size_type N ()
const {
46 template<
class M_ = M,
47 std::enable_if_t<std::is_same_v<M_,M> and not std::is_const_v<M_>,
int> = 0>
48 constexpr value_type& operator[] (size_type row) {
49 return matrix_[row][col_];
52 constexpr const value_type& operator[] (size_type row)
const {
53 return matrix_[row][col_];
64 struct FieldTraits< Impl::ColumnVectorView<M> >
81 template<
class K,
int ROWS,
int COLS = ROWS >
class FieldMatrix;
84 template<
class K,
int ROWS,
int COLS >
97 typedef typename container_type::size_type
size_type;
100 template<
class K,
int ROWS,
int COLS >
115 template<
class K,
int ROWS,
int COLS>
118 std::array< FieldVector<K,COLS>, ROWS > _data;
123 constexpr static int rows = ROWS;
125 constexpr static int cols = COLS;
141 assert(l.size() ==
rows);
142 for(std::size_t i=0; i<std::min(static_cast<std::size_t>(ROWS), l.size()); ++i)
143 _data[i] = std::data(l)[i];
147 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
153 using Base::operator=;
167 template <
typename T,
int rows,
int cols>
174 for(
int i = 0; i < ROWS; ++i )
175 for(
int j = 0; j < COLS; ++j )
176 AT[j][i] = (*
this)[i][j];
181 template <
class OtherScalar>
189 result[i][j] = matrixA[i][j] + matrixB[i][j];
195 template <
class OtherScalar>
203 result[i][j] = matrixA[i][j] - matrixB[i][j];
209 template <
class Scalar,
210 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
217 result[i][j] = matrix[i][j] * scalar;
223 template <
class Scalar,
224 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
231 result[i][j] = scalar * matrix[i][j];
237 template <
class Scalar,
238 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
245 result[i][j] = matrix[i][j] / scalar;
252 template <
class OtherScalar,
int otherCols>
263 result[i][j] += matrixA[i][k] * matrixB[k][j];
275 template <
class OtherMatrix, std::enable_if_t<
276 Impl::IsStaticSizeMatrix_v<OtherMatrix>
277 and not Impl::IsFieldMatrix_v<OtherMatrix>
280 const OtherMatrix& matrixB)
284 for (std::size_t j=0; j<
rows; ++j)
285 matrixB.
mtv(matrixA[j], result[j]);
295 template <
class OtherMatrix, std::enable_if_t<
296 Impl::IsStaticSizeMatrix_v<OtherMatrix>
297 and not Impl::IsFieldMatrix_v<OtherMatrix>
304 for (std::size_t j=0; j<
cols; ++j)
306 auto B_j = Impl::ColumnVectorView(matrixB, j);
307 auto result_j = Impl::ColumnVectorView(result, j);
308 matrixA.mv(B_j, result_j);
323 C[i][j] +=
M[i][k]*(*
this)[k][j];
332 template <
int r,
int c>
335 static_assert(r == c,
"Cannot rightmultiply with non-square matrix");
336 static_assert(r ==
cols,
"Size mismatch");
343 (*
this)[i][j] += C[i][k]*
M[k][j];
358 C[i][j] += (*
this)[i][k]*
M[k][j];
385 class FieldMatrix<K,1,1> :
public DenseMatrix< FieldMatrix<K,1,1> >
387 FieldVector<K,1> _data;
388 typedef DenseMatrix< FieldMatrix<K,1,1> > Base;
408 constexpr static int rows = 1;
411 constexpr static int cols = 1;
422 std::copy_n(l.begin(), std::min(
static_cast< std::size_t
>( 1 ), l.size()), &_data);
426 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
432 using Base::operator=;
441 template <
class OtherScalar>
443 const FieldMatrix<OtherScalar,1,1>& matrixB)
445 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] + matrixB[0][0]};
450 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
452 const Scalar& scalar)
454 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] + scalar};
459 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
460 friend auto operator+ (
const Scalar& scalar,
463 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar + matrix[0][0]};
467 template <
class OtherScalar>
469 const FieldMatrix<OtherScalar,1,1>& matrixB)
471 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] - matrixB[0][0]};
476 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
478 const Scalar& scalar)
480 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] - scalar};
485 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
486 friend auto operator- (
const Scalar& scalar,
489 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar - matrix[0][0]};
494 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
497 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] * scalar};
502 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
505 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {scalar * matrix[0][0]};
510 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
513 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] / scalar};
520 template <
class OtherScalar,
int otherCols>
522 const FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
524 FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,otherCols> result;
526 for (
size_type j = 0; j < matrixB.mat_cols(); ++j)
527 result[0][j] = matrixA[0][0] * matrixB[0][j];
538 template <
class OtherMatrix, std::enable_if_t<
539 Impl::IsStaticSizeMatrix_v<OtherMatrix>
540 and not Impl::IsFieldMatrix_v<OtherMatrix>
541 and (OtherMatrix::rows==1)
544 const OtherMatrix& matrixB)
548 for (std::size_t j=0; j<
rows; ++j)
549 matrixB.mtv(matrixA[j], result[j]);
559 template <
class OtherMatrix, std::enable_if_t<
560 Impl::IsStaticSizeMatrix_v<OtherMatrix>
561 and not Impl::IsFieldMatrix_v<OtherMatrix>
562 and (OtherMatrix::cols==1)
564 friend auto operator* (
const OtherMatrix& matrixA,
569 for (std::size_t j=0; j<
cols; ++j)
571 auto B_j = Impl::ColumnVectorView(matrixB, j);
572 auto result_j = Impl::ColumnVectorView(result, j);
573 matrixA.mv(B_j, result_j);
582 FieldMatrix<K,l,1> C;
584 C[j][0] =
M[j][0]*(*
this)[0][0];
599 FieldMatrix<K,1,l> C;
602 C[0][j] =
M[0][j]*_data[0];
652 operator const K& ()
const {
return _data[0]; }
658 std::ostream&
operator<< (std::ostream& s,
const FieldMatrix<K,1,1>& a)
666 namespace FMatrixHelp {
669 template <
typename K>
673 inverse[0][0] = real_type(1.0)/matrix[0][0];
678 template <
typename K>
686 template <
typename K>
691 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
692 K det_1 = real_type(1.0)/det;
693 inverse[0][0] = matrix[1][1] * det_1;
694 inverse[0][1] = - matrix[0][1] * det_1;
695 inverse[1][0] = - matrix[1][0] * det_1;
696 inverse[1][1] = matrix[0][0] * det_1;
702 template <
typename K>
707 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
708 K det_1 = real_type(1.0)/det;
709 inverse[0][0] = matrix[1][1] * det_1;
710 inverse[1][0] = - matrix[0][1] * det_1;
711 inverse[0][1] = - matrix[1][0] * det_1;
712 inverse[1][1] = matrix[0][0] * det_1;
717 template <
typename K>
722 K t4 = matrix[0][0] * matrix[1][1];
723 K t6 = matrix[0][0] * matrix[1][2];
724 K t8 = matrix[0][1] * matrix[1][0];
725 K t10 = matrix[0][2] * matrix[1][0];
726 K t12 = matrix[0][1] * matrix[2][0];
727 K t14 = matrix[0][2] * matrix[2][0];
729 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
730 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
731 K t17 = real_type(1.0)/det;
733 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
734 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
735 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
736 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
737 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
738 inverse[1][2] = -(t6-t10) * t17;
739 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
740 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
741 inverse[2][2] = (t4-t8) * t17;
747 template <
typename K>
752 K t4 = matrix[0][0] * matrix[1][1];
753 K t6 = matrix[0][0] * matrix[1][2];
754 K t8 = matrix[0][1] * matrix[1][0];
755 K t10 = matrix[0][2] * matrix[1][0];
756 K t12 = matrix[0][1] * matrix[2][0];
757 K t14 = matrix[0][2] * matrix[2][0];
759 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
760 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
761 K t17 = real_type(1.0)/det;
763 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
764 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
765 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
766 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
767 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
768 inverse[2][1] = -(t6-t10) * t17;
769 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
770 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
771 inverse[2][2] = (t4-t8) * t17;
777 template<
class K,
int m,
int n,
int p >
784 for( size_type i = 0; i < m; ++i )
786 for( size_type j = 0; j < p; ++j )
788 ret[ i ][ j ] = K( 0 );
789 for( size_type k = 0; k < n; ++k )
790 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
796 template <
typename K,
int rows,
int cols>
801 for(size_type i=0; i<cols; i++)
802 for(size_type j=0; j<cols; j++)
805 for(size_type k=0; k<rows; k++)
806 ret[i][j]+=matrix[k][i]*matrix[k][j];
813 template <
typename K,
int rows,
int cols>
818 for(size_type i=0; i<cols; ++i)
821 for(size_type j=0; j<rows; ++j)
822 ret[i] += matrix[j][i]*x[j];
827 template <
typename K,
int rows,
int cols>
831 multAssign(matrix,x,ret);
836 template <
typename K,
int rows,
int cols>
A few common exception classes.
Compute type of the result of an arithmetic operation involving two different number types.
Various precision settings for calculations with FieldMatrix and FieldVector.
Implements a matrix constructed from a given type representing a field and a compile-time given numbe...
Macro for wrapping boundary checks.
Eigenvalue computations for the FieldMatrix class.
Implements a vector constructed from a given type representing a field and a compile-time given size.
Traits for type conversions and type information.
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition boundschecking.hh:30
std::ostream & operator<<(std::ostream &s, const bigunsignedint< k > &x)
Definition bigunsignedint.hh:278
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition simd/interface.hh:235
Dune namespace.
Definition alignedallocator.hh:13
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition densematrix.hh:1169
static FieldVector< K, cols > multTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x)
calculates ret = matrix^T * x
Definition fmatrix.hh:837
static K invertMatrix_retTransposed(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition fmatrix.hh:679
static void multMatrix(const FieldMatrix< K, m, n > &A, const FieldMatrix< K, n, p > &B, FieldMatrix< K, m, p > &ret)
calculates ret = A * B
Definition fmatrix.hh:778
static K invertMatrix(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition fmatrix.hh:670
static FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition fmatrix.hh:828
static void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition fmatrix.hh:797
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:814
A dense n x m matrix.
Definition densematrix.hh:140
void mtv(const X &x, Y &y) const
y = A^T x
Definition densematrix.hh:387
constexpr size_type M() const
number of columns
Definition densematrix.hh:703
FieldMatrix< K, ROWS, COLS > & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition densematrix.hh:645
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
derived_type & operator-=(const DenseMatrix< Other > &x)
vector space subtraction
Definition densematrix.hh:312
static constexpr int blocklevel
The number of block levels we contain. This is the leaf, that is, 1.
Definition densematrix.hh:178
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition densematrix.hh:169
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
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition densematrix.hh:172
derived_type & operator+=(const DenseMatrix< Other > &x)
vector space addition
Definition densematrix.hh:289
A dense n x m matrix.
Definition fmatrix.hh:117
static constexpr size_type mat_cols()
Definition fmatrix.hh:366
constexpr FieldMatrix()=default
Default constructor.
Base::const_row_reference const_row_reference
Definition fmatrix.hh:131
FieldMatrix & operator=(const FieldMatrix< T, ROWS, COLS > &x)
copy assignment from FieldMatrix over a different field
Definition fmatrix.hh:160
FieldMatrix< K, rows, l > rightmultiplyany(const FieldMatrix< K, cols, l > &M) const
Multiplies M from the right to this matrix, this matrix is not modified.
Definition fmatrix.hh:350
Base::row_type row_type
Definition fmatrix.hh:128
Base::size_type size_type
Definition fmatrix.hh:127
FieldMatrix< K, l, cols > leftmultiplyany(const FieldMatrix< K, l, rows > &M) const
Multiplies M from the left to this matrix, this matrix is not modified.
Definition fmatrix.hh:315
FieldMatrix & rightmultiply(const FieldMatrix< K, r, c > &M)
Multiplies M from the right to this matrix.
Definition fmatrix.hh:333
FieldMatrix(T const &rhs)
Definition fmatrix.hh:148
friend auto operator*(const FieldMatrix &matrix, Scalar scalar)
vector space multiplication with scalar
Definition fmatrix.hh:211
FieldMatrix & operator=(FieldMatrix< T, rows, cols > const &)=delete
no copy assignment from FieldMatrix of different size
Base::row_reference row_reference
Definition fmatrix.hh:130
constexpr FieldMatrix(std::initializer_list< Dune::FieldVector< K, cols > > const &l)
Constructor initializing the matrix from a list of vector.
Definition fmatrix.hh:140
row_reference mat_access(size_type i)
Definition fmatrix.hh:368
static constexpr int rows
The number of rows.
Definition fmatrix.hh:123
FieldMatrix< K, COLS, ROWS > transposed() const
Return transposed of the matrix as FieldMatrix.
Definition fmatrix.hh:171
static constexpr size_type mat_rows()
Definition fmatrix.hh:365
static constexpr int cols
The number of columns.
Definition fmatrix.hh:125
friend auto operator/(const FieldMatrix &matrix, Scalar scalar)
vector space division by scalar
Definition fmatrix.hh:239
friend auto operator+(const FieldMatrix &matrixA, const FieldMatrix< OtherScalar, ROWS, COLS > &matrixB)
vector space addition – two-argument version
Definition fmatrix.hh:182
FieldMatrix & operator=(const FieldMatrix &)=default
copy assignment operator
friend auto operator-(const FieldMatrix &matrixA, const FieldMatrix< OtherScalar, ROWS, COLS > &matrixB)
vector space subtraction – two-argument version
Definition fmatrix.hh:196
const_row_reference mat_access(size_type i) const
Definition fmatrix.hh:374
vector space out of a tensor product of fields.
Definition fvector.hh:91
std::array< row_type, ROWS > container_type
Definition fmatrix.hh:95
FieldMatrix< K, ROWS, COLS > derived_type
Definition fmatrix.hh:87
K value_type
Definition fmatrix.hh:96
const row_type & const_row_reference
Definition fmatrix.hh:93
FieldVector< K, COLS > row_type
Definition fmatrix.hh:90
container_type::size_type size_type
Definition fmatrix.hh:97
row_type & row_reference
Definition fmatrix.hh:92
FieldTraits< K >::field_type field_type
Definition fmatrix.hh:103
FieldTraits< K >::real_type real_type
Definition fmatrix.hh:104
T field_type
export the type representing the field
Definition ftraits.hh:28
T real_type
export the type representing the real type of the field
Definition ftraits.hh:30
Definition matvectraits.hh:31
decltype(std::declval< T1 >()+std::declval< T2 >()) PromotedType
Definition promotiontraits.hh:28