Go to the documentation of this file.
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];
283 #ifndef DOXYGEN // hide specialization
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>
constexpr size_type mat_cols() const
Definition: fmatrix.hh:268
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
A few common exception classes.
K value_type
Definition: fmatrix.hh:48
T real_type
export the type representing the real type of the field
Definition: ftraits.hh:28
FieldMatrix(std::initializer_list< Dune::FieldVector< K, cols > > const &l)
Constructor initializing the matrix from a list of vector.
Definition: fmatrix.hh:95
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
container_type::size_type size_type
Definition: fmatrix.hh:49
const typedef row_type & const_row_reference
Definition: fmatrix.hh:45
Compute type of the result of an arithmetic operation involving two different number types.
auto min(const AlignedNumber< T, align > &a, const AlignedNumber< T, align > &b)
Definition: debugalign.hh:434
FieldTraits< K >::real_type real_type
Definition: fmatrix.hh:56
row_reference mat_access(size_type i)
Definition: fmatrix.hh:270
Base::size_type size_type
Definition: fmatrix.hh:82
Macro for wrapping boundary checks.
#define DUNE_UNUSED_PARAMETER(parm)
A macro to mark intentionally unused function parameters with.
Definition: unused.hh:25
const_row_reference mat_access(size_type i) const
Definition: fmatrix.hh:276
Eigenvalue computations for the FieldMatrix class.
FieldTraits< K >::field_type field_type
Definition: fmatrix.hh:55
static void multTransposedMatrix(const FieldMatrix< K, rows, cols > &matrix, FieldMatrix< K, cols, cols > &ret)
calculates ret= A_t*A
Definition: fmatrix.hh:655
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition: densematrix.hh:672
@ cols
The number of columns.
Definition: fmatrix.hh:79
A dense n x m matrix.
Definition: densematrix.hh:41
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
T field_type
export the type representing the field
Definition: ftraits.hh:26
Implements a vector constructed from a given type representing a field and a compile-time given size.
A dense n x m matrix.
Definition: densematrix.hh:28
FieldMatrix & rightmultiply(const FieldMatrix< K, r, c > &M)
Multiplies M from the right to this matrix.
Definition: fmatrix.hh:235
Traits for type conversions and type information.
FieldMatrix(T const &rhs)
Definition: fmatrix.hh:104
friend auto operator+(const FieldMatrix &matrixA, const FieldMatrix< OtherScalar, ROWS, COLS > &matrixB)
vector space addition – two-argument version
Definition: fmatrix.hh:128
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition: densematrix.hh:1198
@ rows
The number of rows.
Definition: fmatrix.hh:77
FieldMatrix & operator=(const FieldMatrix &)=default
copy assignment operator
std::ostream & operator<<(std::ostream &s, const bigunsignedint< k > &x)
Definition: bigunsignedint.hh:273
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface)
Definition: densematrix.hh:193
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
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
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &)
Definition: densematrix.hh:196
Definition: ftraits.hh:23
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 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
Traits::size_type size_type
The type used for the index access and size operation.
Definition: densematrix.hh:190
Base::row_type row_type
Definition: fmatrix.hh:83
row_type & row_reference
Definition: fmatrix.hh:44
friend auto operator/(const FieldMatrix &matrix, Scalar scalar)
vector space division by scalar
Definition: fmatrix.hh:185
friend auto operator*(const FieldMatrix &matrix, Scalar scalar)
vector space multiplication with scalar
Definition: fmatrix.hh:157
Definition: matvectraits.hh:29
Various precision settings for calculations with FieldMatrix and FieldVector.
constexpr FieldMatrix()=default
Default constructor.
FieldMatrix< K, ROWS, COLS > derived_type
Definition: fmatrix.hh:39
Base::row_reference row_reference
Definition: fmatrix.hh:85
std::array< row_type, ROWS > container_type
Definition: fmatrix.hh:47
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
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition: boundschecking.hh:28
constexpr size_type mat_rows() const
Definition: fmatrix.hh:267
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
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: simd/interface.hh:233
Implements a matrix constructed from a given type representing a field and a compile-time given numbe...
friend auto operator-(const FieldMatrix &matrixA, const FieldMatrix< OtherScalar, ROWS, COLS > &matrixB)
vector space subtraction – two-argument version
Definition: fmatrix.hh:142
FieldVector< K, COLS > row_type
Definition: fmatrix.hh:42
FieldMatrix & operator=(const FieldMatrix< T, ROWS, COLS > &x)
copy assignment from FieldMatrix over a different field
Definition: fmatrix.hh:116
Dune namespace.
Definition: alignedallocator.hh:13
Base::const_row_reference const_row_reference
Definition: fmatrix.hh:86