5#ifndef DUNE_ISTL_GSETC_HH
6#define DUNE_ISTL_GSETC_HH
14#include <dune/common/hybridutilities.hh>
68 template<
int I, WithDiagType diag, WithRelaxType relax>
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);
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)
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& )
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)
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>
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());
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)
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>
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);
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>
569 typename... MultiTypeVectorArgs,
581 typename... MultiTypeVectorArgs,
593 typename... MultiTypeVectorArgs,
605 typename... MultiTypeVectorArgs,
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> )
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:568
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:484
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:540
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:64
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:513
void bltsolve(const M &A, X &v, const Y &d)
block lower triangular solve
Definition: gsetc.hh:175
WithDiagType
Definition: gsetc.hh:49
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
WithRelaxType
Definition: gsetc.hh:54
void bdsolve(const M &A, X &v, const Y &d)
block diagonal solve, no relaxation
Definition: gsetc.hh:337
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
@ nodiag
Definition: gsetc.hh:51
@ withdiag
Definition: gsetc.hh:50
@ norelax
Definition: gsetc.hh:56
@ withrelax
Definition: gsetc.hh:55
Definition: allocator.hh:11
A Vector class to support different block types.
Definition: multitypeblockvector.hh:59
A Matrix class to support different block types.
Definition: multitypeblockmatrix.hh:46
compile-time parameter for block recursion depth
Definition: gsetc.hh:45
@ recursion_level
Definition: gsetc.hh:46
static void butsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:90
static void bltsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:71
static void bltsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:114
static void butsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:120
static void bltsolve(const M &A, X &v, const Y &d, const K &)
Definition: gsetc.hh:129
static void butsolve(const M &A, X &v, const Y &d, const K &)
Definition: gsetc.hh:134
static void bltsolve(const M &, X &v, const Y &d, const K &w)
Definition: gsetc.hh:142
static void butsolve(const M &, X &v, const Y &d, const K &w)
Definition: gsetc.hh:148
static void bltsolve(const M &, X &v, const Y &d, const K &)
Definition: gsetc.hh:157
static void butsolve(const M &, X &v, const Y &d, const K &)
Definition: gsetc.hh:162
static void bdsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:296
static void bdsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:316
static void bdsolve(const M &A, X &v, const Y &d, const K &)
Definition: gsetc.hh:325
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:461
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:418
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:504
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:378
static void dbgs(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:545
static void dbjac(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:560
static void bsorf(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:550
static void bsorb(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:555
static void dbgs(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:571
static void dbjac(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:608
static void bsorf(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:583
static void bsorb(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:595