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;
84 return (
this==&other);
89 typedef ContainerWrapperIterator<const WrapperType, reference, reference>
Iterator;
100 return Iterator(WrapperType(
this),0);
106 return Iterator(WrapperType(
this),n);
113 return Iterator(WrapperType(
this),n-1);
120 return Iterator(WrapperType(
this),-1);
125 typedef ContainerWrapperIterator<const WrapperType, const_reference, const_reference>
ConstIterator;
140 ConstIterator
end ()
const 207 return p_==other.
scalar();
213 return p_!=other.
scalar();
219 template<
class X,
class Y>
220 void mv (
const X& x, Y& y)
const 222 #ifdef DUNE_FMatrix_WITH_CHECKING 223 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
224 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
226 for (size_type i=0; i<n; ++i)
231 template<
class X,
class Y>
232 void mtv (
const X& x, Y& y)
const 238 template<
class X,
class Y>
239 void umv (
const X& x, Y& y)
const 241 #ifdef DUNE_FMatrix_WITH_CHECKING 242 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
243 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
245 for (size_type i=0; i<n; ++i)
250 template<
class X,
class Y>
251 void umtv (
const X& x, Y& y)
const 253 #ifdef DUNE_FMatrix_WITH_CHECKING 254 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
255 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
257 for (size_type i=0; i<n; ++i)
262 template<
class X,
class Y>
263 void umhv (
const X& x, Y& y)
const 265 #ifdef DUNE_FMatrix_WITH_CHECKING 266 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
267 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
269 for (size_type i=0; i<n; i++)
270 y[i] += conjugateComplex(p_)*x[i];
274 template<
class X,
class Y>
275 void mmv (
const X& x, Y& y)
const 277 #ifdef DUNE_FMatrix_WITH_CHECKING 278 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
279 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
281 for (size_type i=0; i<n; ++i)
286 template<
class X,
class Y>
287 void mmtv (
const X& x, Y& y)
const 289 #ifdef DUNE_FMatrix_WITH_CHECKING 290 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
291 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
293 for (size_type i=0; i<n; ++i)
298 template<
class X,
class Y>
299 void mmhv (
const X& x, Y& y)
const 301 #ifdef DUNE_FMatrix_WITH_CHECKING 302 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
303 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
305 for (size_type i=0; i<n; i++)
306 y[i] -= conjugateComplex(p_)*x[i];
310 template<
class X,
class Y>
311 void usmv (
const K& alpha,
const X& x, Y& y)
const 313 #ifdef DUNE_FMatrix_WITH_CHECKING 314 if (x.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
315 if (y.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
317 for (size_type i=0; i<n; i++)
318 y[i] += alpha * p_ * x[i];
322 template<
class X,
class Y>
323 void usmtv (
const K& alpha,
const X& x, Y& y)
const 325 #ifdef DUNE_FMatrix_WITH_CHECKING 326 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
327 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
329 for (size_type i=0; i<n; i++)
330 y[i] += alpha * p_ * x[i];
334 template<
class X,
class Y>
335 void usmhv (
const K& alpha,
const X& x, Y& y)
const 337 #ifdef DUNE_FMatrix_WITH_CHECKING 338 if (x.N()!=
N()) DUNE_THROW(FMatrixError,
"index out of range");
339 if (y.N()!=
M()) DUNE_THROW(FMatrixError,
"index out of range");
341 for (size_type i=0; i<n; i++)
342 y[i] += alpha * conjugateComplex(p_) * x[i];
350 return fvmeta::sqrt(n*p_*p_);
368 return fvmeta::absreal(p_);
378 for (
int i=0; i<n; i++)
391 return std::pow(p_,n);
411 bool exists (size_type i, size_type j)
const 413 #ifdef DUNE_FMatrix_WITH_CHECKING 414 if (i<0 || i>=n) DUNE_THROW(FMatrixError,
"row index out of range");
415 if (j<0 || j>=n) DUNE_THROW(FMatrixError,
"column index out of range");
423 friend std::ostream& operator<< (std::ostream& s, const ScaledIdentityMatrix<K,n>& a)
425 for (size_type i=0; i<n; i++) {
426 for (size_type j=0; j<n; j++)
427 s << ((i==j) ? a.p_ : 0) <<
" ";
436 return reference(const_cast<K*>(&p_), i);
477 template <
class DenseMatrix,
class field,
int N>
479 static void apply(DenseMatrix& denseMatrix,
481 assert(denseMatrix.M() ==
N);
482 assert(denseMatrix.N() ==
N);
483 denseMatrix = field(0);
484 for (
int i = 0; i <
N; ++i)
485 denseMatrix[i][i] = rhs.
scalar();
Iterator begin()
begin iterator
Definition: scaledidmatrix.hh:98
const K & diagonal(size_type) const
Get const reference to diagonal entry.
Definition: scaledidmatrix.hh:446
void usmv(const K &alpha, const X &x, Y &y) const
y += alpha A x
Definition: scaledidmatrix.hh:311
reference operator[](size_type i)
Return reference object as row replacement.
Definition: scaledidmatrix.hh:434
Definition: allocator.hh:7
The number of block levels we contain. This is 1.
Definition: scaledidmatrix.hh:46
size_type M() const
number of blocks in column direction
Definition: scaledidmatrix.hh:403
void solve(V &x, const V &b) const
Solve system A x = b.
Definition: scaledidmatrix.hh:376
ContainerWrapperIterator< const WrapperType, reference, reference > Iterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:89
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:129
void umv(const X &x, Y &y) const
y += A x
Definition: scaledidmatrix.hh:239
ConstIterator beforeBegin() const
Definition: scaledidmatrix.hh:154
ScaledIdentityMatrix & operator*=(const K &k)
vector space multiplication with scalar
Definition: scaledidmatrix.hh:189
K & scalar()
Get reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:466
static void apply(DenseMatrix &denseMatrix, ScaledIdentityMatrix< field, N > const &rhs)
Definition: scaledidmatrix.hh:479
void mmv(const X &x, Y &y) const
y -= A x
Definition: scaledidmatrix.hh:275
ScaledIdentityMatrix(const K &k)
Constructor initializing the whole matrix with a scalar.
Definition: scaledidmatrix.hh:70
DiagonalRowVectorConst< K, n > const_row_type
Definition: scaledidmatrix.hh:52
row_type::Iterator ColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:95
Iterator iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:91
Iterator beforeBegin()
Definition: scaledidmatrix.hh:118
ConstIterator beforeEnd() const
Definition: scaledidmatrix.hh:147
const_row_type const_reference
Definition: scaledidmatrix.hh:53
Iterator beforeEnd()
Definition: scaledidmatrix.hh:111
ScaledIdentityMatrix & operator+=(const ScaledIdentityMatrix &y)
vector space addition
Definition: scaledidmatrix.hh:162
K field_type
export the type representing the field
Definition: scaledidmatrix.hh:35
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: scaledidmatrix.hh:348
Iterator end()
end iterator
Definition: scaledidmatrix.hh:104
The number of columns.
Definition: scaledidmatrix.hh:60
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: scaledidmatrix.hh:411
K determinant() const
calculates the determinant of this matrix
Definition: scaledidmatrix.hh:390
ConstIterator begin() const
begin iterator
Definition: scaledidmatrix.hh:134
A multiple of the identity matrix of static size.
Definition: scaledidmatrix.hh:27
const K & scalar() const
Get const reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:459
ScaledIdentityMatrix & operator/=(const K &k)
vector space division by scalar
Definition: scaledidmatrix.hh:196
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: scaledidmatrix.hh:299
bool identical(const ScaledIdentityMatrix< K, n > &other) const
Definition: scaledidmatrix.hh:82
void usmhv(const K &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: scaledidmatrix.hh:335
The number of rows.
Definition: scaledidmatrix.hh:58
K & diagonal(size_type)
Get reference to diagonal entry.
Definition: scaledidmatrix.hh:452
const_row_type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:131
void invert()
Compute inverse.
Definition: scaledidmatrix.hh:384
Iterator RowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:93
ScaledIdentityMatrix & operator=(const K &k)
Definition: scaledidmatrix.hh:75
const_reference operator[](size_type i) const
Return const_reference object as row replacement.
Definition: scaledidmatrix.hh:440
K block_type
export the type representing the components
Definition: scaledidmatrix.hh:38
ScaledIdentityMatrix()
Default constructor.
Definition: scaledidmatrix.hh:66
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: scaledidmatrix.hh:287
ContainerWrapperIterator< const WrapperType, const_reference, const_reference > ConstIterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:125
void mv(const X &x, Y &y) const
y = A x
Definition: scaledidmatrix.hh:220
ConstIterator const_iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:127
row_type reference
Definition: scaledidmatrix.hh:51
void umhv(const X &x, Y &y) const
y += A^H x
Definition: scaledidmatrix.hh:263
void mtv(const X &x, Y &y) const
y = A^T x
Definition: scaledidmatrix.hh:232
std::size_t size_type
The type used for the index access and size operations.
Definition: scaledidmatrix.hh:41
ConstIterator end() const
end iterator
Definition: scaledidmatrix.hh:140
bool operator==(const ScaledIdentityMatrix &other) const
comparison operator
Definition: scaledidmatrix.hh:205
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: scaledidmatrix.hh:354
DiagonalRowVector< K, n > row_type
Each row is implemented by a field vector.
Definition: scaledidmatrix.hh:50
void usmtv(const K &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: scaledidmatrix.hh:323
size_type N() const
number of blocks in row direction
Definition: scaledidmatrix.hh:397
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: scaledidmatrix.hh:366
ScaledIdentityMatrix & operator-=(const ScaledIdentityMatrix &y)
vector space subtraction
Definition: scaledidmatrix.hh:169
FieldTraits< field_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: scaledidmatrix.hh:360
void umtv(const X &x, Y &y) const
y += A^T x
Definition: scaledidmatrix.hh:251
bool operator!=(const ScaledIdentityMatrix &other) const
incomparison operator
Definition: scaledidmatrix.hh:211