DUNE PDELab (git)

scaledidmatrix.hh
Go to the documentation of this file.
1// SPDX-FileCopyrightText: Copyright © DUNE Project contributors, see file LICENSE.md in module root
2// SPDX-License-Identifier: LicenseRef-GPL-2.0-only-with-DUNE-exception
3// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
4// vi: set et ts=4 sw=2 sts=2:
5#ifndef DUNE_ISTL_SCALEDIDMATRIX_HH
6#define DUNE_ISTL_SCALEDIDMATRIX_HH
7
14#include <cmath>
15#include <cstddef>
16#include <complex>
17#include <iostream>
22
23namespace Dune {
24
28 template<class K, int n>
30 {
31 typedef DiagonalMatrixWrapper< ScaledIdentityMatrix<K,n> > WrapperType;
32
33 public:
34 //===== type definitions and constants
35
37 typedef K field_type;
38
40 typedef K block_type;
41
43 typedef std::size_t size_type;
44
46 typedef DiagonalRowVector<K,n> row_type;
47 typedef row_type reference;
48 typedef DiagonalRowVectorConst<K,n> const_row_type;
49 typedef const_row_type const_reference;
50
52 enum {
54 rows = n,
56 cols = n
57 };
58
59 //===== constructors
63
67 : p_(k)
68 {}
69
70 //===== assignment from scalar
71 ScaledIdentityMatrix& operator= (const K& k)
72 {
73 p_ = k;
74 return *this;
75 }
76
77 // check if matrix is identical to other matrix (not only identical values)
78 bool identical(const ScaledIdentityMatrix<K,n>& other) const
79 {
80 return (this==&other);
81 }
82
83 //===== iterator interface to rows of the matrix
92
95 {
96 return Iterator(WrapperType(this),0);
97 }
98
101 {
102 return Iterator(WrapperType(this),n);
103 }
104
108 {
109 return Iterator(WrapperType(this),n-1);
110 }
111
115 {
116 return Iterator(WrapperType(this),-1);
117 }
118
119
128
131 {
132 return ConstIterator(WrapperType(this),0);
133 }
134
137 {
138 return ConstIterator(WrapperType(this),n);
139 }
140
144 {
145 return ConstIterator(WrapperType(this),n-1);
146 }
147
151 {
152 return ConstIterator(WrapperType(this),-1);
153 }
154
155 //===== vector space arithmetic
156
159 {
160 p_ += y.p_;
161 return *this;
162 }
163
166 {
167 p_ -= y.p_;
168 return *this;
169 }
170
173 {
174 p_ += k;
175 return *this;
176 }
177
180 {
181 p_ -= k;
182 return *this;
183 }
186 {
187 p_ *= k;
188 return *this;
189 }
190
193 {
194 p_ /= k;
195 return *this;
196 }
197
198 //===== binary operators
199
201 template <class Scalar,
202 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
203 friend auto operator* ( const ScaledIdentityMatrix& matrix, Scalar scalar)
204 {
206 }
207
209 template <class Scalar,
210 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
211 friend auto operator* (Scalar scalar, const ScaledIdentityMatrix& matrix)
212 {
214 }
215
216 //===== comparison ops
217
219 bool operator==(const ScaledIdentityMatrix& other) const
220 {
221 return p_==other.scalar();
222 }
223
225 bool operator!=(const ScaledIdentityMatrix& other) const
226 {
227 return p_!=other.scalar();
228 }
229
230 //===== linear maps
231
233 template<class X, class Y>
234 void mv (const X& x, Y& y) const
235 {
236#ifdef DUNE_FMatrix_WITH_CHECKING
237 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
238 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
239#endif
240 for (size_type i=0; i<n; ++i)
241 y[i] = p_ * x[i];
242 }
243
245 template<class X, class Y>
246 void mtv (const X& x, Y& y) const
247 {
248 mv(x, y);
249 }
250
252 template<class X, class Y>
253 void umv (const X& x, Y& y) const
254 {
255#ifdef DUNE_FMatrix_WITH_CHECKING
256 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
257 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
258#endif
259 for (size_type i=0; i<n; ++i)
260 y[i] += p_ * x[i];
261 }
262
264 template<class X, class Y>
265 void umtv (const X& x, Y& y) const
266 {
267#ifdef DUNE_FMatrix_WITH_CHECKING
268 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
269 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
270#endif
271 for (size_type i=0; i<n; ++i)
272 y[i] += p_ * x[i];
273 }
274
276 template<class X, class Y>
277 void umhv (const X& x, Y& y) const
278 {
279#ifdef DUNE_FMatrix_WITH_CHECKING
280 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
281 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
282#endif
283 for (size_type i=0; i<n; i++)
284 y[i] += conjugateComplex(p_)*x[i];
285 }
286
288 template<class X, class Y>
289 void mmv (const X& x, Y& y) const
290 {
291#ifdef DUNE_FMatrix_WITH_CHECKING
292 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
293 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
294#endif
295 for (size_type i=0; i<n; ++i)
296 y[i] -= p_ * x[i];
297 }
298
300 template<class X, class Y>
301 void mmtv (const X& x, Y& y) const
302 {
303#ifdef DUNE_FMatrix_WITH_CHECKING
304 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
305 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
306#endif
307 for (size_type i=0; i<n; ++i)
308 y[i] -= p_ * x[i];
309 }
310
312 template<class X, class Y>
313 void mmhv (const X& x, Y& y) const
314 {
315#ifdef DUNE_FMatrix_WITH_CHECKING
316 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
317 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
318#endif
319 for (size_type i=0; i<n; i++)
320 y[i] -= conjugateComplex(p_)*x[i];
321 }
322
324 template<class X, class Y>
325 void usmv (const K& alpha, const X& x, Y& y) const
326 {
327#ifdef DUNE_FMatrix_WITH_CHECKING
328 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
329 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
330#endif
331 for (size_type i=0; i<n; i++)
332 y[i] += alpha * p_ * x[i];
333 }
334
336 template<class X, class Y>
337 void usmtv (const K& alpha, const X& x, Y& y) const
338 {
339#ifdef DUNE_FMatrix_WITH_CHECKING
340 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
341 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
342#endif
343 for (size_type i=0; i<n; i++)
344 y[i] += alpha * p_ * x[i];
345 }
346
348 template<class X, class Y>
349 void usmhv (const K& alpha, const X& x, Y& y) const
350 {
351#ifdef DUNE_FMatrix_WITH_CHECKING
352 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
353 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
354#endif
355 for (size_type i=0; i<n; i++)
356 y[i] += alpha * conjugateComplex(p_) * x[i];
357 }
358
359 //===== norms
360
362 typename FieldTraits<field_type>::real_type frobenius_norm () const
363 {
364 return fvmeta::sqrt(n*p_*p_);
365 }
366
368 typename FieldTraits<field_type>::real_type frobenius_norm2 () const
369 {
370 return n*p_*p_;
371 }
372
374 typename FieldTraits<field_type>::real_type infinity_norm () const
375 {
376 return std::abs(p_);
377 }
378
380 typename FieldTraits<field_type>::real_type infinity_norm_real () const
381 {
382 return fvmeta::absreal(p_);
383 }
384
385 //===== solve
386
389 template<class V>
390 void solve (V& x, const V& b) const
391 {
392 for (int i=0; i<n; i++)
393 x[i] = b[i]/p_;
394 }
395
398 void invert()
399 {
400 p_ = 1/p_;
401 }
402
404 K determinant () const {
405 return std::pow(p_,n);
406 }
407
408 //===== sizes
409
411 size_type N () const
412 {
413 return n;
414 }
415
417 size_type M () const
418 {
419 return n;
420 }
421
422 //===== query
423
425 bool exists (size_type i, size_type j) const
426 {
427#ifdef DUNE_FMatrix_WITH_CHECKING
428 if (i<0 || i>=n) DUNE_THROW(FMatrixError,"row index out of range");
429 if (j<0 || j>=n) DUNE_THROW(FMatrixError,"column index out of range");
430#endif
431 return i==j;
432 }
433
434 //===== conversion operator
435
437 friend std::ostream& operator<< (std::ostream& s, const ScaledIdentityMatrix<K,n>& a)
438 {
439 for (size_type i=0; i<n; i++) {
440 for (size_type j=0; j<n; j++)
441 s << ((i==j) ? a.p_ : 0) << " ";
442 s << std::endl;
443 }
444 return s;
445 }
446
449 {
450 return reference(const_cast<K*>(&p_), i);
451 }
452
454 const_reference operator[](size_type i) const
455 {
456 return const_reference(const_cast<K*>(&p_), i);
457 }
458
460 const K& diagonal(size_type /*i*/) const
461 {
462 return p_;
463 }
464
467 {
468 return p_;
469 }
470
473 const K& scalar() const
474 {
475 return p_;
476 }
477
481 {
482 return p_;
483 }
484
485 private:
486 // the data, very simply a single number
487 K p_;
488
489 };
490
491 template <class DenseMatrix, class field, int N>
492 struct DenseMatrixAssigner<DenseMatrix, ScaledIdentityMatrix<field, N>> {
493 static void apply(DenseMatrix& denseMatrix,
494 ScaledIdentityMatrix<field, N> const& rhs) {
495 assert(denseMatrix.M() == N);
496 assert(denseMatrix.N() == N);
497 denseMatrix = field(0);
498 for (int i = 0; i < N; ++i)
499 denseMatrix[i][i] = rhs.scalar();
500 }
501 };
502
503 template<class K, int n>
504 struct FieldTraits< ScaledIdentityMatrix<K, n> >
505 {
506 using field_type = typename ScaledIdentityMatrix<K, n>::field_type;
507 using real_type = typename FieldTraits<field_type>::real_type;
508 };
509
510} // end namespace
511
512#endif
Iterator class for sparse vector-like containers.
Definition: diagonalmatrix.hh:997
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:126
A multiple of the identity matrix of static size.
Definition: scaledidmatrix.hh:30
void usmhv(const K &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: scaledidmatrix.hh:349
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: scaledidmatrix.hh:301
ScaledIdentityMatrix & operator-=(const ScaledIdentityMatrix &y)
vector space subtraction
Definition: scaledidmatrix.hh:165
const_row_type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:127
ConstIterator end() const
end iterator
Definition: scaledidmatrix.hh:136
Iterator beforeBegin()
Definition: scaledidmatrix.hh:114
bool operator!=(const ScaledIdentityMatrix &other) const
incomparison operator
Definition: scaledidmatrix.hh:225
void mmv(const X &x, Y &y) const
y -= A x
Definition: scaledidmatrix.hh:289
std::size_t size_type
The type used for the index access and size operations.
Definition: scaledidmatrix.hh:43
void usmv(const K &alpha, const X &x, Y &y) const
y += alpha A x
Definition: scaledidmatrix.hh:325
row_type::Iterator ColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:91
void mv(const X &x, Y &y) const
y = A x
Definition: scaledidmatrix.hh:234
void umtv(const X &x, Y &y) const
y += A^T x
Definition: scaledidmatrix.hh:265
void umhv(const X &x, Y &y) const
y += A^H x
Definition: scaledidmatrix.hh:277
DiagonalRowVector< K, n > row_type
Each row is implemented by a field vector.
Definition: scaledidmatrix.hh:46
ContainerWrapperIterator< const WrapperType, reference, reference > Iterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:85
Iterator beforeEnd()
Definition: scaledidmatrix.hh:107
K determinant() const
calculates the determinant of this matrix
Definition: scaledidmatrix.hh:404
K field_type
export the type representing the field
Definition: scaledidmatrix.hh:37
void usmtv(const K &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: scaledidmatrix.hh:337
Iterator end()
end iterator
Definition: scaledidmatrix.hh:100
Iterator iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:87
const K & scalar() const
Get const reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:473
void umv(const X &x, Y &y) const
y += A x
Definition: scaledidmatrix.hh:253
@ rows
The number of rows.
Definition: scaledidmatrix.hh:54
@ cols
The number of columns.
Definition: scaledidmatrix.hh:56
const K & diagonal(size_type) const
Get const reference to diagonal entry.
Definition: scaledidmatrix.hh:460
ContainerWrapperIterator< const WrapperType, const_reference, const_reference > ConstIterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:121
K & diagonal(size_type)
Get reference to diagonal entry.
Definition: scaledidmatrix.hh:466
void solve(V &x, const V &b) const
Solve system A x = b.
Definition: scaledidmatrix.hh:390
bool exists(size_type i, size_type j) const
return true when (i,j) is in pattern
Definition: scaledidmatrix.hh:425
Iterator RowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:89
ConstIterator const_iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:123
ScaledIdentityMatrix()
Default constructor.
Definition: scaledidmatrix.hh:62
bool operator==(const ScaledIdentityMatrix &other) const
comparison operator
Definition: scaledidmatrix.hh:219
ConstIterator beforeBegin() const
Definition: scaledidmatrix.hh:150
ScaledIdentityMatrix & operator/=(const K &k)
vector space division by scalar
Definition: scaledidmatrix.hh:192
friend std::ostream & operator<<(std::ostream &s, const ScaledIdentityMatrix< K, n > &a)
Sends the matrix to an output stream.
Definition: scaledidmatrix.hh:437
FieldTraits< field_type >::real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: scaledidmatrix.hh:368
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: scaledidmatrix.hh:362
FieldTraits< field_type >::real_type infinity_norm() const
infinity norm (row sum norm, how to generalize for blocks?)
Definition: scaledidmatrix.hh:374
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:125
ConstIterator beforeEnd() const
Definition: scaledidmatrix.hh:143
size_type M() const
number of blocks in column direction
Definition: scaledidmatrix.hh:417
const_reference operator[](size_type i) const
Return const_reference object as row replacement.
Definition: scaledidmatrix.hh:454
ScaledIdentityMatrix(const K &k)
Constructor initializing the whole matrix with a scalar.
Definition: scaledidmatrix.hh:66
friend auto operator*(const ScaledIdentityMatrix &matrix, Scalar scalar)
vector space multiplication with scalar
Definition: scaledidmatrix.hh:203
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: scaledidmatrix.hh:380
ScaledIdentityMatrix & operator*=(const K &k)
vector space multiplication with scalar
Definition: scaledidmatrix.hh:185
void invert()
Compute inverse.
Definition: scaledidmatrix.hh:398
Iterator begin()
begin iterator
Definition: scaledidmatrix.hh:94
K & scalar()
Get reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:480
K block_type
export the type representing the components
Definition: scaledidmatrix.hh:40
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: scaledidmatrix.hh:313
ScaledIdentityMatrix & operator+=(const ScaledIdentityMatrix &y)
vector space addition
Definition: scaledidmatrix.hh:158
void mtv(const X &x, Y &y) const
y = A^T x
Definition: scaledidmatrix.hh:246
reference operator[](size_type i)
Return reference object as row replacement.
Definition: scaledidmatrix.hh:448
ConstIterator begin() const
begin iterator
Definition: scaledidmatrix.hh:130
size_type N() const
number of blocks in row direction
Definition: scaledidmatrix.hh:411
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:218
typename Overloads::ScalarType< std::decay_t< V > >::type Scalar
Element type of some SIMD type.
Definition: interface.hh:235
Dune namespace.
Definition: alignedallocator.hh:13
K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:164
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)