3 #ifndef DUNE_FMATRIX_HH
4 #define DUNE_FMATRIX_HH
10 #include <initializer_list>
34 template<
class K,
int ROWS,
int COLS = ROWS >
class FieldMatrix;
36 template<
class K,
int ROWS,
int COLS >
49 typedef typename container_type::size_type
size_type;
52 template<
class K,
int ROWS,
int COLS >
67 template<
class K,
int ROWS,
int COLS>
70 std::array< FieldVector<K,COLS>, ROWS > _data;
96 assert(l.size() ==
rows);
97 std::copy_n(l.begin(),
std::min(
static_cast<std::size_t
>(ROWS),
103 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
109 using Base::operator=;
123 template <
typename T,
int rows,
int cols>
127 template <
class OtherScalar>
135 result[i][j] = matrixA[i][j] + matrixB[i][j];
141 template <
class OtherScalar>
149 result[i][j] = matrixA[i][j] - matrixB[i][j];
156 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
163 result[i][j] = matrix[i][j] * scalar;
170 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
177 result[i][j] = scalar * matrix[i][j];
184 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
191 result[i][j] = matrix[i][j] / scalar;
198 template <
class OtherScalar,
int otherCols>
209 result[i][j] += matrixA[i][k] * matrixB[k][j];
225 C[i][j] +=
M[i][k]*(*
this)[k][j];
234 template <
int r,
int c>
237 static_assert(r == c,
"Cannot rightmultiply with non-square matrix");
238 static_assert(r ==
cols,
"Size mismatch");
245 (*
this)[i][j] += C[i][k]*
M[k][j];
260 C[i][j] += (*
this)[i][k]*
M[k][j];
287 class FieldMatrix<K,1,1> :
public DenseMatrix< FieldMatrix<K,1,1> >
289 FieldVector<K,1> _data;
290 typedef DenseMatrix< FieldMatrix<K,1,1> > Base;
330 std::copy_n(l.begin(),
std::min(
static_cast< std::size_t
>( 1 ), l.size()), &_data);
334 typename = std::enable_if_t<HasDenseMatrixAssigner<FieldMatrix, T>::value>>
340 using Base::operator=;
343 template <
class OtherScalar>
345 const FieldMatrix<OtherScalar,1,1>& matrixB)
347 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] + matrixB[0][0]};
352 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
356 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] + scalar};
361 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
365 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar + matrix[0][0]};
369 template <
class OtherScalar>
371 const FieldMatrix<OtherScalar,1,1>& matrixB)
373 return FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,1>{matrixA[0][0] - matrixB[0][0]};
378 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
382 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1>{matrix[0][0] - scalar};
387 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
391 return FieldMatrix<typename PromotionTraits<Scalar,K>::PromotedType,1,1>{scalar - matrix[0][0]};
396 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
399 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] * scalar};
404 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
407 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {scalar * matrix[0][0]};
412 std::enable_if_t<IsNumber<Scalar>::value,
int> = 0>
415 return FieldMatrix<typename PromotionTraits<K,Scalar>::PromotedType,1,1> {matrix[0][0] / scalar};
422 template <
class OtherScalar,
int otherCols>
424 const FieldMatrix<OtherScalar, 1, otherCols>& matrixB)
426 FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,1,otherCols> result;
428 for (
size_type j = 0; j < matrixB.mat_cols(); ++j)
429 result[0][j] = matrixA[0][0] * matrixB[0][j];
438 FieldMatrix<K,l,1> C;
440 C[j][0] =
M[j][0]*(*
this)[0][0];
455 FieldMatrix<K,1,l> C;
458 C[0][j] =
M[0][j]*_data[0];
508 operator const K& ()
const {
return _data[0]; }
514 std::ostream&
operator<< (std::ostream& s,
const FieldMatrix<K,1,1>& a)
522 namespace FMatrixHelp {
525 template <
typename K>
529 inverse[0][0] = real_type(1.0)/matrix[0][0];
534 template <
typename K>
542 template <
typename K>
547 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
548 K det_1 = real_type(1.0)/det;
549 inverse[0][0] = matrix[1][1] * det_1;
550 inverse[0][1] = - matrix[0][1] * det_1;
551 inverse[1][0] = - matrix[1][0] * det_1;
552 inverse[1][1] = matrix[0][0] * det_1;
558 template <
typename K>
563 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
564 K det_1 = real_type(1.0)/det;
565 inverse[0][0] = matrix[1][1] * det_1;
566 inverse[1][0] = - matrix[0][1] * det_1;
567 inverse[0][1] = - matrix[1][0] * det_1;
568 inverse[1][1] = matrix[0][0] * det_1;
573 template <
typename K>
578 K t4 = matrix[0][0] * matrix[1][1];
579 K t6 = matrix[0][0] * matrix[1][2];
580 K t8 = matrix[0][1] * matrix[1][0];
581 K t10 = matrix[0][2] * matrix[1][0];
582 K t12 = matrix[0][1] * matrix[2][0];
583 K t14 = matrix[0][2] * matrix[2][0];
585 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
586 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
587 K t17 = real_type(1.0)/det;
589 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
590 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
591 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
592 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
593 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
594 inverse[1][2] = -(t6-t10) * t17;
595 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
596 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
597 inverse[2][2] = (t4-t8) * t17;
603 template <
typename K>
608 K t4 = matrix[0][0] * matrix[1][1];
609 K t6 = matrix[0][0] * matrix[1][2];
610 K t8 = matrix[0][1] * matrix[1][0];
611 K t10 = matrix[0][2] * matrix[1][0];
612 K t12 = matrix[0][1] * matrix[2][0];
613 K t14 = matrix[0][2] * matrix[2][0];
615 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
616 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
617 K t17 = real_type(1.0)/det;
619 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
620 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
621 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
622 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
623 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
624 inverse[2][1] = -(t6-t10) * t17;
625 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
626 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
627 inverse[2][2] = (t4-t8) * t17;
633 template<
class K,
int m,
int n,
int p >
640 for( size_type i = 0; i < m; ++i )
642 for( size_type j = 0; j < p; ++j )
644 ret[ i ][ j ] = K( 0 );
645 for( size_type k = 0; k < n; ++k )
646 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
652 template <
typename K,
int rows,
int cols>
657 for(size_type i=0; i<cols; i++)
658 for(size_type j=0; j<cols; j++)
661 for(size_type k=0; k<rows; k++)
662 ret[i][j]+=matrix[k][i]*matrix[k][j];
669 template <
typename K,
int rows,
int cols>
674 for(size_type i=0; i<cols; ++i)
677 for(size_type j=0; j<rows; ++j)
678 ret[i] += matrix[j][i]*x[j];
683 template <
typename K,
int rows,
int cols>
692 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:28
std::ostream & operator<<(std::ostream &s, const bigunsignedint< k > &x)
Definition: bigunsignedint.hh:273
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: simd/interface.hh:233
Dune namespace.
Definition: alignedallocator.hh:11
auto min(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:434
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition: densematrix.hh:1197
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:535
static FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition: fmatrix.hh:684
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:634
static K invertMatrix(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition: fmatrix.hh:526
static FieldVector< K, cols > multTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x)
calculates ret = matrix^T * x
Definition: fmatrix.hh:693
static void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition: fmatrix.hh:653
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:670
A dense n x m matrix.
Definition: densematrix.hh:165
derived_type operator-() const
Matrix negation.
Definition: densematrix.hh:326
FieldMatrix< K, ROWS, COLS > & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:673
constexpr size_type M() const
number of columns
Definition: densematrix.hh:731
derived_type & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:357
derived_type & operator-=(const DenseMatrix< Other > &x)
vector space subtraction
Definition: densematrix.hh:340
derived_type & operator+=(const DenseMatrix< Other > &x)
vector space addition
Definition: densematrix.hh:317
derived_type & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:349
@ blocklevel
The number of block levels we contain. This is 1.
Definition: densematrix.hh:205
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:194
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:191
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:200
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:197
A dense n x m matrix.
Definition: fmatrix.hh:69
static constexpr size_type mat_cols()
Definition: fmatrix.hh:268
constexpr FieldMatrix()=default
Default constructor.
Base::const_row_reference const_row_reference
Definition: fmatrix.hh:86
Base::row_type row_type
Definition: fmatrix.hh:83
Base::size_type size_type
Definition: fmatrix.hh:82
FieldMatrix & operator=(const FieldMatrix< T, ROWS, COLS > &x)
copy assignment from FieldMatrix over a different field
Definition: fmatrix.hh:116
@ rows
The number of rows.
Definition: fmatrix.hh:77
@ cols
The number of columns.
Definition: fmatrix.hh:79
FieldMatrix(T const &rhs)
Definition: fmatrix.hh:104
FieldMatrix & operator=(FieldMatrix< T, rows, cols > const &)=delete
no copy assignment from FieldMatrix of different size
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:217
friend auto operator*(const FieldMatrix &matrix, Scalar scalar)
vector space multiplication with scalar
Definition: fmatrix.hh:157
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:252
Base::row_reference row_reference
Definition: fmatrix.hh:85
FieldMatrix & operator=(const FieldMatrix &)=default
copy assignment operator
row_reference mat_access(size_type i)
Definition: fmatrix.hh:270
static constexpr size_type mat_rows()
Definition: fmatrix.hh:267
FieldMatrix(std::initializer_list< Dune::FieldVector< K, cols > > const &l)
Constructor initializing the matrix from a list of vector.
Definition: fmatrix.hh:95
friend auto operator/(const FieldMatrix &matrix, Scalar scalar)
vector space division by scalar
Definition: fmatrix.hh:185
FieldMatrix & rightmultiply(const FieldMatrix< K, r, c > &M)
Multiplies M from the right to this matrix.
Definition: fmatrix.hh:235
friend auto operator+(const FieldMatrix &matrixA, const FieldMatrix< OtherScalar, ROWS, COLS > &matrixB)
vector space addition – two-argument version
Definition: fmatrix.hh:128
const_row_reference mat_access(size_type i) const
Definition: fmatrix.hh:276
std::array< row_type, ROWS > container_type
Definition: fmatrix.hh:47
FieldMatrix< K, ROWS, COLS > derived_type
Definition: fmatrix.hh:39
K value_type
Definition: fmatrix.hh:48
const row_type & const_row_reference
Definition: fmatrix.hh:45
FieldVector< K, COLS > row_type
Definition: fmatrix.hh:42
container_type::size_type size_type
Definition: fmatrix.hh:49
row_type & row_reference
Definition: fmatrix.hh:44
FieldTraits< K >::field_type field_type
Definition: fmatrix.hh:55
FieldTraits< K >::real_type real_type
Definition: fmatrix.hh:56
Definition: ftraits.hh:24
T field_type
export the type representing the field
Definition: ftraits.hh:26
T real_type
export the type representing the real type of the field
Definition: ftraits.hh:28
Definition: matvectraits.hh:29