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];
211 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
218 result[i][j] = matrix[i][j] * scalar;
225 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
232 result[i][j] = scalar * matrix[i][j];
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>
455 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] + scalar};
460 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
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>
481 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] - scalar};
486 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
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>
837 template <
typename K,
int rows,
int cols>
Traits for type conversions and type information.
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 vector constructed from a given type representing a field and a compile-time given size.
Eigenvalue computations for the FieldMatrix class.
A few common exception classes.
Implements a matrix constructed from a given type representing a field and a compile-time given numbe...
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
auto min(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:447
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
derived_type operator-() const
Matrix negation.
Definition: densematrix.hh:298
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
const_row_reference mat_access(size_type i) const
Definition: fmatrix.hh:375
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
Definition: ftraits.hh:26
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