5#ifndef DUNE_ISTL_GSETC_HH
6#define DUNE_ISTL_GSETC_HH
14#include <dune/common/hybridutilities.hh>
16#include "multitypeblockvector.hh"
17#include "multitypeblockmatrix.hh"
19#include "istlexception.hh"
46 enum {recursion_level = l};
68 template<
int I, WithDiagType diag, WithRelaxType relax>
69 struct algmeta_btsolve {
70 template<
class M,
class X,
class Y,
class K>
71 static void bltsolve (
const M& A, X& v,
const Y& d,
const K& w)
74 typedef typename M::ConstRowIterator rowiterator;
75 typedef typename M::ConstColIterator coliterator;
76 typedef typename Y::block_type bblock;
79 rowiterator endi=A.end();
80 for (rowiterator i=A.begin(); i!=endi; ++i)
82 bblock rhs(d[i.index()]);
84 for (j=(*i).begin(); j.index()<i.index(); ++j)
85 (*j).mmv(v[j.index()],rhs);
89 template<
class M,
class X,
class Y,
class K>
90 static void butsolve (
const M& A, X& v,
const Y& d,
const K& w)
93 typedef typename M::ConstRowIterator rowiterator;
94 typedef typename M::ConstColIterator coliterator;
95 typedef typename Y::block_type bblock;
98 rowiterator rendi=A.beforeBegin();
99 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
101 bblock rhs(d[i.index()]);
103 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
104 (*j).mmv(v[j.index()],rhs);
112 struct algmeta_btsolve<0,withdiag,withrelax> {
113 template<
class M,
class X,
class Y,
class K>
114 static void bltsolve (
const M& A, X& v,
const Y& d,
const K& w)
119 template<
class M,
class X,
class Y,
class K>
120 static void butsolve (
const M& A, X& v,
const Y& d,
const K& w)
127 struct algmeta_btsolve<0,withdiag,norelax> {
128 template<
class M,
class X,
class Y,
class K>
129 static void bltsolve (
const M& A, X& v,
const Y& d,
const K& )
133 template<
class M,
class X,
class Y,
class K>
134 static void butsolve (
const M& A, X& v,
const Y& d,
const K& )
140 struct algmeta_btsolve<0,nodiag,withrelax> {
141 template<
class M,
class X,
class Y,
class K>
142 static void bltsolve (
const M& , X& v,
const Y& d,
const K& w)
147 template<
class M,
class X,
class Y,
class K>
148 static void butsolve (
const M& , X& v,
const Y& d,
const K& w)
155 struct algmeta_btsolve<0,nodiag,norelax> {
156 template<
class M,
class X,
class Y,
class K>
157 static void bltsolve (
const M& , X& v,
const Y& d,
const K& )
161 template<
class M,
class X,
class Y,
class K>
162 static void butsolve (
const M& , X& v,
const Y& d,
const K& )
174 template<
class M,
class X,
class Y>
177 typename X::field_type w=1;
181 template<
class M,
class X,
class Y,
class K>
182 void bltsolve (
const M& A, X& v,
const Y& d,
const K& w)
187 template<
class M,
class X,
class Y>
190 typename X::field_type w=1;
194 template<
class M,
class X,
class Y,
class K>
195 void ubltsolve (
const M& A, X& v,
const Y& d,
const K& w)
201 template<
class M,
class X,
class Y>
204 typename X::field_type w=1;
208 template<
class M,
class X,
class Y,
class K>
209 void butsolve (
const M& A, X& v,
const Y& d,
const K& w)
214 template<
class M,
class X,
class Y>
217 typename X::field_type w=1;
221 template<
class M,
class X,
class Y,
class K>
222 void ubutsolve (
const M& A, X& v,
const Y& d,
const K& w)
230 template<
class M,
class X,
class Y,
int l>
233 typename X::field_type w=1;
237 template<
class M,
class X,
class Y,
class K,
int l>
243 template<
class M,
class X,
class Y,
int l>
246 typename X::field_type w=1;
250 template<
class M,
class X,
class Y,
class K,
int l>
257 template<
class M,
class X,
class Y,
int l>
260 typename X::field_type w=1;
264 template<
class M,
class X,
class Y,
class K,
int l>
270 template<
class M,
class X,
class Y,
int l>
273 typename X::field_type w=1;
277 template<
class M,
class X,
class Y,
class K,
int l>
293 template<
int I, WithRelaxType relax>
294 struct algmeta_bdsolve {
295 template<
class M,
class X,
class Y,
class K>
296 static void bdsolve (
const M& A, X& v,
const Y& d,
const K& w)
299 typedef typename M::ConstRowIterator rowiterator;
300 typedef typename M::ConstColIterator coliterator;
303 rowiterator rendi=A.beforeBegin();
304 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
306 coliterator ii=(*i).find(i.index());
314 struct algmeta_bdsolve<0,withrelax> {
315 template<
class M,
class X,
class Y,
class K>
316 static void bdsolve (
const M& A, X& v,
const Y& d,
const K& w)
323 struct algmeta_bdsolve<0,norelax> {
324 template<
class M,
class X,
class Y,
class K>
325 static void bdsolve (
const M& A, X& v,
const Y& d,
const K& )
336 template<
class M,
class X,
class Y>
339 typename X::field_type w=1;
343 template<
class M,
class X,
class Y,
class K>
344 void bdsolve (
const M& A, X& v,
const Y& d,
const K& w)
352 template<
class M,
class X,
class Y,
int l>
355 typename X::field_type w=1;
359 template<
class M,
class X,
class Y,
class K,
int l>
374 template<
int I,
typename M>
375 struct algmeta_itsteps {
377 template<
class X,
class Y,
class K>
378 static void dbgs (
const M& A, X& x,
const Y& b,
const K& w)
380 typedef typename M::ConstRowIterator rowiterator;
381 typedef typename M::ConstColIterator coliterator;
382 typedef typename Y::block_type bblock;
387 rowiterator endi=A.end();
388 for (rowiterator i=A.begin(); i!=endi; ++i)
391 coliterator endj=(*i).end();
392 coliterator j=(*i).begin();
393 if constexpr (IsNumber<typename M::block_type>())
395 for (; j.index()<i.index(); ++j)
396 rhs -= (*j) * x[j.index()];
397 coliterator diag=j++;
398 for (; j != endj; ++j)
399 rhs -= (*j) * x[j.index()];
400 x[i.index()] = rhs / (*diag);
404 for (; j.index()<i.index(); ++j)
405 (*j).mmv(x[j.index()],rhs);
406 coliterator diag=j++;
407 for (; j != endj; ++j)
408 (*j).mmv(x[j.index()],rhs);
417 template<
class X,
class Y,
class K>
418 static void bsorf (
const M& A, X& x,
const Y& b,
const K& w)
420 typedef typename M::ConstRowIterator rowiterator;
421 typedef typename M::ConstColIterator coliterator;
422 typedef typename Y::block_type bblock;
423 typedef typename X::block_type xblock;
428 if(A.begin()!=A.end())
431 rowiterator endi=A.end();
432 for (rowiterator i=A.begin(); i!=endi; ++i)
435 coliterator endj=(*i).end();
436 coliterator j=(*i).begin();
437 if constexpr (IsNumber<typename M::block_type>())
439 for (; j.index()<i.index(); ++j)
440 rhs -= (*j) * x[j.index()];
443 rhs -= (*j) * x[j.index()];
449 for (; j.index()<i.index(); ++j)
450 (*j).mmv(x[j.index()],rhs);
453 (*j).mmv(x[j.index()],rhs);
455 x[i.index()].axpy(w,v);
460 template<
class X,
class Y,
class K>
461 static void bsorb (
const M& A, X& x,
const Y& b,
const K& w)
463 typedef typename M::ConstRowIterator rowiterator;
464 typedef typename M::ConstColIterator coliterator;
465 typedef typename Y::block_type bblock;
466 typedef typename X::block_type xblock;
471 if(A.begin()!=A.end())
474 rowiterator endi=A.beforeBegin();
475 for (rowiterator i=A.beforeEnd(); i!=endi; --i)
478 coliterator endj=(*i).end();
479 coliterator j=(*i).begin();
480 if constexpr (IsNumber<typename M::block_type>())
482 for (; j.index()<i.index(); ++j)
483 rhs -= (*j) * x[j.index()];
486 rhs -= (*j) * x[j.index()];
492 for (; j.index()<i.index(); ++j)
493 j->mmv(x[j.index()],rhs);
496 j->mmv(x[j.index()],rhs);
498 x[i.index()].axpy(w,v);
503 template<
class X,
class Y,
class K>
504 static void dbjac (
const M& A, X& x,
const Y& b,
const K& w)
506 typedef typename M::ConstRowIterator rowiterator;
507 typedef typename M::ConstColIterator coliterator;
508 typedef typename Y::block_type bblock;
513 rowiterator endi=A.end();
514 for (rowiterator i=A.begin(); i!=endi; ++i)
517 coliterator endj=(*i).end();
518 coliterator j=(*i).begin();
519 if constexpr (IsNumber<typename M::block_type>())
521 for (; j.index()<i.index(); ++j)
522 rhs -= (*j) * x[j.index()];
525 rhs -= (*j) * x[j.index()];
526 v[i.index()] = rhs / (*diag);
530 for (; j.index()<i.index(); ++j)
531 j->mmv(x[j.index()],rhs);
534 j->mmv(x[j.index()],rhs);
543 struct algmeta_itsteps<0,M> {
544 template<
class X,
class Y,
class K>
545 static void dbgs (
const M& A, X& x,
const Y& b,
const K& )
549 template<
class X,
class Y,
class K>
550 static void bsorf (
const M& A, X& x,
const Y& b,
const K& )
554 template<
class X,
class Y,
class K>
555 static void bsorb (
const M& A, X& x,
const Y& b,
const K& )
559 template<
class X,
class Y,
class K>
560 static void dbjac (
const M& A, X& x,
const Y& b,
const K& )
566 template<
int I,
typename T1,
typename... MultiTypeMatrixArgs>
567 struct algmeta_itsteps<I,MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>> {
569 typename... MultiTypeVectorArgs,
571 static void dbgs (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
572 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
573 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
581 typename... MultiTypeVectorArgs,
583 static void bsorf (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
584 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
585 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
593 typename... MultiTypeVectorArgs,
595 static void bsorb (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
596 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
597 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
605 typename... MultiTypeVectorArgs,
608 static void dbjac (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
609 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
610 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
621 template<
class M,
class X,
class Y,
class K>
622 void dbgs (
const M& A, X& x,
const Y& b,
const K& w)
627 template<
class M,
class X,
class Y,
class K,
int l>
628 void dbgs (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
633 template<
class M,
class X,
class Y,
class K>
634 void bsorf (
const M& A, X& x,
const Y& b,
const K& w)
639 template<
class M,
class X,
class Y,
class K,
int l>
640 void bsorf (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
645 template<
class M,
class X,
class Y,
class K>
646 void bsorb (
const M& A, X& x,
const Y& b,
const K& w)
651 template<
class M,
class X,
class Y,
class K,
int l>
652 void bsorb (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
654 algmeta_itsteps<l,typename std::remove_cv<M>::type>
::bsorb(A,x,b,w);
657 template<
class M,
class X,
class Y,
class K>
658 void dbjac (
const M& A, X& x,
const Y& b,
const K& w)
663 template<
class M,
class X,
class Y,
class K,
int l>
664 void dbjac (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:584
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
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:84
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:529
void bltsolve(const M &A, X &v, const Y &d)
block lower triangular solve
Definition: gsetc.hh:175
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:646
void ubltsolve(const M &A, X &v, const Y &d)
unit block lower triangular solve
Definition: gsetc.hh:188
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, BL< l >)
SOR step.
Definition: gsetc.hh:640
void dbjac(const M &A, X &x, const Y &b, const K &w, BL< l >)
Jacobi step.
Definition: gsetc.hh:664
void dbgs(const M &A, X &x, const Y &b, const K &w, BL< l >)
GS step.
Definition: gsetc.hh:628
void bltsolve(const M &A, X &v, const Y &d, const K &w, BL< l >)
relaxed block lower triangular solve
Definition: gsetc.hh:238
void butsolve(const M &A, X &v, const Y &d, const K &w, BL< l > bl)
relaxed block upper triangular solve
Definition: gsetc.hh:265
void bdsolve(const M &A, X &v, const Y &d)
block diagonal solve, no relaxation
Definition: gsetc.hh:337
void bsorb(const M &A, X &x, const Y &b, const K &w, BL< l >)
Backward SOR step.
Definition: gsetc.hh:652
void bdsolve(const M &A, X &v, const Y &d, const K &w, BL< l >)
block diagonal solve, with relaxation
Definition: gsetc.hh:360
void butsolve(const M &A, X &v, const Y &d)
block upper triangular solve
Definition: gsetc.hh:202
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:634
void ubutsolve(const M &A, X &v, const Y &d)
unit block upper triangular solve
Definition: gsetc.hh:215
Dune namespace.
Definition: alignedallocator.hh:13
compile-time parameter for block recursion depth
Definition: gsetc.hh:45