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();
393 for (; j.index()<i.index(); ++j)
394 rhs -=
id(*j) * x[j.index()];
395 coliterator diag=j++;
396 for (; j != endj; ++j)
397 rhs -=
id(*j) * x[j.index()];
398 x[i.index()] = rhs /
id(*diag);
401 for (; j.index()<i.index(); ++j)
402 id(*j).mmv(x[j.index()],rhs);
403 coliterator diag=j++;
404 for (; j != endj; ++j)
405 id(*j).mmv(x[j.index()],rhs);
414 template<
class X,
class Y,
class K>
415 static void bsorf (
const M& A, X& x,
const Y& b,
const K& w)
417 typedef typename M::ConstRowIterator rowiterator;
418 typedef typename M::ConstColIterator coliterator;
419 typedef typename Y::block_type bblock;
420 typedef typename X::block_type xblock;
425 if(A.begin()!=A.end())
428 rowiterator endi=A.end();
429 for (rowiterator i=A.begin(); i!=endi; ++i)
432 coliterator endj=(*i).end();
433 coliterator j=(*i).begin();
436 for (; j.index()<i.index(); ++j)
437 rhs -=
id(*j) * x[j.index()];
440 rhs -=
id(*j) * x[j.index()];
442 id(x[i.index()]) += w*
id(v);
445 for (; j.index()<i.index(); ++j)
446 id(*j).mmv(x[j.index()],rhs);
449 id(*j).mmv(x[j.index()],rhs);
451 id(x[i.index()]).axpy(w,v);
456 template<
class X,
class Y,
class K>
457 static void bsorb (
const M& A, X& x,
const Y& b,
const K& w)
459 typedef typename M::ConstRowIterator rowiterator;
460 typedef typename M::ConstColIterator coliterator;
461 typedef typename Y::block_type bblock;
462 typedef typename X::block_type xblock;
467 if(A.begin()!=A.end())
470 rowiterator endi=A.beforeBegin();
471 for (rowiterator i=A.beforeEnd(); i!=endi; --i)
474 coliterator endj=(*i).end();
475 coliterator j=(*i).begin();
478 for (; j.index()<i.index(); ++j)
479 rhs -=
id(*j) * x[j.index()];
482 rhs -=
id(*j) * x[j.index()];
484 x[i.index()] += w*
id(v);
487 for (; j.index()<i.index(); ++j)
488 id(*j).mmv(x[j.index()],rhs);
491 id(*j).mmv(x[j.index()],rhs);
493 id(x[i.index()]).axpy(w,v);
498 template<
class X,
class Y,
class K>
499 static void dbjac (
const M& A, X& x,
const Y& b,
const K& w)
501 typedef typename M::ConstRowIterator rowiterator;
502 typedef typename M::ConstColIterator coliterator;
503 typedef typename Y::block_type bblock;
508 rowiterator endi=A.end();
509 for (rowiterator i=A.begin(); i!=endi; ++i)
512 coliterator endj=(*i).end();
513 coliterator j=(*i).begin();
516 for (; j.index()<i.index(); ++j)
517 rhs -=
id(*j) * x[j.index()];
520 rhs -=
id(*j) * x[j.index()];
521 v[i.index()] = rhs /
id(*diag);
524 for (; j.index()<i.index(); ++j)
525 id(*j).mmv(x[j.index()],rhs);
528 id(*j).mmv(x[j.index()],rhs);
537 struct algmeta_itsteps<0,M> {
538 template<
class X,
class Y,
class K>
539 static void dbgs (
const M& A, X& x,
const Y& b,
const K& )
543 template<
class X,
class Y,
class K>
544 static void bsorf (
const M& A, X& x,
const Y& b,
const K& )
548 template<
class X,
class Y,
class K>
549 static void bsorb (
const M& A, X& x,
const Y& b,
const K& )
553 template<
class X,
class Y,
class K>
554 static void dbjac (
const M& A, X& x,
const Y& b,
const K& )
560 template<
int I,
typename T1,
typename... MultiTypeMatrixArgs>
561 struct algmeta_itsteps<I,MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>> {
563 typename... MultiTypeVectorArgs,
565 static void dbgs (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
566 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
567 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
575 typename... MultiTypeVectorArgs,
577 static void bsorf (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
578 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
579 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
587 typename... MultiTypeVectorArgs,
589 static void bsorb (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
590 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
591 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
599 typename... MultiTypeVectorArgs,
602 static void dbjac (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
603 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
604 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
615 template<
class M,
class X,
class Y,
class K>
616 void dbgs (
const M& A, X& x,
const Y& b,
const K& w)
621 template<
class M,
class X,
class Y,
class K,
int l>
622 void dbgs (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
627 template<
class M,
class X,
class Y,
class K>
628 void bsorf (
const M& A, X& x,
const Y& b,
const K& w)
633 template<
class M,
class X,
class Y,
class K,
int l>
634 void bsorf (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
639 template<
class M,
class X,
class Y,
class K>
640 void bsorb (
const M& A, X& x,
const Y& b,
const K& w)
645 template<
class M,
class X,
class Y,
class K,
int l>
646 void bsorb (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
648 algmeta_itsteps<l,typename std::remove_cv<M>::type>
::bsorb(A,x,b,w);
651 template<
class M,
class X,
class Y,
class K>
652 void dbjac (
const M& A, X& x,
const Y& b,
const K& w)
657 template<
class M,
class X,
class Y,
class K,
int l>
658 void dbjac (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
decltype(auto) ifElse(const Condition &condition, IfFunc &&ifFunc, ElseFunc &&elseFunc)
A conditional expression.
Definition: hybridutilities.hh:355
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:640
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:652
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:616
void bsorf(const M &A, X &x, const Y &b, const K &w, BL< l >)
SOR step.
Definition: gsetc.hh:634
void dbjac(const M &A, X &x, const Y &b, const K &w, BL< l >)
Jacobi step.
Definition: gsetc.hh:658
void dbgs(const M &A, X &x, const Y &b, const K &w, BL< l >)
GS step.
Definition: gsetc.hh:622
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:646
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:628
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:557
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:473
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:529
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:58
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:502
Dune namespace.
Definition: alignedallocator.hh:14
compile-time parameter for block recursion depth
Definition: gsetc.hh:43