Dune Core Modules (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>
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 [[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;
49 typedef row_type reference;
50 typedef DiagonalRowVectorConst<K,n> const_row_type;
51 typedef const_row_type const_reference;
52
54 enum {
56 rows = n,
58 cols = n
59 };
60
61 //===== constructors
65
69 : p_(k)
70 {}
71
72 //===== assignment from scalar
73 ScaledIdentityMatrix& operator= (const K& k)
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
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
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
438 const_reference operator[](size_type i) const
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
451 {
452 return p_;
453 }
454
457 const K& scalar() const
458 {
459 return p_;
460 }
461
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
Iterator class for sparse vector-like containers.
Definition: diagonalmatrix.hh:977
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:151
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
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
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
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
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
@ rows
The number of rows.
Definition: scaledidmatrix.hh:56
@ cols
The number of columns.
Definition: scaledidmatrix.hh:58
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
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:216
Dune namespace.
Definition: alignedallocator.hh:11
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:161
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)