dune-istl  2.7.0
scaledidmatrix.hh
Go to the documentation of this file.
1 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2 // vi: set et ts=4 sw=2 sts=2:
3 #ifndef DUNE_ISTL_SCALEDIDMATRIX_HH
4 #define DUNE_ISTL_SCALEDIDMATRIX_HH
5 
12 #include <cmath>
13 #include <cstddef>
14 #include <complex>
15 #include <iostream>
16 #include <dune/common/exceptions.hh>
17 #include <dune/common/fmatrix.hh>
18 #include <dune/common/diagonalmatrix.hh>
19 #include <dune/common/ftraits.hh>
20 
21 namespace Dune {
22 
26  template<class K, int n>
28  {
29  typedef DiagonalMatrixWrapper< ScaledIdentityMatrix<K,n> > WrapperType;
30 
31  public:
32  //===== type definitions and constants
33 
35  typedef K field_type;
36 
38  typedef K block_type;
39 
41  typedef std::size_t size_type;
42 
44  enum {
47  };
48 
50  typedef DiagonalRowVector<K,n> row_type;
52  typedef DiagonalRowVectorConst<K,n> const_row_type;
54 
56  enum {
58  rows = n,
60  cols = n
61  };
62 
63  //===== constructors
67 
70  ScaledIdentityMatrix (const K& k)
71  : p_(k)
72  {}
73 
74  //===== assignment from scalar
76  {
77  p_ = k;
78  return *this;
79  }
80 
81  // check if matrix is identical to other matrix (not only identical values)
82  bool identical(const ScaledIdentityMatrix<K,n>& other) const
83  {
84  return (this==&other);
85  }
86 
87  //===== iterator interface to rows of the matrix
89  typedef ContainerWrapperIterator<const WrapperType, reference, reference> Iterator;
91  typedef Iterator iterator;
95  typedef typename row_type::Iterator ColIterator;
96 
99  {
100  return Iterator(WrapperType(this),0);
101  }
102 
105  {
106  return Iterator(WrapperType(this),n);
107  }
108 
112  {
113  return Iterator(WrapperType(this),n-1);
114  }
115 
119  {
120  return Iterator(WrapperType(this),-1);
121  }
122 
123 
125  typedef ContainerWrapperIterator<const WrapperType, const_reference, const_reference> ConstIterator;
131  typedef typename const_row_type::ConstIterator ConstColIterator;
132 
135  {
136  return ConstIterator(WrapperType(this),0);
137  }
138 
141  {
142  return ConstIterator(WrapperType(this),n);
143  }
144 
148  {
149  return ConstIterator(WrapperType(this),n-1);
150  }
151 
155  {
156  return ConstIterator(WrapperType(this),-1);
157  }
158 
159  //===== vector space arithmetic
160 
163  {
164  p_ += y.p_;
165  return *this;
166  }
167 
170  {
171  p_ -= y.p_;
172  return *this;
173  }
174 
177  {
178  p_ += k;
179  return *this;
180  }
181 
184  {
185  p_ -= k;
186  return *this;
187  }
190  {
191  p_ *= k;
192  return *this;
193  }
194 
197  {
198  p_ /= k;
199  return *this;
200  }
201 
202  //===== comparison ops
203 
205  bool operator==(const ScaledIdentityMatrix& other) const
206  {
207  return p_==other.scalar();
208  }
209 
211  bool operator!=(const ScaledIdentityMatrix& other) const
212  {
213  return p_!=other.scalar();
214  }
215 
216  //===== linear maps
217 
219  template<class X, class Y>
220  void mv (const X& x, Y& y) const
221  {
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");
225 #endif
226  for (size_type i=0; i<n; ++i)
227  y[i] = p_ * x[i];
228  }
229 
231  template<class X, class Y>
232  void mtv (const X& x, Y& y) const
233  {
234  mv(x, y);
235  }
236 
238  template<class X, class Y>
239  void umv (const X& x, Y& y) const
240  {
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");
244 #endif
245  for (size_type i=0; i<n; ++i)
246  y[i] += p_ * x[i];
247  }
248 
250  template<class X, class Y>
251  void umtv (const X& x, Y& y) const
252  {
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");
256 #endif
257  for (size_type i=0; i<n; ++i)
258  y[i] += p_ * x[i];
259  }
260 
262  template<class X, class Y>
263  void umhv (const X& x, Y& y) const
264  {
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");
268 #endif
269  for (size_type i=0; i<n; i++)
270  y[i] += conjugateComplex(p_)*x[i];
271  }
272 
274  template<class X, class Y>
275  void mmv (const X& x, Y& y) const
276  {
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");
280 #endif
281  for (size_type i=0; i<n; ++i)
282  y[i] -= p_ * x[i];
283  }
284 
286  template<class X, class Y>
287  void mmtv (const X& x, Y& y) const
288  {
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");
292 #endif
293  for (size_type i=0; i<n; ++i)
294  y[i] -= p_ * x[i];
295  }
296 
298  template<class X, class Y>
299  void mmhv (const X& x, Y& y) const
300  {
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");
304 #endif
305  for (size_type i=0; i<n; i++)
306  y[i] -= conjugateComplex(p_)*x[i];
307  }
308 
310  template<class X, class Y>
311  void usmv (const K& alpha, const X& x, Y& y) const
312  {
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");
316 #endif
317  for (size_type i=0; i<n; i++)
318  y[i] += alpha * p_ * x[i];
319  }
320 
322  template<class X, class Y>
323  void usmtv (const K& alpha, const X& x, Y& y) const
324  {
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");
328 #endif
329  for (size_type i=0; i<n; i++)
330  y[i] += alpha * p_ * x[i];
331  }
332 
334  template<class X, class Y>
335  void usmhv (const K& alpha, const X& x, Y& y) const
336  {
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");
340 #endif
341  for (size_type i=0; i<n; i++)
342  y[i] += alpha * conjugateComplex(p_) * x[i];
343  }
344 
345  //===== norms
346 
348  typename FieldTraits<field_type>::real_type frobenius_norm () const
349  {
350  return fvmeta::sqrt(n*p_*p_);
351  }
352 
354  typename FieldTraits<field_type>::real_type frobenius_norm2 () const
355  {
356  return n*p_*p_;
357  }
358 
360  typename FieldTraits<field_type>::real_type infinity_norm () const
361  {
362  return std::abs(p_);
363  }
364 
366  typename FieldTraits<field_type>::real_type infinity_norm_real () const
367  {
368  return fvmeta::absreal(p_);
369  }
370 
371  //===== solve
372 
375  template<class V>
376  void solve (V& x, const V& b) const
377  {
378  for (int i=0; i<n; i++)
379  x[i] = b[i]/p_;
380  }
381 
384  void invert()
385  {
386  p_ = 1/p_;
387  }
388 
390  K determinant () const {
391  return std::pow(p_,n);
392  }
393 
394  //===== sizes
395 
397  size_type N () const
398  {
399  return n;
400  }
401 
403  size_type M () const
404  {
405  return n;
406  }
407 
408  //===== query
409 
411  bool exists (size_type i, size_type j) const
412  {
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");
416 #endif
417  return i==j;
418  }
419 
420  //===== conversion operator
421 
423  friend std::ostream& operator<< (std::ostream& s, const ScaledIdentityMatrix<K,n>& a)
424  {
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) << " ";
428  s << std::endl;
429  }
430  return s;
431  }
432 
435  {
436  return reference(const_cast<K*>(&p_), i);
437  }
438 
441  {
442  return const_reference(const_cast<K*>(&p_), i);
443  }
444 
446  const K& diagonal(size_type /*i*/) const
447  {
448  return p_;
449  }
450 
452  K& diagonal(size_type /*i*/)
453  {
454  return p_;
455  }
456 
459  const K& scalar() const
460  {
461  return p_;
462  }
463 
466  K& scalar()
467  {
468  return p_;
469  }
470 
471  private:
472  // the data, very simply a single number
473  K p_;
474 
475  };
476 
477  template <class DenseMatrix, class field, int N>
478  struct DenseMatrixAssigner<DenseMatrix, ScaledIdentityMatrix<field, N>> {
479  static void apply(DenseMatrix& denseMatrix,
480  ScaledIdentityMatrix<field, N> const& rhs) {
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();
486  }
487  };
488 } // end namespace
489 
490 #endif
Dune::ScaledIdentityMatrix::mmhv
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: scaledidmatrix.hh:299
Dune::ScaledIdentityMatrix::const_reference
const_row_type const_reference
Definition: scaledidmatrix.hh:53
Dune::ScaledIdentityMatrix::identical
bool identical(const ScaledIdentityMatrix< K, n > &other) const
Definition: scaledidmatrix.hh:82
Dune::ScaledIdentityMatrix::begin
ConstIterator begin() const
begin iterator
Definition: scaledidmatrix.hh:134
Dune::ScaledIdentityMatrix::row_type
DiagonalRowVector< K, n > row_type
Each row is implemented by a field vector.
Definition: scaledidmatrix.hh:50
Dune::ScaledIdentityMatrix::const_row_type
DiagonalRowVectorConst< K, n > const_row_type
Definition: scaledidmatrix.hh:52
Dune::ScaledIdentityMatrix::umtv
void umtv(const X &x, Y &y) const
y += A^T x
Definition: scaledidmatrix.hh:251
Dune::ScaledIdentityMatrix::beforeEnd
Iterator beforeEnd()
Definition: scaledidmatrix.hh:111
Dune::ScaledIdentityMatrix::operator+=
ScaledIdentityMatrix & operator+=(const ScaledIdentityMatrix &y)
vector space addition
Definition: scaledidmatrix.hh:162
Dune::ScaledIdentityMatrix::RowIterator
Iterator RowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:93
Dune::ScaledIdentityMatrix::mtv
void mtv(const X &x, Y &y) const
y = A^T x
Definition: scaledidmatrix.hh:232
Dune::ScaledIdentityMatrix::beforeBegin
ConstIterator beforeBegin() const
Definition: scaledidmatrix.hh:154
Dune::ScaledIdentityMatrix::usmv
void usmv(const K &alpha, const X &x, Y &y) const
y += alpha A x
Definition: scaledidmatrix.hh:311
Dune::ScaledIdentityMatrix::operator[]
reference operator[](size_type i)
Return reference object as row replacement.
Definition: scaledidmatrix.hh:434
Dune::ScaledIdentityMatrix::end
Iterator end()
end iterator
Definition: scaledidmatrix.hh:104
Dune::ScaledIdentityMatrix::ScaledIdentityMatrix
ScaledIdentityMatrix(const K &k)
Constructor initializing the whole matrix with a scalar.
Definition: scaledidmatrix.hh:70
Dune::ScaledIdentityMatrix::operator[]
const_reference operator[](size_type i) const
Return const_reference object as row replacement.
Definition: scaledidmatrix.hh:440
Dune::ScaledIdentityMatrix::beforeBegin
Iterator beforeBegin()
Definition: scaledidmatrix.hh:118
Dune::ScaledIdentityMatrix::ConstRowIterator
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:129
Dune::ScaledIdentityMatrix::operator=
ScaledIdentityMatrix & operator=(const K &k)
Definition: scaledidmatrix.hh:75
Dune::ScaledIdentityMatrix::mmtv
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: scaledidmatrix.hh:287
Dune::ScaledIdentityMatrix::operator*=
ScaledIdentityMatrix & operator*=(const K &k)
vector space multiplication with scalar
Definition: scaledidmatrix.hh:189
Dune::ScaledIdentityMatrix::scalar
const K & scalar() const
Get const reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:459
Dune::ScaledIdentityMatrix
A multiple of the identity matrix of static size.
Definition: scaledidmatrix.hh:27
Dune::ScaledIdentityMatrix::umv
void umv(const X &x, Y &y) const
y += A x
Definition: scaledidmatrix.hh:239
Dune::ScaledIdentityMatrix::ScaledIdentityMatrix
ScaledIdentityMatrix()
Default constructor.
Definition: scaledidmatrix.hh:66
Dune::ScaledIdentityMatrix::exists
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: scaledidmatrix.hh:411
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:423
Dune::ScaledIdentityMatrix::infinity_norm
FieldTraits< field_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: scaledidmatrix.hh:360
Dune::ScaledIdentityMatrix::umhv
void umhv(const X &x, Y &y) const
y += A^H x
Definition: scaledidmatrix.hh:263
Dune::ScaledIdentityMatrix::operator!=
bool operator!=(const ScaledIdentityMatrix &other) const
incomparison operator
Definition: scaledidmatrix.hh:211
Dune::ScaledIdentityMatrix::cols
@ cols
The number of columns.
Definition: scaledidmatrix.hh:60
Dune::ScaledIdentityMatrix::N
size_type N() const
number of blocks in row direction
Definition: scaledidmatrix.hh:397
Dune::ScaledIdentityMatrix::ColIterator
row_type::Iterator ColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:95
Dune::ScaledIdentityMatrix::frobenius_norm
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: scaledidmatrix.hh:348
Dune::ScaledIdentityMatrix::mv
void mv(const X &x, Y &y) const
y = A x
Definition: scaledidmatrix.hh:220
Dune::ScaledIdentityMatrix::beforeEnd
ConstIterator beforeEnd() const
Definition: scaledidmatrix.hh:147
Dune::ScaledIdentityMatrix::operator/=
ScaledIdentityMatrix & operator/=(const K &k)
vector space division by scalar
Definition: scaledidmatrix.hh:196
Dune::ScaledIdentityMatrix::block_type
K block_type
export the type representing the components
Definition: scaledidmatrix.hh:38
Dune::ScaledIdentityMatrix::Iterator
ContainerWrapperIterator< const WrapperType, reference, reference > Iterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:89
Dune::ScaledIdentityMatrix::operator==
bool operator==(const ScaledIdentityMatrix &other) const
comparison operator
Definition: scaledidmatrix.hh:205
Dune::ScaledIdentityMatrix::determinant
K determinant() const
calculates the determinant of this matrix
Definition: scaledidmatrix.hh:390
Dune
Definition: allocator.hh:7
Dune::ScaledIdentityMatrix::ConstColIterator
const_row_type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:131
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:366
Dune::ScaledIdentityMatrix::field_type
K field_type
export the type representing the field
Definition: scaledidmatrix.hh:35
Dune::ScaledIdentityMatrix::reference
row_type reference
Definition: scaledidmatrix.hh:51
Dune::ScaledIdentityMatrix::frobenius_norm2
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: scaledidmatrix.hh:354
Dune::ScaledIdentityMatrix::begin
Iterator begin()
begin iterator
Definition: scaledidmatrix.hh:98
Dune::ScaledIdentityMatrix::blocklevel
@ blocklevel
The number of block levels we contain. This is 1.
Definition: scaledidmatrix.hh:46
Dune::ScaledIdentityMatrix::operator-=
ScaledIdentityMatrix & operator-=(const ScaledIdentityMatrix &y)
vector space subtraction
Definition: scaledidmatrix.hh:169
Dune::ScaledIdentityMatrix::iterator
Iterator iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:91
Dune::DenseMatrixAssigner< DenseMatrix, ScaledIdentityMatrix< field, N > >::apply
static void apply(DenseMatrix &denseMatrix, ScaledIdentityMatrix< field, N > const &rhs)
Definition: scaledidmatrix.hh:479
Dune::ScaledIdentityMatrix::usmtv
void usmtv(const K &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: scaledidmatrix.hh:323
Dune::ScaledIdentityMatrix::scalar
K & scalar()
Get reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:466
Dune::ScaledIdentityMatrix::solve
void solve(V &x, const V &b) const
Solve system A x = b.
Definition: scaledidmatrix.hh:376
Dune::ScaledIdentityMatrix::rows
@ rows
The number of rows.
Definition: scaledidmatrix.hh:58
Dune::ScaledIdentityMatrix::end
ConstIterator end() const
end iterator
Definition: scaledidmatrix.hh:140
Dune::ScaledIdentityMatrix::size_type
std::size_t size_type
The type used for the index access and size operations.
Definition: scaledidmatrix.hh:41
Dune::ScaledIdentityMatrix::usmhv
void usmhv(const K &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: scaledidmatrix.hh:335
Dune::ScaledIdentityMatrix::mmv
void mmv(const X &x, Y &y) const
y -= A x
Definition: scaledidmatrix.hh:275
Dune::ScaledIdentityMatrix::M
size_type M() const
number of blocks in column direction
Definition: scaledidmatrix.hh:403
Dune::ScaledIdentityMatrix::invert
void invert()
Compute inverse.
Definition: scaledidmatrix.hh:384
Dune::ScaledIdentityMatrix::ConstIterator
ContainerWrapperIterator< const WrapperType, const_reference, const_reference > ConstIterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:125
Dune::ScaledIdentityMatrix::const_iterator
ConstIterator const_iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:127
Dune::ScaledIdentityMatrix::diagonal
const K & diagonal(size_type) const
Get const reference to diagonal entry.
Definition: scaledidmatrix.hh:446
Dune::ScaledIdentityMatrix::diagonal
K & diagonal(size_type)
Get reference to diagonal entry.
Definition: scaledidmatrix.hh:452