3#ifndef DUNE_ISTL_GSETC_HH
4#define DUNE_ISTL_GSETC_HH
12#include <dune/common/hybridutilities.hh>
14#include "multitypeblockvector.hh"
15#include "multitypeblockmatrix.hh"
17#include "istlexception.hh"
44 enum {recursion_level = l};
66 template<
int I, WithDiagType diag, WithRelaxType relax>
67 struct algmeta_btsolve {
68 template<
class M,
class X,
class Y,
class K>
69 static void bltsolve (
const M& A, X& v,
const Y& d,
const K& w)
72 typedef typename M::ConstRowIterator rowiterator;
73 typedef typename M::ConstColIterator coliterator;
74 typedef typename Y::block_type bblock;
77 rowiterator endi=A.end();
78 for (rowiterator i=A.begin(); i!=endi; ++i)
80 bblock rhs(d[i.index()]);
82 for (j=(*i).begin(); j.index()<i.index(); ++j)
83 (*j).mmv(v[j.index()],rhs);
87 template<
class M,
class X,
class Y,
class K>
88 static void butsolve (
const M& A, X& v,
const Y& d,
const K& w)
91 typedef typename M::ConstRowIterator rowiterator;
92 typedef typename M::ConstColIterator coliterator;
93 typedef typename Y::block_type bblock;
96 rowiterator rendi=A.beforeBegin();
97 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
99 bblock rhs(d[i.index()]);
101 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
102 (*j).mmv(v[j.index()],rhs);
110 struct algmeta_btsolve<0,withdiag,withrelax> {
111 template<
class M,
class X,
class Y,
class K>
112 static void bltsolve (
const M& A, X& v,
const Y& d,
const K& w)
117 template<
class M,
class X,
class Y,
class K>
118 static void butsolve (
const M& A, X& v,
const Y& d,
const K& w)
125 struct algmeta_btsolve<0,withdiag,norelax> {
126 template<
class M,
class X,
class Y,
class K>
127 static void bltsolve (
const M& A, X& v,
const Y& d,
const K& )
131 template<
class M,
class X,
class Y,
class K>
132 static void butsolve (
const M& A, X& v,
const Y& d,
const K& )
138 struct algmeta_btsolve<0,nodiag,withrelax> {
139 template<
class M,
class X,
class Y,
class K>
140 static void bltsolve (
const M& , X& v,
const Y& d,
const K& w)
145 template<
class M,
class X,
class Y,
class K>
146 static void butsolve (
const M& , X& v,
const Y& d,
const K& w)
153 struct algmeta_btsolve<0,nodiag,norelax> {
154 template<
class M,
class X,
class Y,
class K>
155 static void bltsolve (
const M& , X& v,
const Y& d,
const K& )
159 template<
class M,
class X,
class Y,
class K>
160 static void butsolve (
const M& , X& v,
const Y& d,
const K& )
172 template<
class M,
class X,
class Y>
175 typename X::field_type w=1;
179 template<
class M,
class X,
class Y,
class K>
180 void bltsolve (
const M& A, X& v,
const Y& d,
const K& w)
185 template<
class M,
class X,
class Y>
188 typename X::field_type w=1;
192 template<
class M,
class X,
class Y,
class K>
193 void ubltsolve (
const M& A, X& v,
const Y& d,
const K& w)
199 template<
class M,
class X,
class Y>
202 typename X::field_type w=1;
206 template<
class M,
class X,
class Y,
class K>
207 void butsolve (
const M& A, X& v,
const Y& d,
const K& w)
212 template<
class M,
class X,
class Y>
215 typename X::field_type w=1;
219 template<
class M,
class X,
class Y,
class K>
220 void ubutsolve (
const M& A, X& v,
const Y& d,
const K& w)
228 template<
class M,
class X,
class Y,
int l>
231 typename X::field_type w=1;
235 template<
class M,
class X,
class Y,
class K,
int l>
241 template<
class M,
class X,
class Y,
int l>
244 typename X::field_type w=1;
248 template<
class M,
class X,
class Y,
class K,
int l>
255 template<
class M,
class X,
class Y,
int l>
258 typename X::field_type w=1;
262 template<
class M,
class X,
class Y,
class K,
int l>
268 template<
class M,
class X,
class Y,
int l>
271 typename X::field_type w=1;
275 template<
class M,
class X,
class Y,
class K,
int l>
291 template<
int I, WithRelaxType relax>
292 struct algmeta_bdsolve {
293 template<
class M,
class X,
class Y,
class K>
294 static void bdsolve (
const M& A, X& v,
const Y& d,
const K& w)
297 typedef typename M::ConstRowIterator rowiterator;
298 typedef typename M::ConstColIterator coliterator;
301 rowiterator rendi=A.beforeBegin();
302 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
304 coliterator ii=(*i).find(i.index());
312 struct algmeta_bdsolve<0,withrelax> {
313 template<
class M,
class X,
class Y,
class K>
314 static void bdsolve (
const M& A, X& v,
const Y& d,
const K& w)
321 struct algmeta_bdsolve<0,norelax> {
322 template<
class M,
class X,
class Y,
class K>
323 static void bdsolve (
const M& A, X& v,
const Y& d,
const K& )
334 template<
class M,
class X,
class Y>
337 typename X::field_type w=1;
341 template<
class M,
class X,
class Y,
class K>
342 void bdsolve (
const M& A, X& v,
const Y& d,
const K& w)
350 template<
class M,
class X,
class Y,
int l>
353 typename X::field_type w=1;
357 template<
class M,
class X,
class Y,
class K,
int l>
372 template<
int I,
typename M>
373 struct algmeta_itsteps {
375 template<
class X,
class Y,
class K>
376 static void dbgs (
const M& A, X& x,
const Y& b,
const K& w)
378 typedef typename M::ConstRowIterator rowiterator;
379 typedef typename M::ConstColIterator coliterator;
380 typedef typename Y::block_type bblock;
385 rowiterator endi=A.end();
386 for (rowiterator i=A.begin(); i!=endi; ++i)
389 coliterator endj=(*i).end();
390 coliterator j=(*i).begin();
391 if constexpr (IsNumber<typename M::block_type>())
393 for (; j.index()<i.index(); ++j)
394 rhs -= (*j) * x[j.index()];
395 coliterator diag=j++;
396 for (; j != endj; ++j)
397 rhs -= (*j) * x[j.index()];
398 x[i.index()] = rhs / (*diag);
402 for (; j.index()<i.index(); ++j)
403 (*j).mmv(x[j.index()],rhs);
404 coliterator diag=j++;
405 for (; j != endj; ++j)
406 (*j).mmv(x[j.index()],rhs);
415 template<
class X,
class Y,
class K>
416 static void bsorf (
const M& A, X& x,
const Y& b,
const K& w)
418 typedef typename M::ConstRowIterator rowiterator;
419 typedef typename M::ConstColIterator coliterator;
420 typedef typename Y::block_type bblock;
421 typedef typename X::block_type xblock;
426 if(A.begin()!=A.end())
429 rowiterator endi=A.end();
430 for (rowiterator i=A.begin(); i!=endi; ++i)
433 coliterator endj=(*i).end();
434 coliterator j=(*i).begin();
435 if constexpr (IsNumber<typename M::block_type>())
437 for (; j.index()<i.index(); ++j)
438 rhs -= (*j) * x[j.index()];
441 rhs -= (*j) * x[j.index()];
447 for (; j.index()<i.index(); ++j)
448 (*j).mmv(x[j.index()],rhs);
451 (*j).mmv(x[j.index()],rhs);
453 x[i.index()].axpy(w,v);
458 template<
class X,
class Y,
class K>
459 static void bsorb (
const M& A, X& x,
const Y& b,
const K& w)
461 typedef typename M::ConstRowIterator rowiterator;
462 typedef typename M::ConstColIterator coliterator;
463 typedef typename Y::block_type bblock;
464 typedef typename X::block_type xblock;
469 if(A.begin()!=A.end())
472 rowiterator endi=A.beforeBegin();
473 for (rowiterator i=A.beforeEnd(); i!=endi; --i)
476 coliterator endj=(*i).end();
477 coliterator j=(*i).begin();
478 if constexpr (IsNumber<typename M::block_type>())
480 for (; j.index()<i.index(); ++j)
481 rhs -= (*j) * x[j.index()];
484 rhs -= (*j) * x[j.index()];
490 for (; j.index()<i.index(); ++j)
491 j->mmv(x[j.index()],rhs);
494 j->mmv(x[j.index()],rhs);
496 x[i.index()].axpy(w,v);
501 template<
class X,
class Y,
class K>
502 static void dbjac (
const M& A, X& x,
const Y& b,
const K& w)
504 typedef typename M::ConstRowIterator rowiterator;
505 typedef typename M::ConstColIterator coliterator;
506 typedef typename Y::block_type bblock;
511 rowiterator endi=A.end();
512 for (rowiterator i=A.begin(); i!=endi; ++i)
515 coliterator endj=(*i).end();
516 coliterator j=(*i).begin();
517 if constexpr (IsNumber<typename M::block_type>())
519 for (; j.index()<i.index(); ++j)
520 rhs -= (*j) * x[j.index()];
523 rhs -= (*j) * x[j.index()];
524 v[i.index()] = rhs / (*diag);
528 for (; j.index()<i.index(); ++j)
529 j->mmv(x[j.index()],rhs);
532 j->mmv(x[j.index()],rhs);
541 struct algmeta_itsteps<0,M> {
542 template<
class X,
class Y,
class K>
543 static void dbgs (
const M& A, X& x,
const Y& b,
const K& )
547 template<
class X,
class Y,
class K>
548 static void bsorf (
const M& A, X& x,
const Y& b,
const K& )
552 template<
class X,
class Y,
class K>
553 static void bsorb (
const M& A, X& x,
const Y& b,
const K& )
557 template<
class X,
class Y,
class K>
558 static void dbjac (
const M& A, X& x,
const Y& b,
const K& )
564 template<
int I,
typename T1,
typename... MultiTypeMatrixArgs>
565 struct algmeta_itsteps<I,MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>> {
567 typename... MultiTypeVectorArgs,
569 static void dbgs (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
570 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
571 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
579 typename... MultiTypeVectorArgs,
581 static void bsorf (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
582 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
583 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
591 typename... MultiTypeVectorArgs,
593 static void bsorb (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
594 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
595 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
603 typename... MultiTypeVectorArgs,
606 static void dbjac (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
607 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
608 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
619 template<
class M,
class X,
class Y,
class K>
620 void dbgs (
const M& A, X& x,
const Y& b,
const K& w)
625 template<
class M,
class X,
class Y,
class K,
int l>
626 void dbgs (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
631 template<
class M,
class X,
class Y,
class K>
632 void bsorf (
const M& A, X& x,
const Y& b,
const K& w)
637 template<
class M,
class X,
class Y,
class K,
int l>
638 void bsorf (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
643 template<
class M,
class X,
class Y,
class K>
644 void bsorb (
const M& A, X& x,
const Y& b,
const K& w)
649 template<
class M,
class X,
class Y,
class K,
int l>
650 void bsorb (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
652 algmeta_itsteps<l,typename std::remove_cv<M>::type>
::bsorb(A,x,b,w);
655 template<
class M,
class X,
class Y,
class K>
656 void dbjac (
const M& A, X& x,
const Y& b,
const K& w)
661 template<
class M,
class X,
class Y,
class K,
int l>
662 void dbjac (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
void bltsolve(const M &A, X &v, const Y &d)
block lower triangular solve
Definition: gsetc.hh:173
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:644
void ubltsolve(const M &A, X &v, const Y &d)
unit block lower triangular solve
Definition: gsetc.hh:186
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, BL< l >)
SOR step.
Definition: gsetc.hh:638
void dbjac(const M &A, X &x, const Y &b, const K &w, BL< l >)
Jacobi step.
Definition: gsetc.hh:662
void dbgs(const M &A, X &x, const Y &b, const K &w, BL< l >)
GS step.
Definition: gsetc.hh:626
void bltsolve(const M &A, X &v, const Y &d, const K &w, BL< l >)
relaxed block lower triangular solve
Definition: gsetc.hh:236
void butsolve(const M &A, X &v, const Y &d, const K &w, BL< l > bl)
relaxed block upper triangular solve
Definition: gsetc.hh:263
void bdsolve(const M &A, X &v, const Y &d)
block diagonal solve, no relaxation
Definition: gsetc.hh:335
void bsorb(const M &A, X &x, const Y &b, const K &w, BL< l >)
Backward SOR step.
Definition: gsetc.hh:650
void bdsolve(const M &A, X &v, const Y &d, const K &w, BL< l >)
block diagonal solve, with relaxation
Definition: gsetc.hh:358
void butsolve(const M &A, X &v, const Y &d)
block upper triangular solve
Definition: gsetc.hh:200
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:632
void ubutsolve(const M &A, X &v, const Y &d)
unit block upper triangular solve
Definition: gsetc.hh:213
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:566
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
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:62
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:511
Dune namespace.
Definition: alignedallocator.hh:11
compile-time parameter for block recursion depth
Definition: gsetc.hh:43