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>
Macro for wrapping boundary checks.
Implements a matrix constructed from a given type representing a field and a compile-time given numbe...
A few common exception classes.
Eigenvalue computations for the FieldMatrix class.
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.
Compute type of the result of an arithmetic operation involving two different number types.
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