Dune Core Modules (unstable)

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
40 template<class K, int n>
42 {
43 typedef DiagonalMatrixWrapper< ScaledIdentityMatrix<K,n> > WrapperType;
44
45 public:
46 //===== type definitions and constants
47
49 typedef K field_type;
50
52 typedef K block_type;
53
55 typedef std::size_t size_type;
56
65 typedef DiagonalRowVector<K,n> row_type;
66 typedef row_type reference;
67
69 typedef DiagonalRowVectorConst<K,n> const_row_type;
70 typedef const_row_type const_reference;
71
73 static constexpr std::integral_constant<size_type,n> rows = {};
74
76 static constexpr std::integral_constant<size_type,n> cols = {};
77
78 //===== constructors
82
86 : p_(k)
87 {}
88
92 {
93 p_ = k;
94 return *this;
95 }
96
101 bool identical(const ScaledIdentityMatrix<K,n>& other) const
102 {
103 return (this==&other);
104 }
105
106 //===== iterator interface to rows of the matrix
115
118 {
119 return Iterator(WrapperType(this),0);
120 }
121
124 {
125 return Iterator(WrapperType(this),n);
126 }
127
131 {
132 return Iterator(WrapperType(this),n-1);
133 }
134
138 {
139 return Iterator(WrapperType(this),-1);
140 }
141
142
151
154 {
155 return ConstIterator(WrapperType(this),0);
156 }
157
160 {
161 return ConstIterator(WrapperType(this),n);
162 }
163
167 {
168 return ConstIterator(WrapperType(this),n-1);
169 }
170
174 {
175 return ConstIterator(WrapperType(this),-1);
176 }
177
181
184 {
185 p_ += y.p_;
186 return *this;
187 }
188
191 {
192 p_ -= y.p_;
193 return *this;
194 }
195
201 {
202 p_ += k;
203 return *this;
204 }
205
211 {
212 p_ -= k;
213 return *this;
214 }
215
218 {
219 p_ *= k;
220 return *this;
221 }
222
225 {
226 p_ /= k;
227 return *this;
228 }
229
231 template <class Scalar,
232 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
233 friend auto operator* ( const ScaledIdentityMatrix& matrix, Scalar scalar)
234 {
236 }
237
239 template <class Scalar,
240 std::enable_if_t<IsNumber<Scalar>::value, int> = 0>
241 friend auto operator* (Scalar scalar, const ScaledIdentityMatrix& matrix)
242 {
244 }
245
247 template <class OtherScalar>
248 requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
249 friend auto& operator+= (FieldMatrix<OtherScalar,n,n>& fieldMatrix,
250 const ScaledIdentityMatrix& matrix)
251 {
252 for (int i=0; i<n; i++)
253 fieldMatrix[i][i] += matrix.p_;
254
255 return fieldMatrix;
256 }
257
259 template <class OtherScalar>
260 requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
261 friend auto operator+ (const FieldMatrix<OtherScalar,n,n>& fieldMatrix,
262 const ScaledIdentityMatrix& matrix)
263 {
265 Result result = fieldMatrix;
266 result += matrix;
267 return result;
268 }
269
271 template <class OtherScalar>
272 requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
273 friend auto operator+ (const ScaledIdentityMatrix& matrix,
274 const FieldMatrix<OtherScalar,n,n>& fieldMatrix)
275 {
276 return fieldMatrix + matrix;
277 }
278
280 template <class OtherScalar>
281 requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
282 friend auto operator+= (DiagonalMatrix<OtherScalar,n>& diagonalMatrix,
283 const ScaledIdentityMatrix& matrix)
284 {
285 for (std::size_t i=0; i<n; i++)
286 diagonalMatrix.diagonal(i) += matrix.p_;
287
288 return diagonalMatrix;
289 }
290
292 template <class OtherScalar>
293 requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
294 friend auto operator+ (const DiagonalMatrix<OtherScalar,n>& diagonalMatrix,
295 const ScaledIdentityMatrix& matrix)
296 {
298 Result result = diagonalMatrix;
299 result += matrix;
300 return result;
301 }
302
304 template <class OtherScalar>
305 requires requires(K k, OtherScalar otherScalar) { k + otherScalar; }
306 friend auto operator+ (const ScaledIdentityMatrix& matrix,
307 const DiagonalMatrix<OtherScalar,n>& diagonalMatrix)
308 {
309 return diagonalMatrix + matrix;
310 }
312
313 //===== comparison ops
314
316 bool operator==(const ScaledIdentityMatrix& other) const
317 {
318 return p_==other.scalar();
319 }
320
322 bool operator!=(const ScaledIdentityMatrix& other) const
323 {
324 return p_!=other.scalar();
325 }
326
330
332 template<class X, class Y>
333 void mv (const X& x, Y& y) const
334 {
335#ifdef DUNE_FMatrix_WITH_CHECKING
336 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
337 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
338#endif
339 for (size_type i=0; i<n; ++i)
340 y[i] = p_ * x[i];
341 }
342
344 template<class X, class Y>
345 void mtv (const X& x, Y& y) const
346 {
347 mv(x, y);
348 }
349
351 template<class X, class Y>
352 void umv (const X& x, Y& y) const
353 {
354#ifdef DUNE_FMatrix_WITH_CHECKING
355 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
356 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
357#endif
358 for (size_type i=0; i<n; ++i)
359 y[i] += p_ * x[i];
360 }
361
363 template<class X, class Y>
364 void umtv (const X& x, Y& y) const
365 {
366#ifdef DUNE_FMatrix_WITH_CHECKING
367 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
368 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
369#endif
370 for (size_type i=0; i<n; ++i)
371 y[i] += p_ * x[i];
372 }
373
375 template<class X, class Y>
376 void umhv (const X& x, Y& y) const
377 {
378#ifdef DUNE_FMatrix_WITH_CHECKING
379 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
380 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
381#endif
382 for (size_type i=0; i<n; i++)
383 y[i] += conjugateComplex(p_)*x[i];
384 }
385
387 template<class X, class Y>
388 void mmv (const X& x, Y& y) const
389 {
390#ifdef DUNE_FMatrix_WITH_CHECKING
391 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
392 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
393#endif
394 for (size_type i=0; i<n; ++i)
395 y[i] -= p_ * x[i];
396 }
397
399 template<class X, class Y>
400 void mmtv (const X& x, Y& y) const
401 {
402#ifdef DUNE_FMatrix_WITH_CHECKING
403 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
404 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
405#endif
406 for (size_type i=0; i<n; ++i)
407 y[i] -= p_ * x[i];
408 }
409
411 template<class X, class Y>
412 void mmhv (const X& x, Y& y) const
413 {
414#ifdef DUNE_FMatrix_WITH_CHECKING
415 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
416 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
417#endif
418 for (size_type i=0; i<n; i++)
419 y[i] -= conjugateComplex(p_)*x[i];
420 }
421
423 template<class X, class Y>
424 void usmv (const K& alpha, const X& x, Y& y) const
425 {
426#ifdef DUNE_FMatrix_WITH_CHECKING
427 if (x.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
428 if (y.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
429#endif
430 for (size_type i=0; i<n; i++)
431 y[i] += alpha * p_ * x[i];
432 }
433
435 template<class X, class Y>
436 void usmtv (const K& alpha, const X& x, Y& y) const
437 {
438#ifdef DUNE_FMatrix_WITH_CHECKING
439 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
440 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
441#endif
442 for (size_type i=0; i<n; i++)
443 y[i] += alpha * p_ * x[i];
444 }
445
447 template<class X, class Y>
448 void usmhv (const K& alpha, const X& x, Y& y) const
449 {
450#ifdef DUNE_FMatrix_WITH_CHECKING
451 if (x.N()!=N()) DUNE_THROW(FMatrixError,"index out of range");
452 if (y.N()!=M()) DUNE_THROW(FMatrixError,"index out of range");
453#endif
454 for (size_type i=0; i<n; i++)
455 y[i] += alpha * conjugateComplex(p_) * x[i];
456 }
457
459
460
464
466 typename FieldTraits<field_type>::real_type frobenius_norm () const
467 {
468 return fvmeta::sqrt(n*p_*p_);
469 }
470
472 typename FieldTraits<field_type>::real_type frobenius_norm2 () const
473 {
474 return n*p_*p_;
475 }
476
482 typename FieldTraits<field_type>::real_type infinity_norm () const
483 {
484 return std::abs(p_);
485 }
486
488 typename FieldTraits<field_type>::real_type infinity_norm_real () const
489 {
490 return fvmeta::absreal(p_);
491 }
492
494
499 template<class V>
500 void solve (V& x, const V& b) const
501 {
502 for (int i=0; i<n; i++)
503 x[i] = b[i]/p_;
504 }
505
508 void invert()
509 {
510 p_ = 1/p_;
511 }
512
514 K determinant () const {
515 return std::pow(p_,n);
516 }
517
518 //===== sizes
519
521 size_type N () const
522 {
523 return n;
524 }
525
527 size_type M () const
528 {
529 return n;
530 }
531
532 //===== query
533
535 bool exists (size_type i, size_type j) const
536 {
537#ifdef DUNE_FMatrix_WITH_CHECKING
538 if (i<0 || i>=n) DUNE_THROW(FMatrixError,"row index out of range");
539 if (j<0 || j>=n) DUNE_THROW(FMatrixError,"column index out of range");
540#endif
541 return i==j;
542 }
543
544 //===== conversion operator
545
547 friend std::ostream& operator<< (std::ostream& s, const ScaledIdentityMatrix<K,n>& a)
548 {
549 for (size_type i=0; i<n; i++) {
550 for (size_type j=0; j<n; j++)
551 s << ((i==j) ? a.p_ : 0) << " ";
552 s << std::endl;
553 }
554 return s;
555 }
556
559 {
560 return reference(const_cast<K*>(&p_), i);
561 }
562
564 const_reference operator[](size_type i) const
565 {
566 return const_reference(const_cast<K*>(&p_), i);
567 }
568
570 const K& diagonal(size_type /*i*/) const
571 {
572 return p_;
573 }
574
577 {
578 return p_;
579 }
580
583 const K& scalar() const
584 {
585 return p_;
586 }
587
591 {
592 return p_;
593 }
594
595 private:
596 // the data, very simply a single number
597 K p_;
598
599 };
600
601 template <class DenseMatrix, class field, int N>
602 struct DenseMatrixAssigner<DenseMatrix, ScaledIdentityMatrix<field, N>> {
603 static void apply(DenseMatrix& denseMatrix,
604 ScaledIdentityMatrix<field, N> const& rhs) {
605 assert(denseMatrix.M() == N);
606 assert(denseMatrix.N() == N);
607 denseMatrix = field(0);
608 for (int i = 0; i < N; ++i)
609 denseMatrix[i][i] = rhs.scalar();
610 }
611 };
612
613 template<class K, int n>
614 struct FieldTraits< ScaledIdentityMatrix<K, n> >
615 {
616 using field_type = typename ScaledIdentityMatrix<K, n>::field_type;
617 using real_type = typename FieldTraits<field_type>::real_type;
618 };
619
620} // end namespace
621
622#endif
Iterator class for sparse vector-like containers.
Definition: diagonalmatrix.hh:1078
A diagonal matrix of static size.
Definition: diagonalmatrix.hh:56
Error thrown if operations of a FieldMatrix fail.
Definition: densematrix.hh:131
A dense n x m matrix.
Definition: fmatrix.hh:117
A multiple of the identity matrix of static size.
Definition: scaledidmatrix.hh:42
void usmhv(const K &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: scaledidmatrix.hh:448
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: scaledidmatrix.hh:400
ScaledIdentityMatrix & operator-=(const ScaledIdentityMatrix &y)
Vector space subtraction.
Definition: scaledidmatrix.hh:190
const_row_type::ConstIterator ConstColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:150
ConstIterator end() const
end iterator
Definition: scaledidmatrix.hh:159
Iterator beforeBegin()
Definition: scaledidmatrix.hh:137
static constexpr std::integral_constant< size_type, n > rows
The number of matrix rows.
Definition: scaledidmatrix.hh:73
bool operator!=(const ScaledIdentityMatrix &other) const
Test for inequality.
Definition: scaledidmatrix.hh:322
void mmv(const X &x, Y &y) const
y -= A x
Definition: scaledidmatrix.hh:388
std::size_t size_type
The type used for the indices and sizes.
Definition: scaledidmatrix.hh:55
void usmv(const K &alpha, const X &x, Y &y) const
y += alpha A x
Definition: scaledidmatrix.hh:424
row_type::Iterator ColIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:114
void mv(const X &x, Y &y) const
y = A x
Definition: scaledidmatrix.hh:333
void umtv(const X &x, Y &y) const
y += A^T x
Definition: scaledidmatrix.hh:364
void umhv(const X &x, Y &y) const
y += A^H x
Definition: scaledidmatrix.hh:376
DiagonalRowVector< K, n > row_type
Type of a single matrix row.
Definition: scaledidmatrix.hh:65
ContainerWrapperIterator< const WrapperType, reference, reference > Iterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:108
Iterator beforeEnd()
Definition: scaledidmatrix.hh:130
K determinant() const
Calculates the determinant of this matrix.
Definition: scaledidmatrix.hh:514
K field_type
The type representing numbers.
Definition: scaledidmatrix.hh:49
void usmtv(const K &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: scaledidmatrix.hh:436
Iterator end()
end iterator
Definition: scaledidmatrix.hh:123
Iterator iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:110
static constexpr std::integral_constant< size_type, n > cols
The number of matrix columns.
Definition: scaledidmatrix.hh:76
const K & scalar() const
Get const reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:583
friend auto operator+(const FieldMatrix< OtherScalar, n, n > &fieldMatrix, const ScaledIdentityMatrix &matrix)
Addition of ScaledIdentityMatrix to FieldMatrix.
Definition: scaledidmatrix.hh:261
void umv(const X &x, Y &y) const
y += A x
Definition: scaledidmatrix.hh:352
const K & diagonal(size_type) const
Get const reference to diagonal entry.
Definition: scaledidmatrix.hh:570
ScaledIdentityMatrix & operator=(const K &k)
Assignment from scalar.
Definition: scaledidmatrix.hh:91
ContainerWrapperIterator< const WrapperType, const_reference, const_reference > ConstIterator
Iterator class for sequential access.
Definition: scaledidmatrix.hh:144
K & diagonal(size_type)
Get reference to diagonal entry.
Definition: scaledidmatrix.hh:576
void solve(V &x, const V &b) const
Solve linear system A x = b.
Definition: scaledidmatrix.hh:500
bool exists(size_type i, size_type j) const
Return true if (i,j) is in the matrix pattern, i.e., if i==j.
Definition: scaledidmatrix.hh:535
Iterator RowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:112
ConstIterator const_iterator
typedef for stl compliant access
Definition: scaledidmatrix.hh:146
ScaledIdentityMatrix()
Default constructor.
Definition: scaledidmatrix.hh:81
bool operator==(const ScaledIdentityMatrix &other) const
Test for equality.
Definition: scaledidmatrix.hh:316
ConstIterator beforeBegin() const
Definition: scaledidmatrix.hh:173
ScaledIdentityMatrix & operator/=(const K &k)
Vector space division by scalar.
Definition: scaledidmatrix.hh:224
friend std::ostream & operator<<(std::ostream &s, const ScaledIdentityMatrix< K, n > &a)
Sends the matrix to an output stream.
Definition: scaledidmatrix.hh:547
FieldTraits< field_type >::real_type frobenius_norm2() const
The square of the Frobenius norm.
Definition: scaledidmatrix.hh:472
FieldTraits< field_type >::real_type frobenius_norm() const
The Frobenius norm.
Definition: scaledidmatrix.hh:466
FieldTraits< field_type >::real_type infinity_norm() const
The row sum norm.
Definition: scaledidmatrix.hh:482
ConstIterator ConstRowIterator
rename the iterators for easier access
Definition: scaledidmatrix.hh:148
ConstIterator beforeEnd() const
Definition: scaledidmatrix.hh:166
size_type M() const
number of blocks in column direction
Definition: scaledidmatrix.hh:527
const_reference operator[](size_type i) const
Return const_reference object as row replacement.
Definition: scaledidmatrix.hh:564
ScaledIdentityMatrix(const K &k)
Constructor initializing the whole matrix with a scalar.
Definition: scaledidmatrix.hh:85
friend auto operator*(const ScaledIdentityMatrix &matrix, Scalar scalar)
Vector space multiplication with scalar.
Definition: scaledidmatrix.hh:233
FieldTraits< field_type >::real_type infinity_norm_real() const
simplified infinity norm (uses Manhattan norm for complex values)
Definition: scaledidmatrix.hh:488
ScaledIdentityMatrix & operator*=(const K &k)
Vector space multiplication with scalar.
Definition: scaledidmatrix.hh:217
bool identical(const ScaledIdentityMatrix< K, n > &other) const
Check if matrix is identical to other matrix.
Definition: scaledidmatrix.hh:101
void invert()
Compute inverse.
Definition: scaledidmatrix.hh:508
Iterator begin()
begin iterator
Definition: scaledidmatrix.hh:117
K & scalar()
Get reference to the scalar diagonal value.
Definition: scaledidmatrix.hh:590
K block_type
The type representing matrix entries (which may be matrices themselves)
Definition: scaledidmatrix.hh:52
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: scaledidmatrix.hh:412
ScaledIdentityMatrix & operator+=(const ScaledIdentityMatrix &y)
Vector space addition.
Definition: scaledidmatrix.hh:183
DiagonalRowVectorConst< K, n > const_row_type
Const type of a single row.
Definition: scaledidmatrix.hh:69
void mtv(const X &x, Y &y) const
y = A^T x
Definition: scaledidmatrix.hh:345
reference operator[](size_type i)
Return reference object as row replacement.
Definition: scaledidmatrix.hh:558
ConstIterator begin() const
begin iterator
Definition: scaledidmatrix.hh:153
size_type N() const
number of blocks in row direction
Definition: scaledidmatrix.hh:521
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)
const K & diagonal(size_type i) const
Get const reference to diagonal entry.
Definition: diagonalmatrix.hh:580
#define DUNE_THROW(E,...)
Definition: exceptions.hh:314
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
constexpr K conjugateComplex(const K &x)
compute conjugate complex of x
Definition: math.hh:146
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden & Uni Heidelberg  |  generated with Hugo v0.111.3 (Sep 21, 22:40, 2025)