dune-istl  2.8.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  [[deprecated("Use free function blockLevel(). Will be removed after 2.8.")]]
45  static constexpr std::size_t blocklevel = 1;
46 
48  typedef DiagonalRowVector<K,n> row_type;
50  typedef DiagonalRowVectorConst<K,n> const_row_type;
52 
54  enum {
56  rows = n,
58  cols = n
59  };
60 
61  //===== constructors
65 
68  ScaledIdentityMatrix (const K& k)
69  : p_(k)
70  {}
71 
72  //===== assignment from scalar
74  {
75  p_ = k;
76  return *this;
77  }
78 
79  // check if matrix is identical to other matrix (not only identical values)
80  bool identical(const ScaledIdentityMatrix<K,n>& other) const
81  {
82  return (this==&other);
83  }
84 
85  //===== iterator interface to rows of the matrix
87  typedef ContainerWrapperIterator<const WrapperType, reference, reference> Iterator;
89  typedef Iterator iterator;
93  typedef typename row_type::Iterator ColIterator;
94 
97  {
98  return Iterator(WrapperType(this),0);
99  }
100 
103  {
104  return Iterator(WrapperType(this),n);
105  }
106 
110  {
111  return Iterator(WrapperType(this),n-1);
112  }
113 
117  {
118  return Iterator(WrapperType(this),-1);
119  }
120 
121 
123  typedef ContainerWrapperIterator<const WrapperType, const_reference, const_reference> ConstIterator;
129  typedef typename const_row_type::ConstIterator ConstColIterator;
130 
133  {
134  return ConstIterator(WrapperType(this),0);
135  }
136 
139  {
140  return ConstIterator(WrapperType(this),n);
141  }
142 
146  {
147  return ConstIterator(WrapperType(this),n-1);
148  }
149 
153  {
154  return ConstIterator(WrapperType(this),-1);
155  }
156 
157  //===== vector space arithmetic
158 
161  {
162  p_ += y.p_;
163  return *this;
164  }
165 
168  {
169  p_ -= y.p_;
170  return *this;
171  }
172 
175  {
176  p_ += k;
177  return *this;
178  }
179 
182  {
183  p_ -= k;
184  return *this;
185  }
188  {
189  p_ *= k;
190  return *this;
191  }
192 
195  {
196  p_ /= k;
197  return *this;
198  }
199 
200  //===== comparison ops
201 
203  bool operator==(const ScaledIdentityMatrix& other) const
204  {
205  return p_==other.scalar();
206  }
207 
209  bool operator!=(const ScaledIdentityMatrix& other) const
210  {
211  return p_!=other.scalar();
212  }
213 
214  //===== linear maps
215 
217  template<class X, class Y>
218  void mv (const X& x, Y& y) const
219  {
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");
223 #endif
224  for (size_type i=0; i<n; ++i)
225  y[i] = p_ * x[i];
226  }
227 
229  template<class X, class Y>
230  void mtv (const X& x, Y& y) const
231  {
232  mv(x, y);
233  }
234 
236  template<class X, class Y>
237  void umv (const X& x, Y& y) const
238  {
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");
242 #endif
243  for (size_type i=0; i<n; ++i)
244  y[i] += p_ * x[i];
245  }
246 
248  template<class X, class Y>
249  void umtv (const X& x, Y& y) const
250  {
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");
254 #endif
255  for (size_type i=0; i<n; ++i)
256  y[i] += p_ * x[i];
257  }
258 
260  template<class X, class Y>
261  void umhv (const X& x, Y& y) const
262  {
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");
266 #endif
267  for (size_type i=0; i<n; i++)
268  y[i] += conjugateComplex(p_)*x[i];
269  }
270 
272  template<class X, class Y>
273  void mmv (const X& x, Y& y) const
274  {
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");
278 #endif
279  for (size_type i=0; i<n; ++i)
280  y[i] -= p_ * x[i];
281  }
282 
284  template<class X, class Y>
285  void mmtv (const X& x, Y& y) const
286  {
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");
290 #endif
291  for (size_type i=0; i<n; ++i)
292  y[i] -= p_ * x[i];
293  }
294 
296  template<class X, class Y>
297  void mmhv (const X& x, Y& y) const
298  {
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");
302 #endif
303  for (size_type i=0; i<n; i++)
304  y[i] -= conjugateComplex(p_)*x[i];
305  }
306 
308  template<class X, class Y>
309  void usmv (const K& alpha, const X& x, Y& y) const
310  {
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");
314 #endif
315  for (size_type i=0; i<n; i++)
316  y[i] += alpha * p_ * x[i];
317  }
318 
320  template<class X, class Y>
321  void usmtv (const K& alpha, const X& x, Y& y) const
322  {
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");
326 #endif
327  for (size_type i=0; i<n; i++)
328  y[i] += alpha * p_ * x[i];
329  }
330 
332  template<class X, class Y>
333  void usmhv (const K& alpha, const X& x, Y& y) const
334  {
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");
338 #endif
339  for (size_type i=0; i<n; i++)
340  y[i] += alpha * conjugateComplex(p_) * x[i];
341  }
342 
343  //===== norms
344 
346  typename FieldTraits<field_type>::real_type frobenius_norm () const
347  {
348  return fvmeta::sqrt(n*p_*p_);
349  }
350 
352  typename FieldTraits<field_type>::real_type frobenius_norm2 () const
353  {
354  return n*p_*p_;
355  }
356 
358  typename FieldTraits<field_type>::real_type infinity_norm () const
359  {
360  return std::abs(p_);
361  }
362 
364  typename FieldTraits<field_type>::real_type infinity_norm_real () const
365  {
366  return fvmeta::absreal(p_);
367  }
368 
369  //===== solve
370 
373  template<class V>
374  void solve (V& x, const V& b) const
375  {
376  for (int i=0; i<n; i++)
377  x[i] = b[i]/p_;
378  }
379 
382  void invert()
383  {
384  p_ = 1/p_;
385  }
386 
388  K determinant () const {
389  return std::pow(p_,n);
390  }
391 
392  //===== sizes
393 
395  size_type N () const
396  {
397  return n;
398  }
399 
401  size_type M () const
402  {
403  return n;
404  }
405 
406  //===== query
407 
409  bool exists (size_type i, size_type j) const
410  {
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");
414 #endif
415  return i==j;
416  }
417 
418  //===== conversion operator
419 
421  friend std::ostream& operator<< (std::ostream& s, const ScaledIdentityMatrix<K,n>& a)
422  {
423  for (size_type i=0; i<n; i++) {
424  for (size_type j=0; j<n; j++)
425  s << ((i==j) ? a.p_ : 0) << " ";
426  s << std::endl;
427  }
428  return s;
429  }
430 
433  {
434  return reference(const_cast<K*>(&p_), i);
435  }
436 
439  {
440  return const_reference(const_cast<K*>(&p_), i);
441  }
442 
444  const K& diagonal(size_type /*i*/) const
445  {
446  return p_;
447  }
448 
450  K& diagonal(size_type /*i*/)
451  {
452  return p_;
453  }
454 
457  const K& scalar() const
458  {
459  return p_;
460  }
461 
464  K& scalar()
465  {
466  return p_;
467  }
468 
469  private:
470  // the data, very simply a single number
471  K p_;
472 
473  };
474 
475  template <class DenseMatrix, class field, int N>
476  struct DenseMatrixAssigner<DenseMatrix, ScaledIdentityMatrix<field, N>> {
477  static void apply(DenseMatrix& denseMatrix,
478  ScaledIdentityMatrix<field, N> const& rhs) {
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();
484  }
485  };
486 } // end namespace
487 
488 #endif
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
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
K & diagonal(size_type)
Get reference to diagonal entry.
Definition: scaledidmatrix.hh:450
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
ScaledIdentityMatrix & operator+=(const ScaledIdentityMatrix &y)
vector space addition
Definition: scaledidmatrix.hh:160
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
@ rows
The number of rows.
Definition: scaledidmatrix.hh:56
@ cols
The number of columns.
Definition: scaledidmatrix.hh:58
ContainerWrapperIterator< const WrapperType, const_reference, const_reference > ConstIterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:123
ScaledIdentityMatrix & operator/=(const K &k)
vector space division by scalar
Definition: scaledidmatrix.hh:194
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
FieldTraits< field_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: scaledidmatrix.hh:358
ConstIterator beforeBegin() const
Definition: scaledidmatrix.hh:152
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: scaledidmatrix.hh:364
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
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
ScaledIdentityMatrix & operator-=(const ScaledIdentityMatrix &y)
vector space subtraction
Definition: scaledidmatrix.hh:167
row_type reference
Definition: scaledidmatrix.hh:49
K block_type
export the type representing the components
Definition: scaledidmatrix.hh:38
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: scaledidmatrix.hh:346
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: scaledidmatrix.hh:297
DiagonalRowVectorConst< K, n > const_row_type
Definition: scaledidmatrix.hh:50
friend std::ostream & operator<<(std::ostream &s, const ScaledIdentityMatrix< K, n > &a)
Sends the matrix to an output stream.
Definition: scaledidmatrix.hh:421
void mtv(const X &x, Y &y) const
y = A^T x
Definition: scaledidmatrix.hh:230
ScaledIdentityMatrix & operator=(const K &k)
Definition: scaledidmatrix.hh:73
ScaledIdentityMatrix & operator*=(const K &k)
vector space multiplication with scalar
Definition: scaledidmatrix.hh:187
reference operator[](size_type i)
Return reference object as row replacement.
Definition: scaledidmatrix.hh:432
const K & diagonal(size_type) const
Get const reference to diagonal entry.
Definition: scaledidmatrix.hh:444
ConstIterator begin() const
begin iterator
Definition: scaledidmatrix.hh:132
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: scaledidmatrix.hh:352
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