DUNE PDELab (git)

multitypeblockmatrix.hh
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_MULTITYPEBLOCKMATRIX_HH
6#define DUNE_ISTL_MULTITYPEBLOCKMATRIX_HH
7
8#include <cmath>
9#include <iostream>
10#include <tuple>
11
12#include <dune/common/hybridutilities.hh>
13
14#include "istlexception.hh"
15
16// forward declaration
17namespace Dune
18{
19 template<typename FirstRow, typename... Args>
20 class MultiTypeBlockMatrix;
21
22 template<int I, int crow, int remain_row>
23 class MultiTypeBlockMatrix_Solver;
24}
25
26#include "gsetc.hh"
27
28namespace Dune {
29
43 template<typename FirstRow, typename... Args>
45 : public std::tuple<FirstRow, Args...>
46 {
47 using ParentType = std::tuple<FirstRow, Args...>;
48 public:
49
51 using ParentType::ParentType;
52
56 typedef MultiTypeBlockMatrix<FirstRow, Args...> type;
57
59 using size_type = std::size_t;
60
66 using field_type = Std::detected_t<std::common_type_t,
67 typename FieldTraits<FirstRow>::field_type, typename FieldTraits<Args>::field_type...>;
68
74 using real_type = Std::detected_t<std::common_type_t,
75 typename FieldTraits<FirstRow>::real_type, typename FieldTraits<Args>::real_type...>;
76
77 // make sure that we have an std::common_type: using an assertion produces a more readable error message
78 // than a compiler template instantiation error
79 static_assert ( sizeof...(Args) == 0 or
80 not (std::is_same_v<field_type, Std::nonesuch> or std::is_same_v<real_type, Std::nonesuch>),
81 "No std::common_type implemented for the given field_type/real_type of the Args. Please provide an implementation and check that a FieldTraits class is present for your type.");
82
84 static constexpr size_type N()
85 {
86 return 1+sizeof...(Args);
87 }
88
90 static constexpr size_type M()
91 {
92 return FirstRow::size();
93 }
94
111 template< size_type index >
112 auto
113 operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable)
114 -> decltype(std::get<index>(*this))
115 {
116 return std::get<index>(*this);
117 }
118
124 template< size_type index >
125 auto
126 operator[] ([[maybe_unused]] const std::integral_constant< size_type, index > indexVariable) const
127 -> decltype(std::get<index>(*this))
128 {
129 return std::get<index>(*this);
130 }
131
135 template<typename T>
136 void operator= (const T& newval) {
137 using namespace Dune::Hybrid;
138 auto size = index_constant<1+sizeof...(Args)>();
139 // Since Dune::Hybrid::size(MultiTypeBlockMatrix) is not implemented,
140 // we cannot use a plain forEach(*this, ...). This could be achieved,
141 // e.g., by implementing a static size() method.
142 forEach(integralRange(size), [&](auto&& i) {
143 (*this)[i] = newval;
144 });
145 }
146
147 //===== vector space arithmetic
148
151 {
152 auto size = index_constant<N()>();
154 (*this)[i] *= k;
155 });
156
157 return *this;
158 }
159
162 {
163 auto size = index_constant<N()>();
165 (*this)[i] /= k;
166 });
167
168 return *this;
169 }
170
171
178 {
179 auto size = index_constant<N()>();
181 (*this)[i] += b[i];
182 });
183
184 return *this;
185 }
186
193 {
194 auto size = index_constant<N()>();
196 (*this)[i] -= b[i];
197 });
198
199 return *this;
200 }
201
204 template<typename X, typename Y>
205 void mv (const X& x, Y& y) const {
206 static_assert(X::size() == M(), "length of x does not match row length");
207 static_assert(Y::size() == N(), "length of y does not match row count");
208 y = 0; //reset y (for mv uses umv)
209 umv(x,y);
210 }
211
214 template<typename X, typename Y>
215 void umv (const X& x, Y& y) const {
216 static_assert(X::size() == M(), "length of x does not match row length");
217 static_assert(Y::size() == N(), "length of y does not match row count");
218 using namespace Dune::Hybrid;
219 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
220 using namespace Dune::Hybrid; // needed for icc, see issue #31
221 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
222 (*this)[i][j].umv(x[j], y[i]);
223 });
224 });
225 }
226
229 template<typename X, typename Y>
230 void mmv (const X& x, Y& y) const {
231 static_assert(X::size() == M(), "length of x does not match row length");
232 static_assert(Y::size() == N(), "length of y does not match row count");
233 using namespace Dune::Hybrid;
234 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
235 using namespace Dune::Hybrid; // needed for icc, see issue #31
236 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
237 (*this)[i][j].mmv(x[j], y[i]);
238 });
239 });
240 }
241
244 template<typename AlphaType, typename X, typename Y>
245 void usmv (const AlphaType& alpha, const X& x, Y& y) const {
246 static_assert(X::size() == M(), "length of x does not match row length");
247 static_assert(Y::size() == N(), "length of y does not match row count");
248 using namespace Dune::Hybrid;
249 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
250 using namespace Dune::Hybrid; // needed for icc, see issue #31
251 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
252 (*this)[i][j].usmv(alpha, x[j], y[i]);
253 });
254 });
255 }
256
259 template<typename X, typename Y>
260 void mtv (const X& x, Y& y) const {
261 static_assert(X::size() == N(), "length of x does not match number of rows");
262 static_assert(Y::size() == M(), "length of y does not match number of columns");
263 y = 0;
264 umtv(x,y);
265 }
266
269 template<typename X, typename Y>
270 void umtv (const X& x, Y& y) const {
271 static_assert(X::size() == N(), "length of x does not match number of rows");
272 static_assert(Y::size() == M(), "length of y does not match number of columns");
273 using namespace Dune::Hybrid;
274 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
275 using namespace Dune::Hybrid; // needed for icc, see issue #31
276 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
277 (*this)[j][i].umtv(x[j], y[i]);
278 });
279 });
280 }
281
284 template<typename X, typename Y>
285 void mmtv (const X& x, Y& y) const {
286 static_assert(X::size() == N(), "length of x does not match number of rows");
287 static_assert(Y::size() == M(), "length of y does not match number of columns");
288 using namespace Dune::Hybrid;
289 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
290 using namespace Dune::Hybrid; // needed for icc, see issue #31
291 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
292 (*this)[j][i].mmtv(x[j], y[i]);
293 });
294 });
295 }
296
299 template<typename X, typename Y>
300 void usmtv (const field_type& alpha, const X& x, Y& y) const {
301 static_assert(X::size() == N(), "length of x does not match number of rows");
302 static_assert(Y::size() == M(), "length of y does not match number of columns");
303 using namespace Dune::Hybrid;
304 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
305 using namespace Dune::Hybrid; // needed for icc, see issue #31
306 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
307 (*this)[j][i].usmtv(alpha, x[j], y[i]);
308 });
309 });
310 }
311
314 template<typename X, typename Y>
315 void umhv (const X& x, Y& y) const {
316 static_assert(X::size() == N(), "length of x does not match number of rows");
317 static_assert(Y::size() == M(), "length of y does not match number of columns");
318 using namespace Dune::Hybrid;
319 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
320 using namespace Dune::Hybrid; // needed for icc, see issue #31
321 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
322 (*this)[j][i].umhv(x[j], y[i]);
323 });
324 });
325 }
326
329 template<typename X, typename Y>
330 void mmhv (const X& x, Y& y) const {
331 static_assert(X::size() == N(), "length of x does not match number of rows");
332 static_assert(Y::size() == M(), "length of y does not match number of columns");
333 using namespace Dune::Hybrid;
334 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
335 using namespace Dune::Hybrid; // needed for icc, see issue #31
336 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
337 (*this)[j][i].mmhv(x[j], y[i]);
338 });
339 });
340 }
341
344 template<typename X, typename Y>
345 void usmhv (const field_type& alpha, const X& x, Y& y) const {
346 static_assert(X::size() == N(), "length of x does not match number of rows");
347 static_assert(Y::size() == M(), "length of y does not match number of columns");
348 using namespace Dune::Hybrid;
349 forEach(integralRange(Hybrid::size(y)), [&](auto&& i) {
350 using namespace Dune::Hybrid; // needed for icc, see issue #31
351 forEach(integralRange(Hybrid::size(x)), [&](auto&& j) {
352 (*this)[j][i].usmhv(alpha, x[j], y[i]);
353 });
354 });
355 }
356
357
358 //===== norms
359
362 {
363 real_type sum=0;
364
365 auto rows = index_constant<N()>();
366 Hybrid::forEach(Hybrid::integralRange(rows), [&](auto&& i) {
368 Hybrid::forEach(Hybrid::integralRange(cols), [&](auto&& j) {
369 sum += (*this)[i][j].frobenius_norm2();
370 });
371 });
372
373 return sum;
374 }
375
378 {
379 return sqrt(frobenius_norm2());
380 }
381
384 {
385 using std::max;
386 real_type norm=0;
387
388 auto rows = index_constant<N()>();
389 Hybrid::forEach(Hybrid::integralRange(rows), [&](auto&& i) {
390 real_type sum=0;
392 Hybrid::forEach(Hybrid::integralRange(cols), [&](auto&& j) {
393 sum += (*this)[i][j].infinity_norm();
394 });
395 norm = max(sum, norm);
396 });
397
398 return norm;
399 }
400
403 {
404 using std::max;
405 real_type norm=0;
406
407 auto rows = index_constant<N()>();
408 Hybrid::forEach(Hybrid::integralRange(rows), [&](auto&& i) {
409 real_type sum=0;
411 Hybrid::forEach(Hybrid::integralRange(cols), [&](auto&& j) {
412 sum += (*this)[i][j].infinity_norm_real();
413 });
414 norm = max(sum, norm);
415 });
416
417 return norm;
418 }
419
420 };
421
422
426 template <typename... Rows>
427 struct FieldTraits< MultiTypeBlockMatrix<Rows...> >
428 {
429 using field_type = typename MultiTypeBlockMatrix<Rows...>::field_type;
430 using real_type = typename MultiTypeBlockMatrix<Rows...>::real_type;
431 };
432
433
439 template<typename T1, typename... Args>
440 std::ostream& operator<< (std::ostream& s, const MultiTypeBlockMatrix<T1,Args...>& m) {
443 using namespace Dune::Hybrid;
444 forEach(integralRange(N), [&](auto&& i) {
445 using namespace Dune::Hybrid; // needed for icc, see issue #31
446 forEach(integralRange(M), [&](auto&& j) {
447 s << "\t(" << i << ", " << j << "): \n" << m[i][j];
448 });
449 });
450 s << std::endl;
451 return s;
452 }
453
454 //make algmeta_itsteps known
455 template<int I, typename M>
456 struct algmeta_itsteps;
457
464 template<int I, int crow, int ccol, int remain_col> //MultiTypeBlockMatrix_Solver_Col: iterating over one row
465 class MultiTypeBlockMatrix_Solver_Col { //calculating b- A[i][j]*x[j]
466 public:
470 template <typename Trhs, typename TVector, typename TMatrix, typename K>
471 static void calc_rhs(const TMatrix& A, TVector& x, TVector& v, Trhs& b, const K& w) {
472 std::get<ccol>( std::get<crow>(A) ).mmv( std::get<ccol>(x), b );
474 }
475
476 };
477 template<int I, int crow, int ccol> //MultiTypeBlockMatrix_Solver_Col recursion end
478 class MultiTypeBlockMatrix_Solver_Col<I,crow,ccol,0> {
479 public:
480 template <typename Trhs, typename TVector, typename TMatrix, typename K>
481 static void calc_rhs(const TMatrix&, TVector&, TVector&, Trhs&, const K&) {}
482 };
483
484
485
492 template<int I, int crow, int remain_row>
494 public:
495
499 template <typename TVector, typename TMatrix, typename K>
500 static void dbgs(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
501 TVector xold(x);
502 xold=x; //store old x values
504 x *= w;
505 x.axpy(1-w,xold); //improve x
506 }
507 template <typename TVector, typename TMatrix, typename K>
508 static void dbgs(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
509 auto rhs = std::get<crow> (b);
510
511 MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
512 //solve on blocklevel I-1
513 using M =
514 typename std::remove_cv<
515 typename std::remove_reference<
516 decltype(std::get<crow>( std::get<crow>(A)))
517 >::type
518 >::type;
519 algmeta_itsteps<I-1,M>::dbgs(std::get<crow>( std::get<crow>(A)), std::get<crow>(x),rhs,w);
521 }
522
523
524
528 template <typename TVector, typename TMatrix, typename K>
529 static void bsorf(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
530 TVector v;
531 v=x; //use latest x values in right side calculation
533
534 }
535 template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
536 static void bsorf(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
537 auto rhs = std::get<crow> (b);
538
539 MultiTypeBlockMatrix_Solver_Col<I,crow,0,TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
540 //solve on blocklevel I-1
541 using M =
542 typename std::remove_cv<
543 typename std::remove_reference<
544 decltype(std::get<crow>( std::get<crow>(A)))
545 >::type
546 >::type;
547 algmeta_itsteps<I-1,M>::bsorf(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
548 std::get<crow>(x).axpy(w,std::get<crow>(v));
550 }
551
555 template <typename TVector, typename TMatrix, typename K>
556 static void bsorb(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
557 TVector v;
558 v=x; //use latest x values in right side calculation
560
561 }
562 template <typename TVector, typename TMatrix, typename K> //recursion over all matrix rows (A)
563 static void bsorb(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
564 auto rhs = std::get<crow> (b);
565
566 MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
567 //solve on blocklevel I-1
568 using M =
569 typename std::remove_cv<
570 typename std::remove_reference<
571 decltype(std::get<crow>( std::get<crow>(A)))
572 >::type
573 >::type;
574 algmeta_itsteps<I-1,M>::bsorb(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
575 std::get<crow>(x).axpy(w,std::get<crow>(v));
577 }
578
579
583 template <typename TVector, typename TMatrix, typename K>
584 static void dbjac(const TMatrix& A, TVector& x, const TVector& b, const K& w) {
585 TVector v(x);
586 v=0; //calc new x in v
588 x.axpy(w,v); //improve x
589 }
590 template <typename TVector, typename TMatrix, typename K>
591 static void dbjac(const TMatrix& A, TVector& x, TVector& v, const TVector& b, const K& w) {
592 auto rhs = std::get<crow> (b);
593
594 MultiTypeBlockMatrix_Solver_Col<I,crow,0, TMatrix::M()>::calc_rhs(A,x,v,rhs,w); // calculate right side of equation
595 //solve on blocklevel I-1
596 using M =
597 typename std::remove_cv<
598 typename std::remove_reference<
599 decltype(std::get<crow>( std::get<crow>(A)))
600 >::type
601 >::type;
602 algmeta_itsteps<I-1,M>::dbjac(std::get<crow>( std::get<crow>(A)), std::get<crow>(v),rhs,w);
604 }
605
606
607
608
609 };
610 template<int I, int crow> //recursion end for remain_row = 0
611 class MultiTypeBlockMatrix_Solver<I,crow,0> {
612 public:
613 template <typename TVector, typename TMatrix, typename K>
614 static void dbgs(const TMatrix&, TVector&, TVector&,
615 const TVector&, const K&) {}
616
617 template <typename TVector, typename TMatrix, typename K>
618 static void bsorf(const TMatrix&, TVector&, TVector&,
619 const TVector&, const K&) {}
620
621 template <typename TVector, typename TMatrix, typename K>
622 static void bsorb(const TMatrix&, TVector&, TVector&,
623 const TVector&, const K&) {}
624
625 template <typename TVector, typename TMatrix, typename K>
626 static void dbjac(const TMatrix&, TVector&, TVector&,
627 const TVector&, const K&) {}
628 };
629
630} // end namespace Dune
631
632namespace std
633{
638 template <size_t i, typename... Args>
639 struct tuple_element<i,Dune::MultiTypeBlockMatrix<Args...> >
640 {
641 using type = typename std::tuple_element<i, std::tuple<Args...> >::type;
642 };
643
648 template <typename... Args>
649 struct tuple_size<Dune::MultiTypeBlockMatrix<Args...> >
650 : std::integral_constant<std::size_t, sizeof...(Args)>
651 {};
652}
653#endif
part of solvers for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:465
solver for MultiTypeBlockVector & MultiTypeBlockMatrix types
Definition: multitypeblockmatrix.hh:493
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:46
std::integral_constant< std::size_t, i > index_constant
An index constant with value i.
Definition: indices.hh:29
typename detected_or< nonesuch, Op, Args... >::type detected_t
Returns Op<Args...> if that is valid; otherwise returns nonesuch.
Definition: type_traits.hh:174
MultiTypeBlockMatrix< FirstRow, Args... > type
Definition: multitypeblockmatrix.hh:56
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:584
real_type infinity_norm_real() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:402
MultiTypeBlockMatrix & operator+=(const MultiTypeBlockMatrix &b)
Add the entries of another matrix to this one.
Definition: multitypeblockmatrix.hh:177
void mv(const X &x, Y &y) const
y = A x
Definition: multitypeblockmatrix.hh:205
void mmtv(const X &x, Y &y) const
y -= A^T x
Definition: multitypeblockmatrix.hh:285
void mmhv(const X &x, Y &y) const
y -= A^H x
Definition: multitypeblockmatrix.hh:330
Std::detected_t< std::common_type_t, typename FieldTraits< FirstRow >::real_type, typename FieldTraits< Args >::real_type... > real_type
The type used for real values.
Definition: multitypeblockmatrix.hh:75
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:500
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:556
void usmhv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^H x
Definition: multitypeblockmatrix.hh:345
void usmv(const AlphaType &alpha, const X &x, Y &y) const
y += alpha A x
Definition: multitypeblockmatrix.hh:245
MultiTypeBlockMatrix & operator/=(const field_type &k)
vector space division by scalar
Definition: multitypeblockmatrix.hh:161
real_type infinity_norm() const
Bastardized version of the infinity-norm / row-sum norm.
Definition: multitypeblockmatrix.hh:383
void umhv(const X &x, Y &y) const
y += A^H x
Definition: multitypeblockmatrix.hh:315
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:84
void usmtv(const field_type &alpha, const X &x, Y &y) const
y += alpha A^T x
Definition: multitypeblockmatrix.hh:300
void operator=(const T &newval)
Definition: multitypeblockmatrix.hh:136
void umv(const X &x, Y &y) const
y += A x
Definition: multitypeblockmatrix.hh:215
static void calc_rhs(const TMatrix &A, TVector &x, TVector &v, Trhs &b, const K &w)
Definition: multitypeblockmatrix.hh:471
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:529
std::size_t size_type
Type used for sizes.
Definition: multitypeblockmatrix.hh:59
real_type frobenius_norm2() const
square of frobenius norm, need for block recursion
Definition: multitypeblockmatrix.hh:361
real_type frobenius_norm() const
frobenius norm: sqrt(sum over squared values of entries)
Definition: multitypeblockmatrix.hh:377
MultiTypeBlockMatrix & operator-=(const MultiTypeBlockMatrix &b)
Subtract the entries of another matrix from this one.
Definition: multitypeblockmatrix.hh:192
auto operator[](const std::integral_constant< size_type, index > indexVariable) -> decltype(std::get< index >(*this))
Random-access operator.
Definition: multitypeblockmatrix.hh:113
void mtv(const X &x, Y &y) const
y = A^T x
Definition: multitypeblockmatrix.hh:260
static constexpr size_type M()
Return the number of matrix columns.
Definition: multitypeblockmatrix.hh:90
void mmv(const X &x, Y &y) const
y -= A x
Definition: multitypeblockmatrix.hh:230
Std::detected_t< std::common_type_t, typename FieldTraits< FirstRow >::field_type, typename FieldTraits< Args >::field_type... > field_type
The type used for scalars.
Definition: multitypeblockmatrix.hh:67
void umtv(const X &x, Y &y) const
y += A^T x
Definition: multitypeblockmatrix.hh:270
MultiTypeBlockMatrix & operator*=(const field_type &k)
vector space multiplication with scalar
Definition: multitypeblockmatrix.hh:150
constexpr void forEach(Range &&range, F &&f)
Range based for loop.
Definition: hybridutilities.hh:256
constexpr auto max
Function object that returns the greater of the given values.
Definition: hybridutilities.hh:484
constexpr auto integralRange(const Begin &begin, const End &end)
Create an integral range.
Definition: hybridutilities.hh:172
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:646
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:658
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:622
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:634
Simple iterative methods like Jacobi, Gauss-Seidel, SOR, SSOR, etc. in a generic way.
Dune namespace.
Definition: alignedallocator.hh:13
constexpr std::integral_constant< std::size_t, sizeof...(II)> size(std::integer_sequence< T, II... >)
Return the size of the sequence.
Definition: integersequence.hh:75
STL namespace.
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 24, 23:30, 2024)