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];
510 operator const K& ()
const {
return _data[0]; }
516 std::ostream&
operator<< (std::ostream& s,
const FieldMatrix<K,1,1>& a)
524 namespace FMatrixHelp {
527 template <
typename K>
531 inverse[0][0] = real_type(1.0)/matrix[0][0];
536 template <
typename K>
544 template <
typename K>
549 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
550 K det_1 = real_type(1.0)/det;
551 inverse[0][0] = matrix[1][1] * det_1;
552 inverse[0][1] = - matrix[0][1] * det_1;
553 inverse[1][0] = - matrix[1][0] * det_1;
554 inverse[1][1] = matrix[0][0] * det_1;
560 template <
typename K>
565 K det = (matrix[0][0]*matrix[1][1] - matrix[0][1]*matrix[1][0]);
566 K det_1 = real_type(1.0)/det;
567 inverse[0][0] = matrix[1][1] * det_1;
568 inverse[1][0] = - matrix[0][1] * det_1;
569 inverse[0][1] = - matrix[1][0] * det_1;
570 inverse[1][1] = matrix[0][0] * det_1;
575 template <
typename K>
580 K t4 = matrix[0][0] * matrix[1][1];
581 K t6 = matrix[0][0] * matrix[1][2];
582 K t8 = matrix[0][1] * matrix[1][0];
583 K t10 = matrix[0][2] * matrix[1][0];
584 K t12 = matrix[0][1] * matrix[2][0];
585 K t14 = matrix[0][2] * matrix[2][0];
587 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
588 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
589 K t17 = real_type(1.0)/det;
591 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
592 inverse[0][1] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
593 inverse[0][2] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
594 inverse[1][0] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
595 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
596 inverse[1][2] = -(t6-t10) * t17;
597 inverse[2][0] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
598 inverse[2][1] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
599 inverse[2][2] = (t4-t8) * t17;
605 template <
typename K>
610 K t4 = matrix[0][0] * matrix[1][1];
611 K t6 = matrix[0][0] * matrix[1][2];
612 K t8 = matrix[0][1] * matrix[1][0];
613 K t10 = matrix[0][2] * matrix[1][0];
614 K t12 = matrix[0][1] * matrix[2][0];
615 K t14 = matrix[0][2] * matrix[2][0];
617 K det = (t4*matrix[2][2]-t6*matrix[2][1]-t8*matrix[2][2]+
618 t10*matrix[2][1]+t12*matrix[1][2]-t14*matrix[1][1]);
619 K t17 = real_type(1.0)/det;
621 inverse[0][0] = (matrix[1][1] * matrix[2][2] - matrix[1][2] * matrix[2][1])*t17;
622 inverse[1][0] = -(matrix[0][1] * matrix[2][2] - matrix[0][2] * matrix[2][1])*t17;
623 inverse[2][0] = (matrix[0][1] * matrix[1][2] - matrix[0][2] * matrix[1][1])*t17;
624 inverse[0][1] = -(matrix[1][0] * matrix[2][2] - matrix[1][2] * matrix[2][0])*t17;
625 inverse[1][1] = (matrix[0][0] * matrix[2][2] - t14) * t17;
626 inverse[2][1] = -(t6-t10) * t17;
627 inverse[0][2] = (matrix[1][0] * matrix[2][1] - matrix[1][1] * matrix[2][0]) * t17;
628 inverse[1][2] = -(matrix[0][0] * matrix[2][1] - t12) * t17;
629 inverse[2][2] = (t4-t8) * t17;
635 template<
class K,
int m,
int n,
int p >
642 for( size_type i = 0; i < m; ++i )
644 for( size_type j = 0; j < p; ++j )
646 ret[ i ][ j ] = K( 0 );
647 for( size_type k = 0; k < n; ++k )
648 ret[ i ][ j ] += A[ i ][ k ] * B[ k ][ j ];
654 template <
typename K,
int rows,
int cols>
659 for(size_type i=0; i<cols; i++)
660 for(size_type j=0; j<cols; j++)
663 for(size_type k=0; k<rows; k++)
664 ret[i][j]+=matrix[k][i]*matrix[k][j];
671 template <
typename K,
int rows,
int cols>
676 for(size_type i=0; i<cols; ++i)
679 for(size_type j=0; j<rows; ++j)
680 ret[i] += matrix[j][i]*x[j];
685 template <
typename K,
int rows,
int cols>
694 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
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
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:14
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:1198
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:537
static FieldVector< K, rows > mult(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, cols > &x)
calculates ret = matrix * x
Definition: fmatrix.hh:686
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:636
static K invertMatrix(const FieldMatrix< K, 1, 1 > &matrix, FieldMatrix< K, 1, 1 > &inverse)
invert scalar without changing the original matrix
Definition: fmatrix.hh:528
static FieldVector< K, cols > multTransposed(const FieldMatrix< K, rows, cols > &matrix, const FieldVector< K, rows > &x)
calculates ret = matrix^T * x
Definition: fmatrix.hh:695
static void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition: fmatrix.hh:655
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:672
A dense n x m matrix.
Definition: densematrix.hh:164
derived_type operator-() const
Matrix negation.
Definition: densematrix.hh:325
FieldMatrix< K, ROWS, COLS > & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:672
size_type M() const
number of columns
Definition: densematrix.hh:730
derived_type & operator/=(const field_type &k)
vector space division by scalar
Definition: densematrix.hh:356
derived_type & operator-=(const DenseMatrix< Other > &x)
vector space subtraction
Definition: densematrix.hh:339
derived_type & operator+=(const DenseMatrix< Other > &x)
vector space addition
Definition: densematrix.hh:316
derived_type & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: densematrix.hh:348
@ blocklevel
The number of block levels we contain. This is 1.
Definition: densematrix.hh:204
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:193
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:190
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:199
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:196
A dense n x m matrix.
Definition: fmatrix.hh:69
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
constexpr size_type mat_cols() const
Definition: fmatrix.hh:268
row_reference mat_access(size_type i)
Definition: fmatrix.hh:270
FieldMatrix(std::initializer_list< Dune::FieldVector< K, cols > > const &l)
Constructor initializing the matrix from a list of vector.
Definition: fmatrix.hh:95
constexpr size_type mat_rows() const
Definition: fmatrix.hh:267
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