3#ifndef DUNE_ISTL_SCALEDIDMATRIX_HH
4#define DUNE_ISTL_SCALEDIDMATRIX_HH
16#include <dune/common/exceptions.hh>
17#include <dune/common/fmatrix.hh>
18#include <dune/common/diagonalmatrix.hh>
19#include <dune/common/ftraits.hh>
26 template<
class K,
int n>
29 typedef DiagonalMatrixWrapper< ScaledIdentityMatrix<K,n> > WrapperType;
44 [[deprecated(
"Use free function blockLevel(). Will be removed after 2.8.")]]
82 return (
this==&other);
87 typedef ContainerWrapperIterator<const WrapperType, reference, reference>
Iterator;
98 return Iterator(WrapperType(
this),0);
104 return Iterator(WrapperType(
this),n);
111 return Iterator(WrapperType(
this),n-1);
118 return Iterator(WrapperType(
this),-1);
123 typedef ContainerWrapperIterator<const WrapperType, const_reference, const_reference>
ConstIterator;
205 return p_==other.
scalar();
211 return p_!=other.
scalar();
217 template<
class X,
class Y>
218 void mv (
const X& x, Y& y)
const
220#ifdef DUNE_FMatrix_WITH_CHECKING
221 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
222 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
229 template<
class X,
class Y>
230 void mtv (
const X& x, Y& y)
const
236 template<
class X,
class Y>
237 void umv (
const X& x, Y& y)
const
239#ifdef DUNE_FMatrix_WITH_CHECKING
240 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
241 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
248 template<
class X,
class Y>
249 void umtv (
const X& x, Y& y)
const
251#ifdef DUNE_FMatrix_WITH_CHECKING
252 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
253 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
260 template<
class X,
class Y>
261 void umhv (
const X& x, Y& y)
const
263#ifdef DUNE_FMatrix_WITH_CHECKING
264 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
265 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
268 y[i] += conjugateComplex(p_)*x[i];
272 template<
class X,
class Y>
273 void mmv (
const X& x, Y& y)
const
275#ifdef DUNE_FMatrix_WITH_CHECKING
276 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
277 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
284 template<
class X,
class Y>
285 void mmtv (
const X& x, Y& y)
const
287#ifdef DUNE_FMatrix_WITH_CHECKING
288 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
289 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
296 template<
class X,
class Y>
297 void mmhv (
const X& x, Y& y)
const
299#ifdef DUNE_FMatrix_WITH_CHECKING
300 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
301 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
304 y[i] -= conjugateComplex(p_)*x[i];
308 template<
class X,
class Y>
309 void usmv (
const K& alpha,
const X& x, Y& y)
const
311#ifdef DUNE_FMatrix_WITH_CHECKING
312 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
313 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
316 y[i] += alpha * p_ * x[i];
320 template<
class X,
class Y>
321 void usmtv (
const K& alpha,
const X& x, Y& y)
const
323#ifdef DUNE_FMatrix_WITH_CHECKING
324 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
325 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
328 y[i] += alpha * p_ * x[i];
332 template<
class X,
class Y>
333 void usmhv (
const K& alpha,
const X& x, Y& y)
const
335#ifdef DUNE_FMatrix_WITH_CHECKING
336 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
337 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
340 y[i] += alpha * conjugateComplex(p_) * x[i];
348 return fvmeta::sqrt(n*p_*p_);
366 return fvmeta::absreal(p_);
376 for (
int i=0; i<n; i++)
389 return std::pow(p_,n);
411#ifdef DUNE_FMatrix_WITH_CHECKING
412 if (i<0 || i>=n) DUNE_THROW(FMatrixError,
"row index out of range");
413 if (j<0 || j>=n) DUNE_THROW(FMatrixError,
"column index out of range");
425 s << ((i==j) ? a.p_ : 0) <<
" ";
434 return reference(
const_cast<K*
>(&p_), i);
475 template <
class DenseMatrix,
class field,
int N>
477 static void apply(DenseMatrix& denseMatrix,
479 assert(denseMatrix.M() == N);
480 assert(denseMatrix.N() == N);
481 denseMatrix = field(0);
482 for (
int i = 0; i < N; ++i)
483 denseMatrix[i][i] = rhs.
scalar();
Definition: allocator.hh:9
A multiple of the identity matrix of static size.
Definition: scaledidmatrix.hh:28
void usmhv(const K &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: scaledidmatrix.hh:333
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: scaledidmatrix.hh:285
ScaledIdentityMatrix & operator-=(const ScaledIdentityMatrix &y)
vector space subtraction
Definition: scaledidmatrix.hh:167
const_row_type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:129
ConstIterator end() const
end iterator
Definition: scaledidmatrix.hh:138
Iterator beforeBegin()
Definition: scaledidmatrix.hh:116
bool operator!=(const ScaledIdentityMatrix &other) const
incomparison operator
Definition: scaledidmatrix.hh:209
void mmv(const X &x, Y &y) const
y -= A x
Definition: scaledidmatrix.hh:273
std::size_t size_type
The type used for the index access and size operations.
Definition: scaledidmatrix.hh:41
void usmv(const K &alpha, const X &x, Y &y) const
y += alpha A x
Definition: scaledidmatrix.hh:309
row_type::Iterator ColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:93
const_row_type const_reference
Definition: scaledidmatrix.hh:51
void mv(const X &x, Y &y) const
y = A x
Definition: scaledidmatrix.hh:218
void umtv(const X &x, Y &y) const
y += A^T x
Definition: scaledidmatrix.hh:249
void umhv(const X &x, Y &y) const
y += A^H x
Definition: scaledidmatrix.hh:261
DiagonalRowVector< K, n > row_type
Each row is implemented by a field vector.
Definition: scaledidmatrix.hh:48
ContainerWrapperIterator< const WrapperType, reference, reference > Iterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:87
Iterator beforeEnd()
Definition: scaledidmatrix.hh:109
K determinant() const
calculates the determinant of this matrix
Definition: scaledidmatrix.hh:388
K field_type
export the type representing the field
Definition: scaledidmatrix.hh:35
void usmtv(const K &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: scaledidmatrix.hh:321
Iterator end()
end iterator
Definition: scaledidmatrix.hh:102
Iterator iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:89
const K & scalar() const
Get const reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:457
void umv(const X &x, Y &y) const
y += A x
Definition: scaledidmatrix.hh:237
static constexpr std::size_t blocklevel
We are at the leaf of the block recursion.
Definition: scaledidmatrix.hh:45
const K & diagonal(size_type) const
Get const reference to diagonal entry.
Definition: scaledidmatrix.hh:444
@ rows
The number of rows.
Definition: scaledidmatrix.hh:56
@ cols
The number of columns.
Definition: scaledidmatrix.hh:58
ScaledIdentityMatrix & operator=(const K &k)
Definition: scaledidmatrix.hh:73
ContainerWrapperIterator< const WrapperType, const_reference, const_reference > ConstIterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:123
K & diagonal(size_type)
Get reference to diagonal entry.
Definition: scaledidmatrix.hh:450
void solve(V &x, const V &b) const
Solve system A x = b.
Definition: scaledidmatrix.hh:374
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: scaledidmatrix.hh:409
Iterator RowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:91
ConstIterator const_iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:125
ScaledIdentityMatrix()
Default constructor.
Definition: scaledidmatrix.hh:64
bool operator==(const ScaledIdentityMatrix &other) const
comparison operator
Definition: scaledidmatrix.hh:203
ConstIterator beforeBegin() const
Definition: scaledidmatrix.hh:152
ScaledIdentityMatrix & operator/=(const K &k)
vector space division by scalar
Definition: scaledidmatrix.hh:194
friend std::ostream & operator<<(std::ostream &s, const ScaledIdentityMatrix< K, n > &a)
Sends the matrix to an output stream.
Definition: scaledidmatrix.hh:421
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: scaledidmatrix.hh:352
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: scaledidmatrix.hh:346
FieldTraits< field_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: scaledidmatrix.hh:358
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:127
ConstIterator beforeEnd() const
Definition: scaledidmatrix.hh:145
size_type M() const
number of blocks in column direction
Definition: scaledidmatrix.hh:401
const_reference operator[](size_type i) const
Return const_reference object as row replacement.
Definition: scaledidmatrix.hh:438
ScaledIdentityMatrix(const K &k)
Constructor initializing the whole matrix with a scalar.
Definition: scaledidmatrix.hh:68
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: scaledidmatrix.hh:364
ScaledIdentityMatrix & operator*=(const K &k)
vector space multiplication with scalar
Definition: scaledidmatrix.hh:187
bool identical(const ScaledIdentityMatrix< K, n > &other) const
Definition: scaledidmatrix.hh:80
void invert()
Compute inverse.
Definition: scaledidmatrix.hh:382
Iterator begin()
begin iterator
Definition: scaledidmatrix.hh:96
K & scalar()
Get reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:464
row_type reference
Definition: scaledidmatrix.hh:49
K block_type
export the type representing the components
Definition: scaledidmatrix.hh:38
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: scaledidmatrix.hh:297
ScaledIdentityMatrix & operator+=(const ScaledIdentityMatrix &y)
vector space addition
Definition: scaledidmatrix.hh:160
DiagonalRowVectorConst< K, n > const_row_type
Definition: scaledidmatrix.hh:50
void mtv(const X &x, Y &y) const
y = A^T x
Definition: scaledidmatrix.hh:230
reference operator[](size_type i)
Return reference object as row replacement.
Definition: scaledidmatrix.hh:432
ConstIterator begin() const
begin iterator
Definition: scaledidmatrix.hh:132
size_type N() const
number of blocks in row direction
Definition: scaledidmatrix.hh:395
static void apply(DenseMatrix &denseMatrix, ScaledIdentityMatrix< field, N > const &rhs)
Definition: scaledidmatrix.hh:477