Loading...
Searching...
No Matches
densematrix.hh
Go to the documentation of this file.
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
// SPDX-FileCopyrightInfo: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
#ifndef DUNE_DENSEMATRIX_HH
#define DUNE_DENSEMATRIX_HH
#include <cmath>
#include <cstddef>
#include <iostream>
#include <iterator>
#include <type_traits>
#include <utility>
#include <vector>
#include <dune/common/boundschecking.hh>
#include <dune/common/classname.hh>
#include <dune/common/exceptions.hh>
#include <dune/common/fvector.hh>
#include <dune/common/math.hh>
#include <dune/common/precision.hh>
#include <dune/common/simd/simd.hh>
#include <dune/common/typetraits.hh>
#include <dune/common/scalarvectorview.hh>
#include <dune/common/std/iterator.hh>
namespace Dune
{
template<typename M> class DenseMatrix;
template<typename M>
struct FieldTraits< DenseMatrix<M> >
{
typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::field_type field_type;
typedef const typename FieldTraits< typename DenseMatVecTraits<M>::value_type >::real_type real_type;
};
template<class K, int N, int M> class FieldMatrix;
template<class K, int N> class FieldVector;
template< class DenseMatrix, class RHS >
struct DenseMatrixAssigner;
#ifndef DOXYGEN
namespace Impl
{
template< class DenseMatrix, class RHS >
class DenseMatrixAssigner
{};
template< class DenseMatrix, class RHS >
requires (Dune::IsNumber< RHS >::value)
class DenseMatrixAssigner<DenseMatrix, RHS>
{
public:
constexpr static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
{
typedef typename DenseMatrix::field_type field_type;
std::fill( denseMatrix.begin(), denseMatrix.end(), static_cast< field_type >( rhs ) );
}
};
template< class DenseMatrix, class RHS >
requires Std::indirectly_copyable<
decltype(std::begin(*std::declval<typename RHS::const_iterator>())),
decltype(std::begin(*std::declval<typename DenseMatrix::iterator>()))>
class DenseMatrixAssigner<DenseMatrix, RHS>
{
public:
constexpr static void apply ( DenseMatrix &denseMatrix, const RHS &rhs )
{
DUNE_ASSERT_BOUNDS(rhs.N() == denseMatrix.N());
DUNE_ASSERT_BOUNDS(rhs.M() == denseMatrix.M());
typename DenseMatrix::iterator tIt = std::begin(denseMatrix);
typename RHS::const_iterator sIt = std::begin(rhs);
for(; sIt != std::end(rhs); ++tIt, ++sIt)
std::copy(std::begin(*sIt), std::end(*sIt), std::begin(*tIt));
}
};
} // namespace Impl
template< class DenseMatrix, class RHS >
struct DenseMatrixAssigner
: public Impl::DenseMatrixAssigner< DenseMatrix, RHS >
{};
namespace Impl
{
template< class DenseMatrix, class RHS >
std::true_type hasDenseMatrixAssigner ( DenseMatrix &, const RHS &, decltype( Dune::DenseMatrixAssigner< DenseMatrix, RHS >::apply( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) ) * = nullptr );
std::false_type hasDenseMatrixAssigner ( ... );
} // namespace Impl
template< class DenseMatrix, class RHS >
struct HasDenseMatrixAssigner
: public decltype( Impl::hasDenseMatrixAssigner( std::declval< DenseMatrix & >(), std::declval< const RHS & >() ) )
{};
#endif // #ifndef DOXYGEN
class FMatrixError : public MathError {};
template<typename MAT>
class DenseMatrix
{
typedef DenseMatVecTraits<MAT> Traits;
// Curiously recurring template pattern
constexpr MAT & asImp() { return static_cast<MAT&>(*this); }
constexpr const MAT & asImp() const { return static_cast<const MAT&>(*this); }
template <class>
friend class DenseMatrix;
public:
//===== type definitions and constants
typedef typename Traits::derived_type derived_type;
typedef typename Traits::value_type value_type;
typedef typename Traits::value_type field_type;
typedef typename Traits::value_type block_type;
typedef typename Traits::size_type size_type;
typedef typename Traits::row_type row_type;
typedef typename Traits::row_reference row_reference;
typedef typename Traits::const_row_reference const_row_reference;
constexpr static int blocklevel = 1;
private:
using simd_index_type = Simd::Rebind<std::size_t, value_type>;
public:
//===== access to components
constexpr row_reference operator[] ( size_type i )
{
return asImp().mat_access(i);
}
constexpr const_row_reference operator[] ( size_type i ) const
{
return asImp().mat_access(i);
}
constexpr size_type size() const
{
return rows();
}
//===== iterator interface to rows of the matrix
typedef DenseIterator<DenseMatrix,row_type,row_reference> Iterator;
typedef Iterator iterator;
typedef Iterator RowIterator;
typedef typename std::remove_reference<row_reference>::type::Iterator ColIterator;
constexpr Iterator begin ()
{
return Iterator(*this,0);
}
constexpr Iterator end ()
{
return Iterator(*this,rows());
}
constexpr Iterator beforeEnd ()
{
return Iterator(*this,rows()-1);
}
constexpr Iterator beforeBegin ()
{
return Iterator(*this,-1);
}
typedef DenseIterator<const DenseMatrix,const row_type,const_row_reference> ConstIterator;
typedef ConstIterator const_iterator;
typedef ConstIterator ConstRowIterator;
typedef typename std::remove_reference<const_row_reference>::type::ConstIterator ConstColIterator;
constexpr ConstIterator begin () const
{
return ConstIterator(*this,0);
}
constexpr ConstIterator end () const
{
return ConstIterator(*this,rows());
}
constexpr ConstIterator beforeEnd () const
{
return ConstIterator(*this,rows()-1);
}
constexpr ConstIterator beforeBegin () const
{
return ConstIterator(*this,-1);
}
//===== assignment
template< class RHS, class = std::enable_if_t< HasDenseMatrixAssigner< MAT, RHS >::value > >
constexpr derived_type &operator= ( const RHS &rhs )
{
DenseMatrixAssigner< MAT, RHS >::apply( asImp(), rhs );
return asImp();
}
//===== vector space arithmetic
template <class Other>
constexpr derived_type &operator+= (const DenseMatrix<Other>& x)
{
DUNE_ASSERT_BOUNDS(rows() == x.rows());
for (size_type i=0; i<rows(); i++)
(*this)[i] += x[i];
return asImp();
}
constexpr derived_type operator- () const
{
MAT result;
using idx_type = typename decltype(result)::size_type;
for (idx_type i = 0; i < rows(); ++i)
for (idx_type j = 0; j < cols(); ++j)
result[i][j] = - asImp()[i][j];
return result;
}
template <class Other>
constexpr derived_type &operator-= (const DenseMatrix<Other>& x)
{
DUNE_ASSERT_BOUNDS(rows() == x.rows());
for (size_type i=0; i<rows(); i++)
(*this)[i] -= x[i];
return asImp();
}
constexpr derived_type &operator*= (const field_type& k)
{
for (size_type i=0; i<rows(); i++)
(*this)[i] *= k;
return asImp();
}
constexpr derived_type &operator/= (const field_type& k)
{
for (size_type i=0; i<rows(); i++)
(*this)[i] /= k;
return asImp();
}
template <class Other>
constexpr derived_type &axpy (const field_type &a, const DenseMatrix<Other> &x )
{
DUNE_ASSERT_BOUNDS(rows() == x.rows());
for( size_type i = 0; i < rows(); ++i )
(*this)[ i ].axpy( a, x[ i ] );
return asImp();
}
template <class Other>
constexpr bool operator== (const DenseMatrix<Other>& x) const
{
DUNE_ASSERT_BOUNDS(rows() == x.rows());
for (size_type i=0; i<rows(); i++)
if ((*this)[i]!=x[i])
return false;
return true;
}
template <class Other>
constexpr bool operator!= (const DenseMatrix<Other>& x) const
{
return !operator==(x);
}
//===== linear maps
template<class X, class Y>
constexpr void mv (const X& x, Y& y) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS((void*)(&x) != (void*)(&y));
DUNE_ASSERT_BOUNDS(xx.N() == M());
DUNE_ASSERT_BOUNDS(yy.N() == N());
using y_field_type = typename FieldTraits<Y>::field_type;
for (size_type i=0; i<rows(); ++i)
{
yy[i] = y_field_type(0);
for (size_type j=0; j<cols(); j++)
yy[i] += (*this)[i][j] * xx[j];
}
}
template< class X, class Y >
constexpr void mtv ( const X &x, Y &y ) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS((void*)(&x) != (void*)(&y));
DUNE_ASSERT_BOUNDS(xx.N() == N());
DUNE_ASSERT_BOUNDS(yy.N() == M());
using y_field_type = typename FieldTraits<Y>::field_type;
for(size_type i = 0; i < cols(); ++i)
{
yy[i] = y_field_type(0);
for(size_type j = 0; j < rows(); ++j)
yy[i] += (*this)[j][i] * xx[j];
}
}
template<class X, class Y>
constexpr void umv (const X& x, Y& y) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS(xx.N() == M());
DUNE_ASSERT_BOUNDS(yy.N() == N());
for (size_type i=0; i<rows(); ++i)
for (size_type j=0; j<cols(); j++)
yy[i] += (*this)[i][j] * xx[j];
}
template<class X, class Y>
constexpr void umtv (const X& x, Y& y) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS(xx.N() == N());
DUNE_ASSERT_BOUNDS(yy.N() == M());
for(size_type i = 0; i<rows(); ++i)
for (size_type j=0; j<cols(); j++)
yy[j] += (*this)[i][j]*xx[i];
}
template<class X, class Y>
constexpr void umhv (const X& x, Y& y) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS(xx.N() == N());
DUNE_ASSERT_BOUNDS(yy.N() == M());
for (size_type i=0; i<rows(); i++)
for (size_type j=0; j<cols(); j++)
yy[j] += conjugateComplex((*this)[i][j])*xx[i];
}
template<class X, class Y>
constexpr void mmv (const X& x, Y& y) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS(xx.N() == M());
DUNE_ASSERT_BOUNDS(yy.N() == N());
for (size_type i=0; i<rows(); i++)
for (size_type j=0; j<cols(); j++)
yy[i] -= (*this)[i][j] * xx[j];
}
template<class X, class Y>
constexpr void mmtv (const X& x, Y& y) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS(xx.N() == N());
DUNE_ASSERT_BOUNDS(yy.N() == M());
for (size_type i=0; i<rows(); i++)
for (size_type j=0; j<cols(); j++)
yy[j] -= (*this)[i][j]*xx[i];
}
template<class X, class Y>
constexpr void mmhv (const X& x, Y& y) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS(xx.N() == N());
DUNE_ASSERT_BOUNDS(yy.N() == M());
for (size_type i=0; i<rows(); i++)
for (size_type j=0; j<cols(); j++)
yy[j] -= conjugateComplex((*this)[i][j])*xx[i];
}
template<class X, class Y>
constexpr void usmv (const typename FieldTraits<Y>::field_type & alpha,
const X& x, Y& y) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS(xx.N() == M());
DUNE_ASSERT_BOUNDS(yy.N() == N());
for (size_type i=0; i<rows(); i++)
for (size_type j=0; j<cols(); j++)
yy[i] += alpha * (*this)[i][j] * xx[j];
}
template<class X, class Y>
constexpr void usmtv (const typename FieldTraits<Y>::field_type & alpha,
const X& x, Y& y) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS(xx.N() == N());
DUNE_ASSERT_BOUNDS(yy.N() == M());
for (size_type i=0; i<rows(); i++)
for (size_type j=0; j<cols(); j++)
yy[j] += alpha*(*this)[i][j]*xx[i];
}
template<class X, class Y>
constexpr void usmhv (const typename FieldTraits<Y>::field_type & alpha,
const X& x, Y& y) const
{
auto&& xx = Impl::asVector(x);
auto&& yy = Impl::asVector(y);
DUNE_ASSERT_BOUNDS(xx.N() == N());
DUNE_ASSERT_BOUNDS(yy.N() == M());
for (size_type i=0; i<rows(); i++)
for (size_type j=0; j<cols(); j++)
yy[j] +=
alpha*conjugateComplex((*this)[i][j])*xx[i];
}
//===== norms
constexpr typename FieldTraits<value_type>::real_type frobenius_norm () const
{
typename FieldTraits<value_type>::real_type sum=(0.0);
for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
return fvmeta::sqrt(sum);
}
constexpr typename FieldTraits<value_type>::real_type frobenius_norm2 () const
{
typename FieldTraits<value_type>::real_type sum=(0.0);
for (size_type i=0; i<rows(); ++i) sum += (*this)[i].two_norm2();
return sum;
}
template <typename vt = value_type,
typename std::enable_if<!HasNaN<vt>::value, int>::type = 0>
constexpr typename FieldTraits<vt>::real_type infinity_norm() const {
using real_type = typename FieldTraits<vt>::real_type;
using std::max;
real_type norm = 0;
for (auto const &x : *this) {
real_type const a = x.one_norm();
norm = max(a, norm);
}
return norm;
}
template <typename vt = value_type,
typename std::enable_if<!HasNaN<vt>::value, int>::type = 0>
constexpr typename FieldTraits<vt>::real_type infinity_norm_real() const {
using real_type = typename FieldTraits<vt>::real_type;
using std::max;
real_type norm = 0;
for (auto const &x : *this) {
real_type const a = x.one_norm_real();
norm = max(a, norm);
}
return norm;
}
template <typename vt = value_type,
typename std::enable_if<HasNaN<vt>::value, int>::type = 0>
constexpr typename FieldTraits<vt>::real_type infinity_norm() const {
using real_type = typename FieldTraits<vt>::real_type;
using std::max;
real_type norm = 0;
real_type is_nan = 1;
for (auto const &x : *this) {
real_type const a = x.one_norm();
norm = max(a, norm);
is_nan += a;
}
return norm * (is_nan / is_nan);
}
template <typename vt = value_type,
typename std::enable_if<HasNaN<vt>::value, int>::type = 0>
constexpr typename FieldTraits<vt>::real_type infinity_norm_real() const {
using real_type = typename FieldTraits<vt>::real_type;
using std::max;
real_type norm = 0;
real_type is_nan = 1;
for (auto const &x : *this) {
real_type const a = x.one_norm_real();
norm = max(a, norm);
is_nan += a;
}
return norm * (is_nan / is_nan);
}
//===== solve
template <class V1, class V2>
void solve (V1& x, const V2& b, bool doPivoting = true) const;
void invert(bool doPivoting = true);
field_type determinant (bool doPivoting = true) const;
template<typename M2>
MAT& leftmultiply (const DenseMatrix<M2>& M)
{
DUNE_ASSERT_BOUNDS(M.rows() == M.cols());
DUNE_ASSERT_BOUNDS(M.rows() == rows());
AutonomousValue<MAT> C(asImp());
for (size_type i=0; i<rows(); i++)
for (size_type j=0; j<cols(); j++) {
(*this)[i][j] = 0;
for (size_type k=0; k<rows(); k++)
(*this)[i][j] += M[i][k]*C[k][j];
}
return asImp();
}
template<typename M2>
MAT& rightmultiply (const DenseMatrix<M2>& M)
{
DUNE_ASSERT_BOUNDS(M.rows() == M.cols());
DUNE_ASSERT_BOUNDS(M.cols() == cols());
AutonomousValue<MAT> C(asImp());
for (size_type i=0; i<rows(); i++)
for (size_type j=0; j<cols(); j++) {
(*this)[i][j] = 0;
for (size_type k=0; k<cols(); k++)
(*this)[i][j] += C[i][k]*M[k][j];
}
return asImp();
}
#if 0
template<int l>
DenseMatrix<K,l,cols> leftmultiplyany (const FieldMatrix<K,l,rows>& M) const
{
FieldMatrix<K,l,cols> C;
for (size_type i=0; i<l; i++) {
for (size_type j=0; j<cols(); j++) {
C[i][j] = 0;
for (size_type k=0; k<rows(); k++)
C[i][j] += M[i][k]*(*this)[k][j];
}
}
return C;
}
template<int l>
FieldMatrix<K,rows,l> rightmultiplyany (const FieldMatrix<K,cols,l>& M) const
{
FieldMatrix<K,rows,l> C;
for (size_type i=0; i<rows(); i++) {
for (size_type j=0; j<l; j++) {
C[i][j] = 0;
for (size_type k=0; k<cols(); k++)
C[i][j] += (*this)[i][k]*M[k][j];
}
}
return C;
}
#endif
//===== sizes
constexpr size_type N () const
{
return rows();
}
constexpr size_type M () const
{
return cols();
}
constexpr size_type rows() const
{
return asImp().mat_rows();
}
constexpr size_type cols() const
{
return asImp().mat_cols();
}
//===== query
constexpr bool exists ([[maybe_unused]] size_type i, [[maybe_unused]] size_type j) const
{
DUNE_ASSERT_BOUNDS(i >= 0 && i < rows());
DUNE_ASSERT_BOUNDS(j >= 0 && j < cols());
return true;
}
protected:
#ifndef DOXYGEN
struct ElimPivot
{
ElimPivot(std::vector<simd_index_type> & pivot);
void swap(std::size_t i, simd_index_type j);
template<typename T>
void operator()(const T&, int, int)
{}
std::vector<simd_index_type> & pivot_;
};
template<typename V>
struct Elim
{
Elim(V& rhs);
void swap(std::size_t i, simd_index_type j);
void operator()(const typename V::field_type& factor, int k, int i);
V* rhs_;
};
struct ElimDet
{
ElimDet(field_type& sign) : sign_(sign)
{ sign_ = 1; }
void swap(std::size_t i, simd_index_type j)
{
sign_ *=
Simd::cond(simd_index_type(i) == j, field_type(1), field_type(-1));
}
void operator()(const field_type&, int, int)
{}
field_type& sign_;
};
#endif // DOXYGEN
template<class Func, class Mask>
static void luDecomposition(DenseMatrix<MAT>& A, Func func,
Mask &nonsingularLanes, bool throwEarly, bool doPivoting);
};
#ifndef DOXYGEN
template<typename MAT>
DenseMatrix<MAT>::ElimPivot::ElimPivot(std::vector<simd_index_type> & pivot)
: pivot_(pivot)
{
for(typename std::vector<size_type>::size_type i=0; i < pivot_.size(); ++i)
pivot_[i]=i;
}
template<typename MAT>
void DenseMatrix<MAT>::ElimPivot::swap(std::size_t i, simd_index_type j)
{
pivot_[i] =
Simd::cond(Simd::Scalar<simd_index_type>(i) == j, pivot_[i], j);
}
template<typename MAT>
template<typename V>
DenseMatrix<MAT>::Elim<V>::Elim(V& rhs)
: rhs_(&rhs)
{}
template<typename MAT>
template<typename V>
void DenseMatrix<MAT>::Elim<V>::swap(std::size_t i, simd_index_type j)
{
using std::swap;
// see the comment in luDecomposition()
for(std::size_t l = 0; l < Simd::lanes(j); ++l)
swap(Simd::lane(l, (*rhs_)[ i ]),
Simd::lane(l, (*rhs_)[Simd::lane(l, j)]));
}
template<typename MAT>
template<typename V>
void DenseMatrix<MAT>::
Elim<V>::operator()(const typename V::field_type& factor, int k, int i)
{
(*rhs_)[k] -= factor*(*rhs_)[i];
}
template<typename MAT>
template<typename Func, class Mask>
inline void DenseMatrix<MAT>::
luDecomposition(DenseMatrix<MAT>& A, Func func, Mask &nonsingularLanes,
bool throwEarly, bool doPivoting)
{
using std::max;
using std::swap;
typedef typename FieldTraits<value_type>::real_type real_type;
// LU decomposition of A in A
for (size_type i=0; i<A.rows(); i++) // loop over all rows
{
real_type pivmax = fvmeta::absreal(A[i][i]);
if (doPivoting)
{
// compute maximum of column
simd_index_type imax=i;
for (size_type k=i+1; k<A.rows(); k++)
{
auto abs = fvmeta::absreal(A[k][i]);
auto mask = abs > pivmax;
pivmax = Simd::cond(mask, abs, pivmax);
imax = Simd::cond(mask, simd_index_type(k), imax);
}
// swap rows
for (size_type j=0; j<A.rows(); j++)
{
// This is a swap operation where the second operand is scattered,
// and on top of that is also extracted from deep within a
// moderately complicated data structure (a DenseMatrix), where we
// can't assume much on the memory layout. On intel processors,
// the only instruction that might help us here is vgather, but it
// is unclear whether that is even faster than a software
// implementation, and we would also need vscatter which does not
// exist. So break vectorization here and do it manually.
for(std::size_t l = 0; l < Simd::lanes(A[i][j]); ++l)
swap(Simd::lane(l, A[i][j]),
Simd::lane(l, A[Simd::lane(l, imax)][j]));
}
func.swap(i, imax); // swap the pivot or rhs
}
// singular ?
nonsingularLanes = nonsingularLanes && (pivmax != real_type(0));
if (throwEarly) {
if(!Simd::allTrue(nonsingularLanes))
DUNE_THROW(FMatrixError, "matrix is singular");
}
else { // !throwEarly
if(!Simd::anyTrue(nonsingularLanes))
return;
}
// eliminate
for (size_type k=i+1; k<A.rows(); k++)
{
// in the simd case, A[i][i] may be close to zero in some lanes. Pray
// that the result is no worse than a quiet NaN.
field_type factor = A[k][i]/A[i][i];
A[k][i] = factor;
for (size_type j=i+1; j<A.rows(); j++)
A[k][j] -= factor*A[i][j];
func(factor, k, i);
}
}
}
template<typename MAT>
template <class V1, class V2>
inline void DenseMatrix<MAT>::solve(V1& x, const V2& b, bool doPivoting) const
{
using real_type = typename FieldTraits<value_type>::real_type;
// never mind those ifs, because they get optimized away
if (rows()!=cols())
DUNE_THROW(FMatrixError, "Can't solve for a " << rows() << "x" << cols() << " matrix!");
if (rows()==1) {
#ifdef DUNE_FMatrix_WITH_CHECKING
if (Simd::anyTrue(fvmeta::absreal((*this)[0][0])
< FMatrixPrecision<>::absolute_limit()))
DUNE_THROW(FMatrixError,"matrix is singular");
#endif
x[0] = b[0]/(*this)[0][0];
}
else if (rows()==2) {
field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
#ifdef DUNE_FMatrix_WITH_CHECKING
if (Simd::anyTrue(fvmeta::absreal(detinv)
< FMatrixPrecision<>::absolute_limit()))
DUNE_THROW(FMatrixError,"matrix is singular");
#endif
detinv = real_type(1.0)/detinv;
x[0] = detinv*((*this)[1][1]*b[0]-(*this)[0][1]*b[1]);
x[1] = detinv*((*this)[0][0]*b[1]-(*this)[1][0]*b[0]);
}
else if (rows()==3) {
field_type d = determinant(doPivoting);
#ifdef DUNE_FMatrix_WITH_CHECKING
if (Simd::anyTrue(fvmeta::absreal(d)
< FMatrixPrecision<>::absolute_limit()))
DUNE_THROW(FMatrixError,"matrix is singular");
#endif
x[0] = (b[0]*(*this)[1][1]*(*this)[2][2] - b[0]*(*this)[2][1]*(*this)[1][2]
- b[1] *(*this)[0][1]*(*this)[2][2] + b[1]*(*this)[2][1]*(*this)[0][2]
+ b[2] *(*this)[0][1]*(*this)[1][2] - b[2]*(*this)[1][1]*(*this)[0][2]) / d;
x[1] = ((*this)[0][0]*b[1]*(*this)[2][2] - (*this)[0][0]*b[2]*(*this)[1][2]
- (*this)[1][0] *b[0]*(*this)[2][2] + (*this)[1][0]*b[2]*(*this)[0][2]
+ (*this)[2][0] *b[0]*(*this)[1][2] - (*this)[2][0]*b[1]*(*this)[0][2]) / d;
x[2] = ((*this)[0][0]*(*this)[1][1]*b[2] - (*this)[0][0]*(*this)[2][1]*b[1]
- (*this)[1][0] *(*this)[0][1]*b[2] + (*this)[1][0]*(*this)[2][1]*b[0]
+ (*this)[2][0] *(*this)[0][1]*b[1] - (*this)[2][0]*(*this)[1][1]*b[0]) / d;
}
else {
V1& rhs = x; // use x to store rhs
rhs = b; // copy data
Elim<V1> elim(rhs);
AutonomousValue<MAT> A(asImp());
Simd::Mask<typename FieldTraits<value_type>::real_type>
nonsingularLanes(true);
AutonomousValue<MAT>::luDecomposition(A, elim, nonsingularLanes, true, doPivoting);
// backsolve
for(int i=rows()-1; i>=0; i--) {
for (size_type j=i+1; j<rows(); j++)
rhs[i] -= A[i][j]*x[j];
x[i] = rhs[i]/A[i][i];
}
}
}
template<typename MAT>
inline void DenseMatrix<MAT>::invert(bool doPivoting)
{
using real_type = typename FieldTraits<MAT>::real_type;
using std::swap;
// never mind those ifs, because they get optimized away
if (rows()!=cols())
DUNE_THROW(FMatrixError, "Can't invert a " << rows() << "x" << cols() << " matrix!");
if (rows()==1) {
#ifdef DUNE_FMatrix_WITH_CHECKING
if (Simd::anyTrue(fvmeta::absreal((*this)[0][0])
< FMatrixPrecision<>::absolute_limit()))
DUNE_THROW(FMatrixError,"matrix is singular");
#endif
(*this)[0][0] = real_type( 1 ) / (*this)[0][0];
}
else if (rows()==2) {
field_type detinv = (*this)[0][0]*(*this)[1][1]-(*this)[0][1]*(*this)[1][0];
#ifdef DUNE_FMatrix_WITH_CHECKING
if (Simd::anyTrue(fvmeta::absreal(detinv)
< FMatrixPrecision<>::absolute_limit()))
DUNE_THROW(FMatrixError,"matrix is singular");
#endif
detinv = real_type( 1 ) / detinv;
field_type temp=(*this)[0][0];
(*this)[0][0] = (*this)[1][1]*detinv;
(*this)[0][1] = -(*this)[0][1]*detinv;
(*this)[1][0] = -(*this)[1][0]*detinv;
(*this)[1][1] = temp*detinv;
}
else if (rows()==3)
{
using K = field_type;
// code generated by maple
K t4 = (*this)[0][0] * (*this)[1][1];
K t6 = (*this)[0][0] * (*this)[1][2];
K t8 = (*this)[0][1] * (*this)[1][0];
K t10 = (*this)[0][2] * (*this)[1][0];
K t12 = (*this)[0][1] * (*this)[2][0];
K t14 = (*this)[0][2] * (*this)[2][0];
K det = (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
K t17 = K(1.0)/det;
K matrix01 = (*this)[0][1];
K matrix00 = (*this)[0][0];
K matrix10 = (*this)[1][0];
K matrix11 = (*this)[1][1];
(*this)[0][0] = ((*this)[1][1] * (*this)[2][2] - (*this)[1][2] * (*this)[2][1])*t17;
(*this)[0][1] = -((*this)[0][1] * (*this)[2][2] - (*this)[0][2] * (*this)[2][1])*t17;
(*this)[0][2] = (matrix01 * (*this)[1][2] - (*this)[0][2] * (*this)[1][1])*t17;
(*this)[1][0] = -((*this)[1][0] * (*this)[2][2] - (*this)[1][2] * (*this)[2][0])*t17;
(*this)[1][1] = (matrix00 * (*this)[2][2] - t14) * t17;
(*this)[1][2] = -(t6-t10) * t17;
(*this)[2][0] = (matrix10 * (*this)[2][1] - matrix11 * (*this)[2][0]) * t17;
(*this)[2][1] = -(matrix00 * (*this)[2][1] - t12) * t17;
(*this)[2][2] = (t4-t8) * t17;
}
else {
using std::swap;
AutonomousValue<MAT> A(asImp());
std::vector<simd_index_type> pivot(rows());
Simd::Mask<typename FieldTraits<value_type>::real_type>
nonsingularLanes(true);
AutonomousValue<MAT>::luDecomposition(A, ElimPivot(pivot), nonsingularLanes, true, doPivoting);
auto& L=A;
auto& U=A;
// initialize inverse
*this=field_type();
for(size_type i=0; i<rows(); ++i)
(*this)[i][i]=1;
// L Y = I; multiple right hand sides
for (size_type i=0; i<rows(); i++)
for (size_type j=0; j<i; j++)
for (size_type k=0; k<rows(); k++)
(*this)[i][k] -= L[i][j]*(*this)[j][k];
// U A^{-1} = Y
for (size_type i=rows(); i>0;) {
--i;
for (size_type k=0; k<rows(); k++) {
for (size_type j=i+1; j<rows(); j++)
(*this)[i][k] -= U[i][j]*(*this)[j][k];
(*this)[i][k] /= U[i][i];
}
}
for(size_type i=rows(); i>0; ) {
--i;
for(std::size_t l = 0; l < Simd::lanes((*this)[0][0]); ++l)
{
std::size_t pi = Simd::lane(l, pivot[i]);
if(i!=pi)
for(size_type j=0; j<rows(); ++j)
swap(Simd::lane(l, (*this)[j][pi]),
Simd::lane(l, (*this)[j][ i]));
}
}
}
}
// implementation of the determinant
template<typename MAT>
inline typename DenseMatrix<MAT>::field_type
DenseMatrix<MAT>::determinant(bool doPivoting) const
{
// never mind those ifs, because they get optimized away
if (rows()!=cols())
DUNE_THROW(FMatrixError, "There is no determinant for a " << rows() << "x" << cols() << " matrix!");
if (rows()==1)
return (*this)[0][0];
if (rows()==2)
return (*this)[0][0]*(*this)[1][1] - (*this)[0][1]*(*this)[1][0];
if (rows()==3) {
// code generated by maple
field_type t4 = (*this)[0][0] * (*this)[1][1];
field_type t6 = (*this)[0][0] * (*this)[1][2];
field_type t8 = (*this)[0][1] * (*this)[1][0];
field_type t10 = (*this)[0][2] * (*this)[1][0];
field_type t12 = (*this)[0][1] * (*this)[2][0];
field_type t14 = (*this)[0][2] * (*this)[2][0];
return (t4*(*this)[2][2]-t6*(*this)[2][1]-t8*(*this)[2][2]+
t10*(*this)[2][1]+t12*(*this)[1][2]-t14*(*this)[1][1]);
}
AutonomousValue<MAT> A(asImp());
field_type det;
Simd::Mask<typename FieldTraits<value_type>::real_type>
nonsingularLanes(true);
AutonomousValue<MAT>::luDecomposition(A, ElimDet(det), nonsingularLanes, false, doPivoting);
det = Simd::cond(nonsingularLanes, det, field_type(0));
for (size_type i = 0; i < rows(); ++i)
det *= A[i][i];
return det;
}
#endif // DOXYGEN
namespace DenseMatrixHelp {
template <typename MAT, typename V1, typename V2>
static inline void multAssign(const DenseMatrix<MAT> &matrix, const DenseVector<V1> & x, DenseVector<V2> & ret)
{
DUNE_ASSERT_BOUNDS(x.size() == matrix.cols());
DUNE_ASSERT_BOUNDS(ret.size() == matrix.rows());
typedef typename DenseMatrix<MAT>::size_type size_type;
for(size_type i=0; i<matrix.rows(); ++i)
{
ret[i] = 0.0;
for(size_type j=0; j<matrix.cols(); ++j)
{
ret[i] += matrix[i][j]*x[j];
}
}
}
#if 0
template <typename K, int rows, int cols>
static inline void multAssignTransposed( const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x, FieldVector<K,cols> & ret)
{
typedef typename FieldMatrix<K,rows,cols>::size_type size_type;
for(size_type i=0; i<cols(); ++i)
{
ret[i] = 0.0;
for(size_type j=0; j<rows(); ++j)
ret[i] += matrix[j][i]*x[j];
}
}
template <typename K, int rows, int cols>
static inline FieldVector<K,rows> mult(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,cols> & x)
{
FieldVector<K,rows> ret;
multAssign(matrix,x,ret);
return ret;
}
template <typename K, int rows, int cols>
static inline FieldVector<K,cols> multTransposed(const FieldMatrix<K,rows,cols> &matrix, const FieldVector<K,rows> & x)
{
FieldVector<K,cols> ret;
multAssignTransposed( matrix, x, ret );
return ret;
}
#endif
} // end namespace DenseMatrixHelp
template<typename MAT>
std::ostream& operator<< (std::ostream& s, const DenseMatrix<MAT>& a)
{
for (typename DenseMatrix<MAT>::size_type i=0; i<a.rows(); i++)
s << a[i] << std::endl;
return s;
}
} // end namespace Dune
#endif
typetraits.hh
Traits for type conversions and type information.
precision.hh
Various precision settings for calculations with FieldMatrix and FieldVector.
exceptions.hh
A few common exception classes.
iterator.hh
classname.hh
A free function to provide the demangled class name of a given object or type as a string.
scalarvectorview.hh
Implements a scalar vector view wrapper around an existing scalar.
boundschecking.hh
Macro for wrapping boundary checks.
fvector.hh
Implements a vector constructed from a given type representing a field and a compile-time given size.
math.hh
Some useful basic math stuff.
DUNE_ASSERT_BOUNDS
#define DUNE_ASSERT_BOUNDS(cond)
If DUNE_CHECK_BOUNDS is defined: check if condition cond holds; otherwise, do nothing.
Definition boundschecking.hh:30
Dune::AutonomousValue
typename AutonomousValueType< T >::type AutonomousValue
Type free of internal references that T can be converted to.
Definition typetraits.hh:566
Dune::operator<<
std::ostream & operator<<(std::ostream &s, const bigunsignedint< k > &x)
Definition bigunsignedint.hh:301
DUNE_THROW
#define DUNE_THROW(E,...)
Definition exceptions.hh:314
Dune::Simd::anyTrue
bool anyTrue(const Mask &mask)
Whether any entry is true.
Definition simd/interface.hh:429
Dune::Simd::cond
V cond(M &&mask, const V &ifTrue, const V &ifFalse)
Like the ?: operator.
Definition simd/interface.hh:386
Dune::Simd::allTrue
bool allTrue(const Mask &mask)
Whether all entries are true.
Definition simd/interface.hh:439
Dune::Simd::Scalar
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition simd/interface.hh:235
Dune::Simd::Mask
Rebind< bool, V > Mask
Mask type type of some SIMD type.
Definition simd/interface.hh:289
Dune::Simd::Rebind
typename Overloads::RebindType< std::decay_t< S >, std::decay_t< V > >::type Rebind
Construct SIMD type with different scalar type.
Definition simd/interface.hh:252
Dune::Simd::lanes
constexpr std::size_t lanes()
Number of lanes in a SIMD type.
Definition simd/interface.hh:305
Dune::Simd::lane
decltype(auto) lane(std::size_t l, V &&v)
Extract an element of a SIMD type.
Definition simd/interface.hh:324
Dune::Simd::Overloads::mask
Mask< V > mask(ADLTag< 0, std::is_same< V, Mask< V > >::value >, const V &v)
implements Simd::mask()
Definition defaults.hh:153
Dune
Dune namespace
Definition alignedallocator.hh:13
Dune::sign
constexpr int sign(const T &val)
Return the sign of the value.
Definition math.hh:162
Dune::swap
void swap(T &v1, T &v2, bool mask)
Definition simd.hh:472
Dune::conjugateComplex
constexpr K conjugateComplex(const K &x)
compute conjugate complex of x
Definition math.hh:146
Dune::DenseMatrixHelp
Definition densematrix.hh:1170
Dune::DenseMatrixHelp::multAssign
static void multAssign(const DenseMatrix< MAT > &matrix, const DenseVector< V1 > &x, DenseVector< V2 > &ret)
calculates ret = matrix * x
Definition densematrix.hh:1174
Dune::Std::abs
constexpr T abs(T t)
Definition cmath.hh:27
Dune::IsNumber
Whether this type acts as a scalar in the context of (hierarchically blocked) containers.
Definition typetraits.hh:194
Dune::DenseMatrix
A dense n x m matrix.
Definition densematrix.hh:145
Dune::DenseMatrix::const_iterator
ConstIterator const_iterator
typedef for stl compliant access
Definition densematrix.hh:249
Dune::DenseMatrix::operator+=
constexpr derived_type & operator+=(const DenseMatrix< Other > &x)
vector space addition
Definition densematrix.hh:294
Dune::DenseMatrix::beforeBegin
constexpr Iterator beforeBegin()
Definition densematrix.hh:241
Dune::DenseMatrix::umhv
constexpr void umhv(const X &x, Y &y) const
y += A^H x
Definition densematrix.hh:437
Dune::DenseMatrix::operator*=
constexpr derived_type & operator*=(const field_type &k)
vector space multiplication with scalar
Definition densematrix.hh:326
Dune::DenseMatrix::usmhv
constexpr void usmhv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition densematrix.hh:517
Dune::DenseMatrix::operator-=
constexpr derived_type & operator-=(const DenseMatrix< Other > &x)
vector space subtraction
Definition densematrix.hh:317
Dune::DenseMatrix::solve
void solve(V1 &x, const V2 &b, bool doPivoting=true) const
Solve system A x = b.
Dune::DenseMatrix::field_type
Traits::value_type field_type
export the type representing the field
Definition densematrix.hh:165
Dune::DenseMatrix::beforeEnd
constexpr ConstIterator beforeEnd() const
Definition densematrix.hh:269
Dune::DenseMatrix::operator!=
constexpr bool operator!=(const DenseMatrix< Other > &x) const
Binary matrix incomparison.
Definition densematrix.hh:363
Dune::DenseMatrix::axpy
constexpr derived_type & axpy(const field_type &a, const DenseMatrix< Other > &x)
vector space axpy operation (*this += a x)
Definition densematrix.hh:343
Dune::DenseMatrix::ConstColIterator
std::remove_reference< const_row_reference >::type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition densematrix.hh:253
Dune::DenseMatrix::invert
void invert(bool doPivoting=true)
Compute inverse.
Dune::DenseMatrix::usmtv
constexpr void usmtv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition densematrix.hh:503
Dune::DenseMatrix::luDecomposition
static void luDecomposition(DenseMatrix< MAT > &A, Func func, Mask &nonsingularLanes, bool throwEarly, bool doPivoting)
do an LU-Decomposition on matrix A
Dune::DenseMatrix::block_type
Traits::value_type block_type
export the type representing the components
Definition densematrix.hh:168
Dune::DenseMatrix::cols
constexpr size_type cols() const
number of columns
Definition densematrix.hh:720
Dune::DenseMatrix::size
constexpr size_type size() const
size method (number of rows)
Definition densematrix.hh:205
Dune::DenseMatrix< DynamicMatrix< K > >::M
constexpr size_type M() const
Definition densematrix.hh:708
Dune::DenseMatrix::mmhv
constexpr void mmhv(const X &x, Y &y) const
y -= A^H x
Definition densematrix.hh:476
Dune::DenseMatrix::rightmultiply
MAT & rightmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the right to this matrix.
Definition densematrix.hh:650
Dune::DenseMatrix::beforeEnd
constexpr Iterator beforeEnd()
Definition densematrix.hh:234
Dune::DenseMatrix::beforeBegin
constexpr ConstIterator beforeBegin() const
Definition densematrix.hh:276
Dune::DenseMatrix::frobenius_norm
constexpr FieldTraits< value_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition densematrix.hh:533
Dune::DenseMatrix::value_type
Traits::value_type value_type
export the type representing the field
Definition densematrix.hh:162
Dune::DenseMatrix::end
constexpr Iterator end()
end iterator
Definition densematrix.hh:227
Dune::DenseMatrix::operator/=
constexpr derived_type & operator/=(const field_type &k)
vector space division by scalar
Definition densematrix.hh:334
Dune::DenseMatrix::operator-
constexpr derived_type operator-() const
Matrix negation.
Definition densematrix.hh:303
Dune::DenseMatrix::rows
constexpr size_type rows() const
number of rows
Definition densematrix.hh:714
Dune::DenseMatrix::mv
constexpr void mv(const X &x, Y &y) const
y = A x
Definition densematrix.hh:373
Dune::DenseMatrix::leftmultiply
MAT & leftmultiply(const DenseMatrix< M2 > &M)
Multiplies M from the left to this matrix.
Definition densematrix.hh:632
Dune::DenseMatrix::mtv
constexpr void mtv(const X &x, Y &y) const
y = A^T x
Definition densematrix.hh:392
Dune::DenseMatrix::infinity_norm
constexpr FieldTraits< vt >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition densematrix.hh:551
Dune::DenseMatrix::operator[]
constexpr row_reference operator[](size_type i)
random access
Definition densematrix.hh:194
Dune::DenseMatrix::umv
constexpr void umv(const X &x, Y &y) const
y += A x
Definition densematrix.hh:411
Dune::DenseMatrix::derived_type
Traits::derived_type derived_type
type of derived matrix class
Definition densematrix.hh:159
Dune::DenseMatrix::exists
constexpr bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition densematrix.hh:728
Dune::DenseMatrix::usmv
constexpr void usmv(const typename FieldTraits< Y >::field_type &alpha, const X &x, Y &y) const
y += alpha A x
Definition densematrix.hh:489
Dune::DenseMatrix::RowIterator
Iterator RowIterator
rename the iterators for easier access
Definition densematrix.hh:216
Dune::DenseMatrix< DynamicMatrix< K > >::blocklevel
static constexpr int blocklevel
Definition densematrix.hh:183
Dune::DenseMatrix::begin
constexpr ConstIterator begin() const
begin iterator
Definition densematrix.hh:256
Dune::DenseMatrix::umtv
constexpr void umtv(const X &x, Y &y) const
y += A^T x
Definition densematrix.hh:424
Dune::DenseMatrix::ConstIterator
DenseIterator< const DenseMatrix, const row_type, const_row_reference > ConstIterator
Iterator class for sequential access.
Definition densematrix.hh:247
Dune::DenseMatrix::mmtv
constexpr void mmtv(const X &x, Y &y) const
y -= A^T x
Definition densematrix.hh:463
Dune::DenseMatrix::end
constexpr ConstIterator end() const
end iterator
Definition densematrix.hh:262
Dune::DenseMatrix::frobenius_norm2
constexpr FieldTraits< value_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition densematrix.hh:541
Dune::DenseMatrix::row_type
Traits::row_type row_type
The type used to represent a row (must fulfill the Dune::DenseVector interface).
Definition densematrix.hh:174
Dune::DenseMatrix< DynamicMatrix< K > >::N
constexpr size_type N() const
Definition densematrix.hh:702
Dune::DenseMatrix::operator==
constexpr bool operator==(const DenseMatrix< Other > &x) const
Binary matrix comparison.
Definition densematrix.hh:353
Dune::DenseMatrix::size_type
Traits::size_type size_type
The type used for the index access and size operation.
Definition densematrix.hh:171
Dune::DenseMatrix::const_row_reference
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:180
Dune::DenseMatrix::ColIterator
std::remove_reference< row_reference >::type::Iterator ColIterator
rename the iterators for easier access
Definition densematrix.hh:218
Dune::DenseMatrix::row_reference
Traits::row_reference row_reference
The type used to represent a reference to a row (usually row_type &).
Definition densematrix.hh:177
Dune::DenseMatrix::DenseMatrix
friend class DenseMatrix
Definition densematrix.hh:153
Dune::DenseMatrix::mmv
constexpr void mmv(const X &x, Y &y) const
y -= A x
Definition densematrix.hh:450
Dune::DenseMatrix::iterator
Iterator iterator
typedef for stl compliant access
Definition densematrix.hh:214
Dune::DenseMatrix::infinity_norm_real
constexpr FieldTraits< vt >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition densematrix.hh:566
Dune::DenseMatrix::ConstRowIterator
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition densematrix.hh:251
Dune::DenseMatrix::begin
constexpr Iterator begin()
begin iterator
Definition densematrix.hh:221
Dune::DenseMatrix::operator=
constexpr derived_type & operator=(const RHS &rhs)
Definition densematrix.hh:284
Dune::DenseMatrix::Iterator
DenseIterator< DenseMatrix, row_type, row_reference > Iterator
Iterator class for sequential access.
Definition densematrix.hh:212
Dune::DenseMatrix::determinant
field_type determinant(bool doPivoting=true) const
calculates the determinant of this matrix
Dune::FieldTraits< DenseMatrix< M > >::real_type
const FieldTraits< typenameDenseMatVecTraits< M >::value_type >::real_type real_type
Definition densematrix.hh:36
Dune::FieldTraits< DenseMatrix< M > >::field_type
const FieldTraits< typenameDenseMatVecTraits< M >::value_type >::field_type field_type
Definition densematrix.hh:35
Dune::FieldMatrix
A dense n x m matrix.
Definition fmatrix.hh:117
Dune::FieldMatrix::size_type
Base::size_type size_type
Definition fmatrix.hh:128
Dune::FieldVector
vector space out of a tensor product of fields.
Definition fvector.hh:97
Dune::DenseMatrixAssigner
you have to specialize this structure for any type that should be assignable to a DenseMatrix
Definition densematrix.hh:61
Dune::FMatrixError
Error thrown if operations of a FieldMatrix fail.
Definition densematrix.hh:131
Dune::DenseVector
Interface for a class of dense vectors over a given field.
Definition densevector.hh:244
Dune::DenseVector::size
constexpr size_type size() const
size method
Definition densevector.hh:351
Dune::DenseIterator
Generic iterator class for dense vector and matrix implementations.
Definition densevector.hh:132
Dune::MathError
Default exception class for mathematical errors.
Definition exceptions.hh:335
Dune::FieldTraits
Definition ftraits.hh:26
Dune::FieldTraits::field_type
T field_type
export the type representing the field
Definition ftraits.hh:28
Dune::FieldTraits::real_type
T real_type
export the type representing the real type of the field
Definition ftraits.hh:30
Dune::DenseMatVecTraits
Definition matvectraits.hh:31
Dune::FMatrixPrecision::absolute_limit
static ctype absolute_limit()
return threshold to declare matrix singular
Definition precision.hh:28
Dune::type
simd.hh
Include file for users of the SIMD abstraction layer.