Go to the documentation of this file.
3 #ifndef DUNE_ISTL_GSETC_HH
4 #define DUNE_ISTL_GSETC_HH
12 #include <dune/common/hybridutilities.hh>
66 template<
int I, WithDiagType diag, WithRelaxType relax>
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);
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)
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& )
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)
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>
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());
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)
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>
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 Hybrid::ifElse(IsNumber<typename M::block_type>(),
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();
434 Hybrid::ifElse(IsNumber<typename M::block_type>(),
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();
476 Hybrid::ifElse(IsNumber<typename M::block_type>(),
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();
514 Hybrid::ifElse(IsNumber<typename M::block_type>(),
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);
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>
563 typename... MultiTypeVectorArgs,
575 typename... MultiTypeVectorArgs,
587 typename... MultiTypeVectorArgs,
599 typename... MultiTypeVectorArgs,
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> )
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> )
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:502
void bdsolve(const M &A, X &v, const Y &d)
block diagonal solve, no relaxation
Definition: gsetc.hh:335
@ nodiag
Definition: gsetc.hh:49
@ norelax
Definition: gsetc.hh:54
static void dbjac(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:602
static void bltsolve(const M &, X &v, const Y &d, const K &)
Definition: gsetc.hh:155
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
static void butsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:118
A Matrix class to support different block types.
Definition: matrixutils.hh:87
static void bsorb(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:549
WithRelaxType
Definition: gsetc.hh:52
static void bsorb(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:457
@ recursion_level
Definition: gsetc.hh:44
WithDiagType
Definition: gsetc.hh:47
static void bdsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:294
static void bdsolve(const M &A, X &v, const Y &d, const K &)
Definition: gsetc.hh:323
static void bltsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:69
static void bltsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:112
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:529
static void dbjac(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:554
static void bsorf(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:577
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:58
A Vector class to support different block types.
Definition: multitypeblockvector.hh:20
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:557
Definition: allocator.hh:7
static void bltsolve(const M &, X &v, const Y &d, const K &w)
Definition: gsetc.hh:140
static void bltsolve(const M &A, X &v, const Y &d, const K &)
Definition: gsetc.hh:127
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:652
static void bsorb(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:589
static void bdsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:314
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:473
void ubutsolve(const M &A, X &v, const Y &d)
unit block upper triangular solve
Definition: gsetc.hh:213
static void bsorf(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:544
void ubltsolve(const M &A, X &v, const Y &d)
unit block lower triangular solve
Definition: gsetc.hh:186
static void dbgs(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:376
void butsolve(const M &A, X &v, const Y &d)
block upper triangular solve
Definition: gsetc.hh:200
compile-time parameter for block recursion depth
Definition: gsetc.hh:43
static void dbgs(const MultiTypeBlockMatrix< T1, MultiTypeMatrixArgs... > &A, MultiTypeBlockVector< MultiTypeVectorArgs... > &x, const MultiTypeBlockVector< MultiTypeVectorArgs... > &b, const K &w)
Definition: gsetc.hh:565
@ withrelax
Definition: gsetc.hh:53
static void butsolve(const M &A, X &v, const Y &d, const K &w)
Definition: gsetc.hh:88
static void dbgs(const M &A, X &x, const Y &b, const K &)
Definition: gsetc.hh:539
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:616
@ withdiag
Definition: gsetc.hh:48
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:628
static void bsorf(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:415
static void butsolve(const M &A, X &v, const Y &d, const K &)
Definition: gsetc.hh:132
static void dbjac(const M &A, X &x, const Y &b, const K &w)
Definition: gsetc.hh:499
static void butsolve(const M &, X &v, const Y &d, const K &w)
Definition: gsetc.hh:146
static void butsolve(const M &, X &v, const Y &d, const K &)
Definition: gsetc.hh:160