Dune Core Modules (2.4.1)

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>
20
21namespace 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 {
46 blocklevel = 1
47 };
48
50 typedef DiagonalRowVector<K,n> row_type;
51 typedef row_type reference;
52 typedef DiagonalRowVectorConst<K,n> const_row_type;
53 typedef const_row_type const_reference;
54
56 enum {
58 rows = n,
60 cols = n
61 };
62
63 //===== constructors
67
71 : p_(k)
72 {}
73
74 //===== assignment from scalar
75 ScaledIdentityMatrix& operator= (const K& k)
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
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
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
440 const_reference operator[](size_type i) const
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
453 {
454 return p_;
455 }
456
459 const K& scalar() const
460 {
461 return p_;
462 }
463
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 M, class K, int n>
478 void istl_assign_to_fmatrix(DenseMatrix<M>& fm, const ScaledIdentityMatrix<K,n>& s)
479 {
480 fm = K();
481 for(int i=0; i<n; ++i)
482 fm[i][i] = s.scalar();
483 }
484
485} // end namespace
486
487#endif
Iterator class for sparse vector-like containers.
Definition: diagonalmatrix.hh:977
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:171
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:335
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: scaledidmatrix.hh:287
ScaledIdentityMatrix & operator-=(const ScaledIdentityMatrix &y)
vector space subtraction
Definition: scaledidmatrix.hh:169
const_row_type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:131
ConstIterator end() const
end iterator
Definition: scaledidmatrix.hh:140
Iterator beforeBegin()
Definition: scaledidmatrix.hh:118
bool operator!=(const ScaledIdentityMatrix &other) const
incomparison operator
Definition: scaledidmatrix.hh:211
void mmv(const X &x, Y &y) const
y -= A x
Definition: scaledidmatrix.hh:275
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:311
row_type::Iterator ColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:95
void mv(const X &x, Y &y) const
y = A x
Definition: scaledidmatrix.hh:220
void umtv(const X &x, Y &y) const
y += A^T x
Definition: scaledidmatrix.hh:251
void umhv(const X &x, Y &y) const
y += A^H x
Definition: scaledidmatrix.hh:263
DiagonalRowVector< K, n > row_type
Each row is implemented by a field vector.
Definition: scaledidmatrix.hh:50
ContainerWrapperIterator< const WrapperType, reference, reference > Iterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:89
Iterator beforeEnd()
Definition: scaledidmatrix.hh:111
K determinant() const
calculates the determinant of this matrix
Definition: scaledidmatrix.hh:390
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:323
Iterator end()
end iterator
Definition: scaledidmatrix.hh:104
Iterator iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:91
const K & scalar() const
Get const reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:459
void umv(const X &x, Y &y) const
y += A x
Definition: scaledidmatrix.hh:239
const K & diagonal(size_type) const
Get const reference to diagonal entry.
Definition: scaledidmatrix.hh:446
ContainerWrapperIterator< const WrapperType, const_reference, const_reference > ConstIterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:125
K & diagonal(size_type)
Get reference to diagonal entry.
Definition: scaledidmatrix.hh:452
@ rows
The number of rows.
Definition: scaledidmatrix.hh:58
@ cols
The number of columns.
Definition: scaledidmatrix.hh:60
void solve(V &x, const V &b) const
Solve system A x = b.
Definition: scaledidmatrix.hh:376
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: scaledidmatrix.hh:411
Iterator RowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:93
ConstIterator const_iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:127
ScaledIdentityMatrix()
Default constructor.
Definition: scaledidmatrix.hh:66
bool operator==(const ScaledIdentityMatrix &other) const
comparison operator
Definition: scaledidmatrix.hh:205
ConstIterator beforeBegin() const
Definition: scaledidmatrix.hh:154
ScaledIdentityMatrix & operator/=(const K &k)
vector space division by scalar
Definition: scaledidmatrix.hh:196
friend std::ostream & operator<<(std::ostream &s, const ScaledIdentityMatrix< K, n > &a)
Sends the matrix to an output stream.
Definition: scaledidmatrix.hh:423
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: scaledidmatrix.hh:354
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: scaledidmatrix.hh:348
FieldTraits< field_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: scaledidmatrix.hh:360
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:129
ConstIterator beforeEnd() const
Definition: scaledidmatrix.hh:147
size_type M() const
number of blocks in column direction
Definition: scaledidmatrix.hh:403
const_reference operator[](size_type i) const
Return const_reference object as row replacement.
Definition: scaledidmatrix.hh:440
ScaledIdentityMatrix(const K &k)
Constructor initializing the whole matrix with a scalar.
Definition: scaledidmatrix.hh:70
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 K &k)
vector space multiplication with scalar
Definition: scaledidmatrix.hh:189
void invert()
Compute inverse.
Definition: scaledidmatrix.hh:384
Iterator begin()
begin iterator
Definition: scaledidmatrix.hh:98
K & scalar()
Get reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:466
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:299
ScaledIdentityMatrix & operator+=(const ScaledIdentityMatrix &y)
vector space addition
Definition: scaledidmatrix.hh:162
void mtv(const X &x, Y &y) const
y = A^T x
Definition: scaledidmatrix.hh:232
reference operator[](size_type i)
Return reference object as row replacement.
Definition: scaledidmatrix.hh:434
ConstIterator begin() const
begin iterator
Definition: scaledidmatrix.hh:134
@ blocklevel
The number of block levels we contain. This is 1.
Definition: scaledidmatrix.hh:46
size_type N() const
number of blocks in row direction
Definition: scaledidmatrix.hh:397
This file implements a quadratic diagonal matrix of fixed size.
A few common exception classes.
Implements a matrix constructed from a given type representing a field and compile-time given number ...
Type traits to determine the type of reals (when working with complex numbers)
#define DUNE_THROW(E, m)
Definition: exceptions.hh:243
Dune namespace.
Definition: alignment.hh:10
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:103
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Dec 22, 23:30, 2024)