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 std::copy_n(l.begin(), std::min(
static_cast<std::size_t
>(ROWS),
148 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
154 using Base::operator=;
168 template <
typename T,
int rows,
int cols>
175 for(
int i = 0; i < ROWS; ++i )
176 for(
int j = 0; j < COLS; ++j )
177 AT[j][i] = (*
this)[i][j];
182 template <
class OtherScalar>
190 result[i][j] = matrixA[i][j] + matrixB[i][j];
196 template <
class OtherScalar>
204 result[i][j] = matrixA[i][j] - matrixB[i][j];
210 template <
class Scalar,
211 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
218 result[i][j] = matrix[i][j] * scalar;
224 template <
class Scalar,
225 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
232 result[i][j] = scalar * matrix[i][j];
238 template <
class Scalar,
239 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
246 result[i][j] = matrix[i][j] / scalar;
253 template <
class OtherScalar,
int otherCols>
264 result[i][j] += matrixA[i][k] * matrixB[k][j];
276 template <
class OtherMatrix, std::enable_if_t<
277 Impl::IsStaticSizeMatrix_v<OtherMatrix>
278 and not Impl::IsFieldMatrix_v<OtherMatrix>
281 const OtherMatrix& matrixB)
285 for (std::size_t j=0; j<
rows; ++j)
286 matrixB.
mtv(matrixA[j], result[j]);
296 template <
class OtherMatrix, std::enable_if_t<
297 Impl::IsStaticSizeMatrix_v<OtherMatrix>
298 and not Impl::IsFieldMatrix_v<OtherMatrix>
305 for (std::size_t j=0; j<
cols; ++j)
307 auto B_j = Impl::ColumnVectorView(matrixB, j);
308 auto result_j = Impl::ColumnVectorView(result, j);
309 matrixA.mv(B_j, result_j);
324 C[i][j] +=
M[i][k]*(*
this)[k][j];
333 template <
int r,
int c>
336 static_assert(r == c,
"Cannot rightmultiply with non-square matrix");
337 static_assert(r ==
cols,
"Size mismatch");
344 (*
this)[i][j] += C[i][k]*
M[k][j];
359 C[i][j] += (*
this)[i][k]*
M[k][j];
386 class FieldMatrix<K,1,1> :
public DenseMatrix< FieldMatrix<K,1,1> >
388 FieldVector<K,1> _data;
389 typedef DenseMatrix< FieldMatrix<K,1,1> > Base;
409 constexpr static int rows = 1;
412 constexpr static int cols = 1;
423 std::copy_n(l.begin(), std::min(
static_cast< std::size_t
>( 1 ), l.size()), &_data);
427 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
433 using Base::operator=;
442 template <
class OtherScalar>
444 const FieldMatrix<OtherScalar,1,1>& matrixB)
446 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] + matrixB[0][0]};
451 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
453 const Scalar& scalar)
455 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] + scalar};
460 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
461 friend auto operator+ (
const Scalar& scalar,
464 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar + matrix[0][0]};
468 template <
class OtherScalar>
470 const FieldMatrix<OtherScalar,1,1>& matrixB)
472 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] - matrixB[0][0]};
477 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
479 const Scalar& scalar)
481 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] - scalar};
486 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
487 friend auto operator- (
const Scalar& scalar,
490 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar - matrix[0][0]};
495 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
498 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] * scalar};
503 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
506 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {scalar * matrix[0][0]};
511 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
514 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] / scalar};
521 template <
class OtherScalar,
int otherCols>
523 const FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
525 FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,otherCols> result;
527 for (
size_type j = 0; j < matrixB.mat_cols(); ++j)
528 result[0][j] = matrixA[0][0] * matrixB[0][j];
539 template <
class OtherMatrix, std::enable_if_t<
540 Impl::IsStaticSizeMatrix_v<OtherMatrix>
541 and not Impl::IsFieldMatrix_v<OtherMatrix>
542 and (OtherMatrix::rows==1)
545 const OtherMatrix& matrixB)
549 for (std::size_t j=0; j<
rows; ++j)
550 matrixB.mtv(matrixA[j], result[j]);
560 template <
class OtherMatrix, std::enable_if_t<
561 Impl::IsStaticSizeMatrix_v<OtherMatrix>
562 and not Impl::IsFieldMatrix_v<OtherMatrix>
563 and (OtherMatrix::cols==1)
565 friend auto operator* (
const OtherMatrix& matrixA,
570 for (std::size_t j=0; j<
cols; ++j)
572 auto B_j = Impl::ColumnVectorView(matrixB, j);
573 auto result_j = Impl::ColumnVectorView(result, j);
574 matrixA.mv(B_j, result_j);
583 FieldMatrix<K,l,1> C;
585 C[j][0] =
M[j][0]*(*
this)[0][0];
600 FieldMatrix<K,1,l> C;
603 C[0][j] =
M[0][j]*_data[0];
653 operator const K& ()
const {
return _data[0]; }
659 std::ostream&
operator<< (std::ostream& s,
const FieldMatrix<K,1,1>& a)
667 namespace FMatrixHelp {
670 template <
typename K>
674 inverse[0][0] = real_type(1.0)/matrix[0][0];
679 template <
typename K>
687 template <
typename K>
692 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
693 K det_1 = real_type(1.0)/det;
694 inverse[0][0] = matrix[1][1] * det_1;
695 inverse[0][1] = - matrix[0][1] * det_1;
696 inverse[1][0] = - matrix[1][0] * det_1;
697 inverse[1][1] = matrix[0][0] * det_1;
703 template <
typename K>
708 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
709 K det_1 = real_type(1.0)/det;
710 inverse[0][0] = matrix[1][1] * det_1;
711 inverse[1][0] = - matrix[0][1] * det_1;
712 inverse[0][1] = - matrix[1][0] * det_1;
713 inverse[1][1] = matrix[0][0] * det_1;
718 template <
typename K>
723 K t4 = matrix[0][0] * matrix[1][1];
724 K t6 = matrix[0][0] * matrix[1][2];
725 K t8 = matrix[0][1] * matrix[1][0];
726 K t10 = matrix[0][2] * matrix[1][0];
727 K t12 = matrix[0][1] * matrix[2][0];
728 K t14 = matrix[0][2] * matrix[2][0];
730 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
731 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
732 K t17 = real_type(1.0)/det;
734 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
735 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
736 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
737 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
738 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
739 inverse[1][2] = -(t6-t10) * t17;
740 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
741 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
742 inverse[2][2] = (t4-t8) * t17;
748 template <
typename K>
753 K t4 = matrix[0][0] * matrix[1][1];
754 K t6 = matrix[0][0] * matrix[1][2];
755 K t8 = matrix[0][1] * matrix[1][0];
756 K t10 = matrix[0][2] * matrix[1][0];
757 K t12 = matrix[0][1] * matrix[2][0];
758 K t14 = matrix[0][2] * matrix[2][0];
760 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
761 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
762 K t17 = real_type(1.0)/det;
764 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
765 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
766 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
767 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
768 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
769 inverse[2][1] = -(t6-t10) * t17;
770 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
771 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
772 inverse[2][2] = (t4-t8) * t17;
778 template<
class K,
int m,
int n,
int p >
785 for( size_type i = 0; i < m; ++i )
787 for( size_type j = 0; j < p; ++j )
789 ret[ i ][ j ] = K( 0 );
790 for( size_type k = 0; k < n; ++k )
791 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
797 template <
typename K,
int rows,
int cols>
802 for(size_type i=0; i<cols; i++)
803 for(size_type j=0; j<cols; j++)
806 for(size_type k=0; k<rows; k++)
807 ret[i][j]+=matrix[k][i]*matrix[k][j];
814 template <
typename K,
int rows,
int cols>
819 for(size_type i=0; i<cols; ++i)
822 for(size_type j=0; j<rows; ++j)
823 ret[i] += matrix[j][i]*x[j];
828 template <
typename K,
int rows,
int cols>
832 multAssign(matrix,x,ret);
837 template <
typename K,
int rows,
int cols>
Various precision settings for calculations with FieldMatrix and FieldVector.
Implements a vector constructed from a given type representing a field and a compile-time given size.
Compute type of the result of an arithmetic operation involving two different number types.
Traits for type conversions and type information.
Eigenvalue computations for the FieldMatrix class.
Implements a matrix constructed from a given type representing a field and a compile-time given numbe...
A few common exception classes.
Macro for wrapping boundary checks.
#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:838
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:680
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:779
static K invertMatrix(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition fmatrix.hh:671
static FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition fmatrix.hh:829
static void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition fmatrix.hh:798
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:815
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:367
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:161
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:351
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:316
FieldMatrix & rightmultiply(const FieldMatrix< K, r, c > &M)
Multiplies M from the right to this matrix.
Definition fmatrix.hh:334
FieldMatrix(T const &rhs)
Definition fmatrix.hh:149
friend auto operator*(const FieldMatrix &matrix, Scalar scalar)
vector space multiplication with scalar
Definition fmatrix.hh:212
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
row_reference mat_access(size_type i)
Definition fmatrix.hh:369
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:172
static constexpr size_type mat_rows()
Definition fmatrix.hh:366
FieldMatrix(std::initializer_list< Dune::FieldVector< K, cols > > const &l)
Constructor initializing the matrix from a list of vector.
Definition fmatrix.hh:140
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:240
friend auto operator+(const FieldMatrix &matrixA, const FieldMatrix< OtherScalar, ROWS, COLS > &matrixB)
vector space addition – two-argument version
Definition fmatrix.hh:183
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:197
const_row_reference mat_access(size_type i) const
Definition fmatrix.hh:375
vector space out of a tensor product of fields.
Definition fvector.hh:95
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