Dune Core Modules (2.8.0)

multitypeblockmatrix.hh
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_MULTITYPEBLOCKMATRIX_HH
4#define DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
5
6#include <cmath>
7#include <iostream>
8#include <tuple>
9
10#include <dune/common/hybridutilities.hh>
11
12#include "istlexception.hh"
13
14// forward declaration
15namespace Dune
16{
17 template<typename FirstRow, typename... Args>
18 class MultiTypeBlockMatrix;
19
20 template<int I, int crow, int remain_row>
21 class MultiTypeBlockMatrix_Solver;
22}
23
24#include "gsetc.hh"
25
26namespace Dune {
27
41 template<typename FirstRow, typename... Args>
43 : public std::tuple<FirstRow, Args...>
44 {
45 using ParentType = std::tuple<FirstRow, Args...>;
46 public:
47
49 using ParentType::ParentType;
50
54 typedef MultiTypeBlockMatrix<FirstRow, Args...> type;
55
57 using size_type = std::size_t;
58
59 typedef typename FirstRow::field_type field_type;
60
62 static constexpr size_type N()
63 {
64 return 1+sizeof...(Args);
65 }
66
72 [[deprecated("Use method 'N' instead")]]
73 static constexpr size_type size()
74 {
75 return 1+sizeof...(Args);
76 }
77
79 static constexpr size_type M()
80 {
81 return FirstRow::size();
82 }
83
100 template< size_type index >
101 auto
102 operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable)
103 -> decltype(std::get<index>(*this))
104 {
105 return std::get<index>(*this);
106 }
107
113 template< size_type index >
114 auto
115 operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable) const
116 -> decltype(std::get<index>(*this))
117 {
118 return std::get<index>(*this);
119 }
120
124 template<typename T>
125 void operator= (const T& newval) {
126 using namespace Dune::Hybrid;
127 auto size = index_constant<1+sizeof...(Args)>();
128 // Since Dune::Hybrid::size(MultiTypeBlockMatrix) is not implemented,
129 // we cannot use a plain forEach(*this, ...). This could be achieved,
130 // e.g., by implementing a static size() method.
131 forEach(integralRange(size), [&](auto&& i) {
132 (*this)[i] = newval;
133 });
134 }
135
136 //===== vector space arithmetic
137
139 MultiTypeBlockMatrix& operator*= (const field_type& k)
140 {
141 auto size = index_constant<N()>();
143 (*this)[i] *= k;
144 });
145
146 return *this;
147 }
148
150 MultiTypeBlockMatrix& operator/= (const field_type& k)
151 {
152 auto size = index_constant<N()>();
154 (*this)[i] /= k;
155 });
156
157 return *this;
158 }
159
160
167 {
168 auto size = index_constant<N()>();
170 (*this)[i] += b[i];
171 });
172
173 return *this;
174 }
175
182 {
183 auto size = index_constant<N()>();
185 (*this)[i] -= b[i];
186 });
187
188 return *this;
189 }
190
193 template<typename X, typename Y>
194 void mv (const X& x, Y& y) const {
195 static_assert(X::size() == M(), "length of x does not match row length");
196 static_assert(Y::size() == N(), "length of y does not match row count");
197 y = 0; //reset y (for mv uses umv)
198 umv(x,y);
199 }
200
203 template<typename X, typename Y>
204 void umv (const X& x, Y& y) const {
205 static_assert(X::size() == M(), "length of x does not match row length");
206 static_assert(Y::size() == N(), "length of y does not match row count");
207 using namespace Dune::Hybrid;
208 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
209 using namespace Dune::Hybrid; // needed for icc, see issue #31
210 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
211 (*this)[i][j].umv(x[j], y[i]);
212 });
213 });
214 }
215
218 template<typename X, typename Y>
219 void mmv (const X& x, Y& y) const {
220 static_assert(X::size() == M(), "length of x does not match row length");
221 static_assert(Y::size() == N(), "length of y does not match row count");
222 using namespace Dune::Hybrid;
223 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
224 using namespace Dune::Hybrid; // needed for icc, see issue #31
225 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
226 (*this)[i][j].mmv(x[j], y[i]);
227 });
228 });
229 }
230
233 template<typename AlphaType, typename X, typename Y>
234 void usmv (const AlphaType& alpha, const X& x, Y& y) const {
235 static_assert(X::size() == M(), "length of x does not match row length");
236 static_assert(Y::size() == N(), "length of y does not match row count");
237 using namespace Dune::Hybrid;
238 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
239 using namespace Dune::Hybrid; // needed for icc, see issue #31
240 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
241 (*this)[i][j].usmv(alpha, x[j], y[i]);
242 });
243 });
244 }
245
248 template<typename X, typename Y>
249 void mtv (const X& x, Y& y) const {
250 static_assert(X::size() == N(), "length of x does not match number of rows");
251 static_assert(Y::size() == M(), "length of y does not match number of columns");
252 y = 0;
253 umtv(x,y);
254 }
255
258 template<typename X, typename Y>
259 void umtv (const X& x, Y& y) const {
260 static_assert(X::size() == N(), "length of x does not match number of rows");
261 static_assert(Y::size() == M(), "length of y does not match number of columns");
262 using namespace Dune::Hybrid;
263 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
264 using namespace Dune::Hybrid; // needed for icc, see issue #31
265 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
266 (*this)[j][i].umtv(x[j], y[i]);
267 });
268 });
269 }
270
273 template<typename X, typename Y>
274 void mmtv (const X& x, Y& y) const {
275 static_assert(X::size() == N(), "length of x does not match number of rows");
276 static_assert(Y::size() == M(), "length of y does not match number of columns");
277 using namespace Dune::Hybrid;
278 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
279 using namespace Dune::Hybrid; // needed for icc, see issue #31
280 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
281 (*this)[j][i].mmtv(x[j], y[i]);
282 });
283 });
284 }
285
288 template<typename X, typename Y>
289 void usmtv (const field_type& alpha, const X& x, Y& y) const {
290 static_assert(X::size() == N(), "length of x does not match number of rows");
291 static_assert(Y::size() == M(), "length of y does not match number of columns");
292 using namespace Dune::Hybrid;
293 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
294 using namespace Dune::Hybrid; // needed for icc, see issue #31
295 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
296 (*this)[j][i].usmtv(alpha, x[j], y[i]);
297 });
298 });
299 }
300
303 template<typename X, typename Y>
304 void umhv (const X& x, Y& y) const {
305 static_assert(X::size() == N(), "length of x does not match number of rows");
306 static_assert(Y::size() == M(), "length of y does not match number of columns");
307 using namespace Dune::Hybrid;
308 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
309 using namespace Dune::Hybrid; // needed for icc, see issue #31
310 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
311 (*this)[j][i].umhv(x[j], y[i]);
312 });
313 });
314 }
315
318 template<typename X, typename Y>
319 void mmhv (const X& x, Y& y) const {
320 static_assert(X::size() == N(), "length of x does not match number of rows");
321 static_assert(Y::size() == M(), "length of y does not match number of columns");
322 using namespace Dune::Hybrid;
323 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
324 using namespace Dune::Hybrid; // needed for icc, see issue #31
325 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
326 (*this)[j][i].mmhv(x[j], y[i]);
327 });
328 });
329 }
330
333 template<typename X, typename Y>
334 void usmhv (const field_type& alpha, const X& x, Y& y) const {
335 static_assert(X::size() == N(), "length of x does not match number of rows");
336 static_assert(Y::size() == M(), "length of y does not match number of columns");
337 using namespace Dune::Hybrid;
338 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
339 using namespace Dune::Hybrid; // needed for icc, see issue #31
340 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
341 (*this)[j][i].usmhv(alpha, x[j], y[i]);
342 });
343 });
344 }
345
346
347 //===== norms
348
350 auto frobenius_norm2 () const
351 {
352 using field_type_00 = typename std::decay_t<decltype((*this)[Indices::_0][Indices::_0])>::field_type;
353 typename FieldTraits<field_type_00>::real_type sum=0;
354
355 auto rows = index_constant<N()>();
356 Hybrid::forEach(Hybrid::integralRange(rows), [&](auto&& i) {
358 Hybrid::forEach(Hybrid::integralRange(cols), [&](auto&& j) {
359 sum += (*this)[i][j].frobenius_norm2();
360 });
361 });
362
363 return sum;
364 }
365
367 typename FieldTraits<field_type>::real_type frobenius_norm () const
368 {
369 return sqrt(frobenius_norm2());
370 }
371
373 auto infinity_norm () const
374 {
375 using field_type_00 = typename std::decay_t<decltype((*this)[Indices::_0][Indices::_0])>::field_type;
376 using std::max;
377 typename FieldTraits<field_type_00>::real_type norm=0;
378
379 auto rows = index_constant<N()>();
380 Hybrid::forEach(Hybrid::integralRange(rows), [&](auto&& i) {
381
382 typename FieldTraits<field_type_00>::real_type sum=0;
384 Hybrid::forEach(Hybrid::integralRange(cols), [&](auto&& j) {
385 sum += (*this)[i][j].infinity_norm();
386 });
387 norm = max(sum, norm);
388 });
389
390 return norm;
391 }
392
394 auto infinity_norm_real () const
395 {
396 using field_type_00 = typename std::decay_t<decltype((*this)[Indices::_0][Indices::_0])>::field_type;
397 using std::max;
398 typename FieldTraits<field_type_00>::real_type norm=0;
399
400 auto rows = index_constant<N()>();
401 Hybrid::forEach(Hybrid::integralRange(rows), [&](auto&& i) {
402
403 typename FieldTraits<field_type_00>::real_type sum=0;
405 Hybrid::forEach(Hybrid::integralRange(cols), [&](auto&& j) {
406 sum += (*this)[i][j].infinity_norm_real();
407 });
408 norm = max(sum, norm);
409 });
410
411 return norm;
412 }
413
414 };
415
421 template<typename T1, typename... Args>
422 std::ostream& operator<< (std::ostream& s, const MultiTypeBlockMatrix<T1,Args...>& m) {
425 using namespace Dune::Hybrid;
426 forEach(integralRange(N), [&](auto&& i) {
427 using namespace Dune::Hybrid; // needed for icc, see issue #31
428 forEach(integralRange(M), [&](auto&& j) {
429 s << "\t(" << i << ", " << j << "): \n" << m[i][j];
430 });
431 });
432 s << std::endl;
433 return s;
434 }
435
436 //make algmeta_itsteps known
437 template<int I, typename M>
438 struct algmeta_itsteps;
439
446 template<int I, int crow, int ccol, int remain_col> //MultiTypeBlockMatrix_Solver_Col: iterating over one row
447 class MultiTypeBlockMatrix_Solver_Col { //calculating b- A[i][j]*x[j]
448 public:
452 template <typename Trhs, typename TVector, typename TMatrix, typename K>
453 static void calc_rhs(const TMatrix& A, TVector& x, TVector& v, Trhs& b, const K& w) {
454 std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
456 }
457
458 };
459 template<int I, int crow, int ccol> //MultiTypeBlockMatrix_Solver_Col recursion end
460 class MultiTypeBlockMatrix_Solver_Col<I,crow,ccol,0> {
461 public:
462 template <typename Trhs, typename TVector, typename TMatrix, typename K>
463 static void calc_rhs(const TMatrix&, TVector&, TVector&, Trhs&, const K&) {}
464 };
465
466
467
474 template<int I, int crow, int remain_row>
476 public:
477
481 template <typename TVector, typename TMatrix, typename K>
482 static void dbgs(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
483 TVector xold(x);
484 xold=x; //store old x values
486 x *= w;
487 x.axpy(1-w,xold); //improve x
488 }
489 template <typename TVector, typename TMatrix, typename K>
490 static void dbgs(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
491 auto rhs = std::get<crow> (b);
492
493 MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
494 //solve on blocklevel I-1
495 using M =
496 typename std::remove_cv<
497 typename std::remove_reference<
498 decltype(std::get<crow>( std::get<crow>(A)))
499 >::type
500 >::type;
501 algmeta_itsteps<I-1,M>::dbgs(std::get<crow>( std::get<crow>(A)), std::get<crow>(x),rhs,w);
503 }
504
505
506
510 template <typename TVector, typename TMatrix, typename K>
511 static void bsorf(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
512 TVector v;
513 v=x; //use latest x values in right side calculation
515
516 }
517 template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
518 static void bsorf(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
519 auto rhs = std::get<crow> (b);
520
521 MultiTypeBlockMatrix_Solver_Col<I,crow,0,TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
522 //solve on blocklevel I-1
523 using M =
524 typename std::remove_cv<
525 typename std::remove_reference<
526 decltype(std::get<crow>( std::get<crow>(A)))
527 >::type
528 >::type;
529 algmeta_itsteps<I-1,M>::bsorf(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
530 std::get<crow>(x).axpy(w,std::get<crow>(v));
532 }
533
537 template <typename TVector, typename TMatrix, typename K>
538 static void bsorb(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
539 TVector v;
540 v=x; //use latest x values in right side calculation
542
543 }
544 template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
545 static void bsorb(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
546 auto rhs = std::get<crow> (b);
547
548 MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
549 //solve on blocklevel I-1
550 using M =
551 typename std::remove_cv<
552 typename std::remove_reference<
553 decltype(std::get<crow>( std::get<crow>(A)))
554 >::type
555 >::type;
556 algmeta_itsteps<I-1,M>::bsorb(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
557 std::get<crow>(x).axpy(w,std::get<crow>(v));
559 }
560
561
565 template <typename TVector, typename TMatrix, typename K>
566 static void dbjac(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
567 TVector v(x);
568 v=0; //calc new x in v
570 x.axpy(w,v); //improve x
571 }
572 template <typename TVector, typename TMatrix, typename K>
573 static void dbjac(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
574 auto rhs = std::get<crow> (b);
575
576 MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
577 //solve on blocklevel I-1
578 using M =
579 typename std::remove_cv<
580 typename std::remove_reference<
581 decltype(std::get<crow>( std::get<crow>(A)))
582 >::type
583 >::type;
584 algmeta_itsteps<I-1,M>::dbjac(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
586 }
587
588
589
590
591 };
592 template<int I, int crow> //recursion end for remain_row = 0
593 class MultiTypeBlockMatrix_Solver<I,crow,0> {
594 public:
595 template <typename TVector, typename TMatrix, typename K>
596 static void dbgs(const TMatrix&, TVector&, TVector&,
597 const TVector&, const K&) {}
598
599 template <typename TVector, typename TMatrix, typename K>
600 static void bsorf(const TMatrix&, TVector&, TVector&,
601 const TVector&, const K&) {}
602
603 template <typename TVector, typename TMatrix, typename K>
604 static void bsorb(const TMatrix&, TVector&, TVector&,
605 const TVector&, const K&) {}
606
607 template <typename TVector, typename TMatrix, typename K>
608 static void dbjac(const TMatrix&, TVector&, TVector&,
609 const TVector&, const K&) {}
610 };
611
612} // end namespace
613
614#endif
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:447
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:475
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:44
constexpr index_constant< 0 > _0
Compile time index with value 0.
Definition: indices.hh:51
std::integral_constant< std::size_t, i > index_constant
An index constant with value i.
Definition: indices.hh:28
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:266
constexpr auto integralRange(const Begin &begin, const End &end)
Create an integral range.
Definition: hybridutilities.hh:182
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:644
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:656
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:620
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:632
MultiTypeBlockMatrix< FirstRow, Args... > type
Definition: multitypeblockmatrix.hh:54
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:566
MultiTypeBlockMatrix & operator+=(const MultiTypeBlockMatrix &b)
Add the entries of another matrix to this one.
Definition: multitypeblockmatrix.hh:166
void mv(const X &x, Y &y) const
y = A x
Definition: multitypeblockmatrix.hh:194
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: multitypeblockmatrix.hh:274
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: multitypeblockmatrix.hh:319
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:482
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:538
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: multitypeblockmatrix.hh:334
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:234
MultiTypeBlockMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: multitypeblockmatrix.hh:150
void umhv(const X &x, Y &y) const
y += A^H x
Definition: multitypeblockmatrix.hh:304
FieldTraits< field_type >::real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: multitypeblockmatrix.hh:367
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:62
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: multitypeblockmatrix.hh:289
void operator=(const T &newval)
Definition: multitypeblockmatrix.hh:125
void umv(const X &x, Y &y) const
y += A x
Definition: multitypeblockmatrix.hh:204
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition: multitypeblockmatrix.hh:453
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:511
static constexpr size_type size()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:73
auto infinity_norm_real() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:394
auto frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: multitypeblockmatrix.hh:350
std::size_t size_type
Type used for sizes.
Definition: multitypeblockmatrix.hh:57
auto infinity_norm() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:373
MultiTypeBlockMatrix & operator-=(const MultiTypeBlockMatrix &b)
Subtract the entries of another matrix from this one.
Definition: multitypeblockmatrix.hh:181
auto operator[](const std::integral_constant< size_type, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:102
void mtv(const X &x, Y &y) const
y = A^T x
Definition: multitypeblockmatrix.hh:249
static constexpr size_type M()
Return the number of matrix columns.
Definition: multitypeblockmatrix.hh:79
void mmv(const X &x, Y &y) const
y -= A x
Definition: multitypeblockmatrix.hh:219
void umtv(const X &x, Y &y) const
y += A^T x
Definition: multitypeblockmatrix.hh:259
MultiTypeBlockMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: multitypeblockmatrix.hh:139
auto max(ADLTag< 0 >, const V &v1, const V &v2)
implements binary Simd::max()
Definition: defaults.hh:79
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignedallocator.hh:11
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)