Loading...
Searching...
No Matches
scaledidmatrix.hh
Go to the documentation of this file.
// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
// vi: set et ts=4 sw=2 sts=2:
#ifndef DUNE_ISTL_SCALEDIDMATRIX_HH
#define DUNE_ISTL_SCALEDIDMATRIX_HH
#include <cmath>
#include <cstddef>
#include <complex>
#include <iostream>
#include <dune/common/exceptions.hh>
#include <dune/common/fmatrix.hh>
#include <dune/common/diagonalmatrix.hh>
#include <dune/common/ftraits.hh>
namespace Dune {
template<class K, int n>
class ScaledIdentityMatrix
{
typedef DiagonalMatrixWrapper< ScaledIdentityMatrix<K,n> > WrapperType;
public:
//===== type definitions and constants
typedef K field_type;
typedef K block_type;
typedef std::size_t size_type;
typedef DiagonalRowVector<K,n> row_type;
typedef row_type reference;
typedef DiagonalRowVectorConst<K,n> const_row_type;
typedef const_row_type const_reference;
static constexpr std::integral_constant<size_type,n> rows = {};
static constexpr std::integral_constant<size_type,n> cols = {};
//===== constructors
ScaledIdentityMatrix () {}
ScaledIdentityMatrix (const K& k)
: p_(k)
{}
ScaledIdentityMatrix& operator= (const K& k)
{
p_ = k;
return *this;
}
bool identical(const ScaledIdentityMatrix<K,n>& other) const
{
return (this==&other);
}
//===== iterator interface to rows of the matrix
typedef ContainerWrapperIterator<const WrapperType, reference, reference> Iterator;
typedef Iterator iterator;
typedef Iterator RowIterator;
typedef typename row_type::Iterator ColIterator;
Iterator begin ()
{
return Iterator(WrapperType(this),0);
}
Iterator end ()
{
return Iterator(WrapperType(this),n);
}
Iterator beforeEnd ()
{
return Iterator(WrapperType(this),n-1);
}
Iterator beforeBegin ()
{
return Iterator(WrapperType(this),-1);
}
typedef ContainerWrapperIterator<const WrapperType, const_reference, const_reference> ConstIterator;
typedef ConstIterator const_iterator;
typedef ConstIterator ConstRowIterator;
typedef typename const_row_type::ConstIterator ConstColIterator;
ConstIterator begin () const
{
return ConstIterator(WrapperType(this),0);
}
ConstIterator end () const
{
return ConstIterator(WrapperType(this),n);
}
ConstIterator beforeEnd() const
{
return ConstIterator(WrapperType(this),n-1);
}
ConstIterator beforeBegin () const
{
return ConstIterator(WrapperType(this),-1);
}
ScaledIdentityMatrix& operator+= (const ScaledIdentityMatrix& y)
{
p_ += y.p_;
return *this;
}
ScaledIdentityMatrix& operator-= (const ScaledIdentityMatrix& y)
{
p_ -= y.p_;
return *this;
}
ScaledIdentityMatrix& operator+= (const K& k)
{
p_ += k;
return *this;
}
ScaledIdentityMatrix& operator-= (const K& k)
{
p_ -= k;
return *this;
}
ScaledIdentityMatrix& operator*= (const K& k)
{
p_ *= k;
return *this;
}
ScaledIdentityMatrix& operator/= (const K& k)
{
p_ /= k;
return *this;
}
template <class Scalar,
std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
friend auto operator* ( const ScaledIdentityMatrix& matrix, Scalar scalar)
{
return ScaledIdentityMatrix<typename PromotionTraits<K,Scalar>::PromotedType, n>{matrix.scalar()*scalar};
}
template <class Scalar,
std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
friend auto operator* (Scalar scalar, const ScaledIdentityMatrix& matrix)
{
return ScaledIdentityMatrix<typename PromotionTraits<Scalar,K>::PromotedType, n>{scalar*matrix.scalar()};
}
template <class OtherScalar>
requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
friend auto& operator+= (FieldMatrix<OtherScalar,n,n>& fieldMatrix,
const ScaledIdentityMatrix& matrix)
{
for (int i=0; i<n; i++)
fieldMatrix[i][i] += matrix.p_;
return fieldMatrix;
}
template <class OtherScalar>
requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
friend auto operator+ (const FieldMatrix<OtherScalar,n,n>& fieldMatrix,
const ScaledIdentityMatrix& matrix)
{
using Result = FieldMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,n,n>;
Result result = fieldMatrix;
result += matrix;
return result;
}
template <class OtherScalar>
requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
friend auto operator+ (const ScaledIdentityMatrix& matrix,
const FieldMatrix<OtherScalar,n,n>& fieldMatrix)
{
return fieldMatrix + matrix;
}
template <class OtherScalar>
requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
friend auto operator+= (DiagonalMatrix<OtherScalar,n>& diagonalMatrix,
const ScaledIdentityMatrix& matrix)
{
for (std::size_t i=0; i<n; i++)
diagonalMatrix.diagonal(i) += matrix.p_;
return diagonalMatrix;
}
template <class OtherScalar>
requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
friend auto operator+ (const DiagonalMatrix<OtherScalar,n>& diagonalMatrix,
const ScaledIdentityMatrix& matrix)
{
using Result = DiagonalMatrix<typename PromotionTraits<K,OtherScalar>::PromotedType,n>;
Result result = diagonalMatrix;
result += matrix;
return result;
}
template <class OtherScalar>
requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
friend auto operator+ (const ScaledIdentityMatrix& matrix,
const DiagonalMatrix<OtherScalar,n>& diagonalMatrix)
{
return diagonalMatrix + matrix;
}
//===== comparison ops
bool operator==(const ScaledIdentityMatrix& other) const
{
return p_==other.scalar();
}
bool operator!=(const ScaledIdentityMatrix& other) const
{
return p_!=other.scalar();
}
template<class X, class Y>
void mv (const X& x, Y& y) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
#endif
for (size_type i=0; i<n; ++i)
y[i] = p_ * x[i];
}
template<class X, class Y>
void mtv (const X& x, Y& y) const
{
mv(x, y);
}
template<class X, class Y>
void umv (const X& x, Y& y) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
#endif
for (size_type i=0; i<n; ++i)
y[i] += p_ * x[i];
}
template<class X, class Y>
void umtv (const X& x, Y& y) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
#endif
for (size_type i=0; i<n; ++i)
y[i] += p_ * x[i];
}
template<class X, class Y>
void umhv (const X& x, Y& y) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
#endif
for (size_type i=0; i<n; i++)
y[i] += conjugateComplex(p_)*x[i];
}
template<class X, class Y>
void mmv (const X& x, Y& y) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
#endif
for (size_type i=0; i<n; ++i)
y[i] -= p_ * x[i];
}
template<class X, class Y>
void mmtv (const X& x, Y& y) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
#endif
for (size_type i=0; i<n; ++i)
y[i] -= p_ * x[i];
}
template<class X, class Y>
void mmhv (const X& x, Y& y) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
#endif
for (size_type i=0; i<n; i++)
y[i] -= conjugateComplex(p_)*x[i];
}
template<class X, class Y>
void usmv (const K& alpha, const X& x, Y& y) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
#endif
for (size_type i=0; i<n; i++)
y[i] += alpha * p_ * x[i];
}
template<class X, class Y>
void usmtv (const K& alpha, const X& x, Y& y) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
#endif
for (size_type i=0; i<n; i++)
y[i] += alpha * p_ * x[i];
}
template<class X, class Y>
void usmhv (const K& alpha, const X& x, Y& y) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
#endif
for (size_type i=0; i<n; i++)
y[i] += alpha * conjugateComplex(p_) * x[i];
}
typename FieldTraits<field_type>::real_type frobenius_norm () const
{
return fvmeta::sqrt(n*p_*p_);
}
typename FieldTraits<field_type>::real_type frobenius_norm2 () const
{
return n*p_*p_;
}
typename FieldTraits<field_type>::real_type infinity_norm () const
{
return std::abs(p_);
}
typename FieldTraits<field_type>::real_type infinity_norm_real () const
{
return fvmeta::absreal(p_);
}
template<class V>
void solve (V& x, const V& b) const
{
for (int i=0; i<n; i++)
x[i] = b[i]/p_;
}
void invert()
{
p_ = 1/p_;
}
K determinant () const {
return std::pow(p_,n);
}
//===== sizes
size_type N () const
{
return n;
}
size_type M () const
{
return n;
}
//===== query
bool exists (size_type i, size_type j) const
{
#ifdef DUNE_FMatrix_WITH_CHECKING
if (i<0 || i>=n) DUNE_THROW(FMatrixError,"row index out of range");
if (j<0 || j>=n) DUNE_THROW(FMatrixError,"column index out of range");
#endif
return i==j;
}
//===== conversion operator
friend std::ostream& operator<< (std::ostream& s, const ScaledIdentityMatrix<K,n>& a)
{
for (size_type i=0; i<n; i++) {
for (size_type j=0; j<n; j++)
s << ((i==j) ? a.p_ : 0) << " ";
s << std::endl;
}
return s;
}
reference operator[](size_type i)
{
return reference(const_cast<K*>(&p_), i);
}
const_reference operator[](size_type i) const
{
return const_reference(const_cast<K*>(&p_), i);
}
const K& diagonal(size_type /*i*/) const
{
return p_;
}
K& diagonal(size_type /*i*/)
{
return p_;
}
const K& scalar() const
{
return p_;
}
K& scalar()
{
return p_;
}
private:
// the data, very simply a single number
K p_;
};
template <class DenseMatrix, class field, int N>
struct DenseMatrixAssigner<DenseMatrix, ScaledIdentityMatrix<field, N>> {
static void apply(DenseMatrix& denseMatrix,
ScaledIdentityMatrix<field, N> const& rhs) {
assert(denseMatrix.M() == N);
assert(denseMatrix.N() == N);
denseMatrix = field(0);
for (int i = 0; i < N; ++i)
denseMatrix[i][i] = rhs.scalar();
}
};
template<class K, int n>
struct FieldTraits< ScaledIdentityMatrix<K, n> >
{
using field_type = typename ScaledIdentityMatrix<K, n>::field_type;
using real_type = typename FieldTraits<field_type>::real_type;
};
} // end namespace
#endif
Dune
Definition allocator.hh:11
Dune::FieldMatrix
Definition matrixutils.hh:27
Dune::ScaledIdentityMatrix
A multiple of the identity matrix of static size.
Definition scaledidmatrix.hh:42
Dune::ScaledIdentityMatrix::usmhv
void usmhv(const K &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition scaledidmatrix.hh:448
Dune::ScaledIdentityMatrix::mmtv
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition scaledidmatrix.hh:400
Dune::ScaledIdentityMatrix::operator-=
ScaledIdentityMatrix & operator-=(const ScaledIdentityMatrix &y)
Vector space subtraction.
Definition scaledidmatrix.hh:190
Dune::ScaledIdentityMatrix::ConstColIterator
const_row_type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition scaledidmatrix.hh:150
Dune::ScaledIdentityMatrix::end
ConstIterator end() const
end iterator
Definition scaledidmatrix.hh:159
Dune::ScaledIdentityMatrix::beforeBegin
Iterator beforeBegin()
Definition scaledidmatrix.hh:137
Dune::ScaledIdentityMatrix::rows
static constexpr std::integral_constant< size_type, n > rows
The number of matrix rows.
Definition scaledidmatrix.hh:73
Dune::ScaledIdentityMatrix::operator!=
bool operator!=(const ScaledIdentityMatrix &other) const
Test for inequality.
Definition scaledidmatrix.hh:322
Dune::ScaledIdentityMatrix::mmv
void mmv(const X &x, Y &y) const
y -= A x
Definition scaledidmatrix.hh:388
Dune::ScaledIdentityMatrix::size_type
std::size_t size_type
The type used for the indices and sizes.
Definition scaledidmatrix.hh:55
Dune::ScaledIdentityMatrix::usmv
void usmv(const K &alpha, const X &x, Y &y) const
y += alpha A x
Definition scaledidmatrix.hh:424
Dune::ScaledIdentityMatrix::ColIterator
row_type::Iterator ColIterator
rename the iterators for easier access
Definition scaledidmatrix.hh:114
Dune::ScaledIdentityMatrix::const_reference
const_row_type const_reference
Definition scaledidmatrix.hh:70
Dune::ScaledIdentityMatrix::mv
void mv(const X &x, Y &y) const
y = A x
Definition scaledidmatrix.hh:333
Dune::ScaledIdentityMatrix::umtv
void umtv(const X &x, Y &y) const
y += A^T x
Definition scaledidmatrix.hh:364
Dune::ScaledIdentityMatrix::umhv
void umhv(const X &x, Y &y) const
y += A^H x
Definition scaledidmatrix.hh:376
Dune::ScaledIdentityMatrix::row_type
DiagonalRowVector< K, n > row_type
Type of a single matrix row.
Definition scaledidmatrix.hh:65
Dune::ScaledIdentityMatrix::Iterator
ContainerWrapperIterator< const WrapperType, reference, reference > Iterator
Iterator class for sequential access.
Definition scaledidmatrix.hh:108
Dune::ScaledIdentityMatrix::beforeEnd
Iterator beforeEnd()
Definition scaledidmatrix.hh:130
Dune::ScaledIdentityMatrix::determinant
K determinant() const
Calculates the determinant of this matrix.
Definition scaledidmatrix.hh:514
Dune::ScaledIdentityMatrix::field_type
K field_type
The type representing numbers.
Definition scaledidmatrix.hh:49
Dune::ScaledIdentityMatrix::usmtv
void usmtv(const K &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition scaledidmatrix.hh:436
Dune::ScaledIdentityMatrix::end
Iterator end()
end iterator
Definition scaledidmatrix.hh:123
Dune::ScaledIdentityMatrix::iterator
Iterator iterator
typedef for stl compliant access
Definition scaledidmatrix.hh:110
Dune::ScaledIdentityMatrix::cols
static constexpr std::integral_constant< size_type, n > cols
The number of matrix columns.
Definition scaledidmatrix.hh:76
Dune::ScaledIdentityMatrix::scalar
const K & scalar() const
Get const reference to the scalar diagonal value.
Definition scaledidmatrix.hh:583
Dune::ScaledIdentityMatrix::operator+
friend auto operator+(const FieldMatrix< OtherScalar, n, n > &fieldMatrix, const ScaledIdentityMatrix &matrix)
Addition of ScaledIdentityMatrix to FieldMatrix.
Definition scaledidmatrix.hh:261
Dune::ScaledIdentityMatrix::umv
void umv(const X &x, Y &y) const
y += A x
Definition scaledidmatrix.hh:352
Dune::ScaledIdentityMatrix::diagonal
const K & diagonal(size_type) const
Get const reference to diagonal entry.
Definition scaledidmatrix.hh:570
Dune::ScaledIdentityMatrix::operator=
ScaledIdentityMatrix & operator=(const K &k)
Assignment from scalar.
Definition scaledidmatrix.hh:91
Dune::ScaledIdentityMatrix::ConstIterator
ContainerWrapperIterator< const WrapperType, const_reference, const_reference > ConstIterator
Iterator class for sequential access.
Definition scaledidmatrix.hh:144
Dune::ScaledIdentityMatrix::diagonal
K & diagonal(size_type)
Get reference to diagonal entry.
Definition scaledidmatrix.hh:576
Dune::ScaledIdentityMatrix::solve
void solve(V &x, const V &b) const
Solve linear system A x = b.
Definition scaledidmatrix.hh:500
Dune::ScaledIdentityMatrix::exists
bool exists(size_type i, size_type j) const
Return true if (i,j) is in the matrix pattern, i.e., if i==j.
Definition scaledidmatrix.hh:535
Dune::ScaledIdentityMatrix::RowIterator
Iterator RowIterator
rename the iterators for easier access
Definition scaledidmatrix.hh:112
Dune::ScaledIdentityMatrix::const_iterator
ConstIterator const_iterator
typedef for stl compliant access
Definition scaledidmatrix.hh:146
Dune::ScaledIdentityMatrix::ScaledIdentityMatrix
ScaledIdentityMatrix()
Default constructor.
Definition scaledidmatrix.hh:81
Dune::ScaledIdentityMatrix::operator==
bool operator==(const ScaledIdentityMatrix &other) const
Test for equality.
Definition scaledidmatrix.hh:316
Dune::ScaledIdentityMatrix::beforeBegin
ConstIterator beforeBegin() const
Definition scaledidmatrix.hh:173
Dune::ScaledIdentityMatrix::operator/=
ScaledIdentityMatrix & operator/=(const K &k)
Vector space division by scalar.
Definition scaledidmatrix.hh:224
Dune::ScaledIdentityMatrix::operator<<
friend std::ostream & operator<<(std::ostream &s, const ScaledIdentityMatrix< K, n > &a)
Sends the matrix to an output stream.
Definition scaledidmatrix.hh:547
Dune::ScaledIdentityMatrix::frobenius_norm2
FieldTraits< field_type >::real_type frobenius_norm2() const
The square of the Frobenius norm.
Definition scaledidmatrix.hh:472
Dune::ScaledIdentityMatrix::frobenius_norm
FieldTraits< field_type >::real_type frobenius_norm() const
The Frobenius norm.
Definition scaledidmatrix.hh:466
Dune::ScaledIdentityMatrix::infinity_norm
FieldTraits< field_type >::real_type infinity_norm() const
The row sum norm.
Definition scaledidmatrix.hh:482
Dune::ScaledIdentityMatrix::ConstRowIterator
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition scaledidmatrix.hh:148
Dune::ScaledIdentityMatrix::beforeEnd
ConstIterator beforeEnd() const
Definition scaledidmatrix.hh:166
Dune::ScaledIdentityMatrix::M
size_type M() const
number of blocks in column direction
Definition scaledidmatrix.hh:527
Dune::ScaledIdentityMatrix::operator[]
const_reference operator[](size_type i) const
Return const_reference object as row replacement.
Definition scaledidmatrix.hh:564
Dune::ScaledIdentityMatrix::ScaledIdentityMatrix
ScaledIdentityMatrix(const K &k)
Constructor initializing the whole matrix with a scalar.
Definition scaledidmatrix.hh:85
Dune::ScaledIdentityMatrix::operator*
friend auto operator*(const ScaledIdentityMatrix &matrix, Scalar scalar)
Vector space multiplication with scalar.
Definition scaledidmatrix.hh:233
Dune::ScaledIdentityMatrix::infinity_norm_real
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition scaledidmatrix.hh:488
Dune::ScaledIdentityMatrix::operator*=
ScaledIdentityMatrix & operator*=(const K &k)
Vector space multiplication with scalar.
Definition scaledidmatrix.hh:217
Dune::ScaledIdentityMatrix::identical
bool identical(const ScaledIdentityMatrix< K, n > &other) const
Check if matrix is identical to other matrix.
Definition scaledidmatrix.hh:101
Dune::ScaledIdentityMatrix::invert
void invert()
Compute inverse.
Definition scaledidmatrix.hh:508
Dune::ScaledIdentityMatrix::begin
Iterator begin()
begin iterator
Definition scaledidmatrix.hh:117
Dune::ScaledIdentityMatrix::scalar
K & scalar()
Get reference to the scalar diagonal value.
Definition scaledidmatrix.hh:590
Dune::ScaledIdentityMatrix::reference
row_type reference
Definition scaledidmatrix.hh:66
Dune::ScaledIdentityMatrix::block_type
K block_type
The type representing matrix entries (which may be matrices themselves).
Definition scaledidmatrix.hh:52
Dune::ScaledIdentityMatrix::mmhv
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition scaledidmatrix.hh:412
Dune::ScaledIdentityMatrix::operator+=
ScaledIdentityMatrix & operator+=(const ScaledIdentityMatrix &y)
Vector space addition.
Definition scaledidmatrix.hh:183
Dune::ScaledIdentityMatrix::const_row_type
DiagonalRowVectorConst< K, n > const_row_type
Const type of a single row.
Definition scaledidmatrix.hh:69
Dune::ScaledIdentityMatrix::mtv
void mtv(const X &x, Y &y) const
y = A^T x
Definition scaledidmatrix.hh:345
Dune::ScaledIdentityMatrix::operator[]
reference operator[](size_type i)
Return reference object as row replacement.
Definition scaledidmatrix.hh:558
Dune::ScaledIdentityMatrix::begin
ConstIterator begin() const
begin iterator
Definition scaledidmatrix.hh:153
Dune::ScaledIdentityMatrix::N
size_type N() const
number of blocks in row direction
Definition scaledidmatrix.hh:521
Dune::DenseMatrixAssigner< DenseMatrix, ScaledIdentityMatrix< field, N > >::apply
static void apply(DenseMatrix &denseMatrix, ScaledIdentityMatrix< field, N > const &rhs)
Definition scaledidmatrix.hh:603
Dune::FieldTraits< ScaledIdentityMatrix< K, n > >::field_type
typename ScaledIdentityMatrix< K, n >::field_type field_type
Definition scaledidmatrix.hh:616
Dune::FieldTraits< ScaledIdentityMatrix< K, n > >::real_type
typename FieldTraits< field_type >::real_type real_type
Definition scaledidmatrix.hh:617