3#ifndef DUNE_ISTL_GSETC_HH
4#define DUNE_ISTL_GSETC_HH
11#include "multitypeblockvector.hh"
12#include "multitypeblockmatrix.hh"
14#include "istlexception.hh"
41 enum {recursion_level = l};
63 template<
int I, WithDiagType diag, WithRelaxType relax>
64 struct algmeta_btsolve {
65 template<
class M,
class X,
class Y,
class K>
66 static void bltsolve (
const M& A, X& v,
const Y& d,
const K& w)
69 typedef typename M::ConstRowIterator rowiterator;
70 typedef typename M::ConstColIterator coliterator;
71 typedef typename Y::block_type bblock;
74 rowiterator endi=A.end();
75 for (rowiterator i=A.begin(); i!=endi; ++i)
77 bblock rhs(d[i.index()]);
79 for (j=(*i).begin(); j.index()<i.index(); ++j)
80 (*j).mmv(v[j.index()],rhs);
84 template<
class M,
class X,
class Y,
class K>
85 static void butsolve (
const M& A, X& v,
const Y& d,
const K& w)
88 typedef typename M::ConstRowIterator rowiterator;
89 typedef typename M::ConstColIterator coliterator;
90 typedef typename Y::block_type bblock;
93 rowiterator rendi=A.beforeBegin();
94 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
96 bblock rhs(d[i.index()]);
98 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
99 (*j).mmv(v[j.index()],rhs);
107 struct algmeta_btsolve<0,withdiag,withrelax> {
108 template<
class M,
class X,
class Y,
class K>
109 static void bltsolve (
const M& A, X& v,
const Y& d,
const K& w)
114 template<
class M,
class X,
class Y,
class K>
115 static void butsolve (
const M& A, X& v,
const Y& d,
const K& w)
122 struct algmeta_btsolve<0,withdiag,norelax> {
123 template<
class M,
class X,
class Y,
class K>
124 static void bltsolve (
const M& A, X& v,
const Y& d,
const K& )
128 template<
class M,
class X,
class Y,
class K>
129 static void butsolve (
const M& A, X& v,
const Y& d,
const K& )
135 struct algmeta_btsolve<0,nodiag,withrelax> {
136 template<
class M,
class X,
class Y,
class K>
137 static void bltsolve (
const M& , X& v,
const Y& d,
const K& w)
142 template<
class M,
class X,
class Y,
class K>
143 static void butsolve (
const M& , X& v,
const Y& d,
const K& w)
150 struct algmeta_btsolve<0,nodiag,norelax> {
151 template<
class M,
class X,
class Y,
class K>
152 static void bltsolve (
const M& , X& v,
const Y& d,
const K& )
156 template<
class M,
class X,
class Y,
class K>
157 static void butsolve (
const M& , X& v,
const Y& d,
const K& )
169 template<
class M,
class X,
class Y>
172 typename X::field_type w=1;
176 template<
class M,
class X,
class Y,
class K>
177 void bltsolve (
const M& A, X& v,
const Y& d,
const K& w)
182 template<
class M,
class X,
class Y>
185 typename X::field_type w=1;
189 template<
class M,
class X,
class Y,
class K>
190 void ubltsolve (
const M& A, X& v,
const Y& d,
const K& w)
196 template<
class M,
class X,
class Y>
199 typename X::field_type w=1;
203 template<
class M,
class X,
class Y,
class K>
204 void butsolve (
const M& A, X& v,
const Y& d,
const K& w)
209 template<
class M,
class X,
class Y>
212 typename X::field_type w=1;
216 template<
class M,
class X,
class Y,
class K>
217 void ubutsolve (
const M& A, X& v,
const Y& d,
const K& w)
225 template<
class M,
class X,
class Y,
int l>
228 typename X::field_type w=1;
232 template<
class M,
class X,
class Y,
class K,
int l>
238 template<
class M,
class X,
class Y,
int l>
241 typename X::field_type w=1;
245 template<
class M,
class X,
class Y,
class K,
int l>
252 template<
class M,
class X,
class Y,
int l>
255 typename X::field_type w=1;
259 template<
class M,
class X,
class Y,
class K,
int l>
265 template<
class M,
class X,
class Y,
int l>
268 typename X::field_type w=1;
272 template<
class M,
class X,
class Y,
class K,
int l>
288 template<
int I, WithRelaxType relax>
289 struct algmeta_bdsolve {
290 template<
class M,
class X,
class Y,
class K>
291 static void bdsolve (
const M& A, X& v,
const Y& d,
const K& w)
294 typedef typename M::ConstRowIterator rowiterator;
295 typedef typename M::ConstColIterator coliterator;
298 rowiterator rendi=A.beforeBegin();
299 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
301 coliterator ii=(*i).find(i.index());
309 struct algmeta_bdsolve<0,withrelax> {
310 template<
class M,
class X,
class Y,
class K>
311 static void bdsolve (
const M& A, X& v,
const Y& d,
const K& w)
318 struct algmeta_bdsolve<0,norelax> {
319 template<
class M,
class X,
class Y,
class K>
320 static void bdsolve (
const M& A, X& v,
const Y& d,
const K& )
331 template<
class M,
class X,
class Y>
334 typename X::field_type w=1;
338 template<
class M,
class X,
class Y,
class K>
339 void bdsolve (
const M& A, X& v,
const Y& d,
const K& w)
347 template<
class M,
class X,
class Y,
int l>
350 typename X::field_type w=1;
354 template<
class M,
class X,
class Y,
class K,
int l>
369 template<
int I,
typename M>
370 struct algmeta_itsteps {
372 template<
class X,
class Y,
class K>
373 static void dbgs (
const M& A, X& x,
const Y& b,
const K& w)
375 typedef typename M::ConstRowIterator rowiterator;
376 typedef typename M::ConstColIterator coliterator;
377 typedef typename Y::block_type bblock;
382 rowiterator endi=A.end();
383 for (rowiterator i=A.begin(); i!=endi; ++i)
386 coliterator endj=(*i).end();
387 coliterator j=(*i).begin();
388 for (; j.index()<i.index(); ++j)
389 (*j).mmv(x[j.index()],rhs);
390 coliterator diag=j++;
391 for (; j != endj; ++j)
392 (*j).mmv(x[j.index()],rhs);
400 template<
class X,
class Y,
class K>
401 static void bsorf (
const M& A, X& x,
const Y& b,
const K& w)
403 typedef typename M::ConstRowIterator rowiterator;
404 typedef typename M::ConstColIterator coliterator;
405 typedef typename Y::block_type bblock;
406 typedef typename X::block_type xblock;
411 if(A.begin()!=A.end())
414 rowiterator endi=A.end();
415 for (rowiterator i=A.begin(); i!=endi; ++i)
418 coliterator endj=(*i).end();
419 coliterator j=(*i).begin();
420 for (; j.index()<i.index(); ++j)
421 (*j).mmv(x[j.index()],rhs);
424 (*j).mmv(x[j.index()],rhs);
426 x[i.index()].axpy(w,v);
430 template<
class X,
class Y,
class K>
431 static void bsorb (
const M& A, X& x,
const Y& b,
const K& w)
433 typedef typename M::ConstRowIterator rowiterator;
434 typedef typename M::ConstColIterator coliterator;
435 typedef typename Y::block_type bblock;
436 typedef typename X::block_type xblock;
441 if(A.begin()!=A.end())
444 rowiterator endi=A.beforeBegin();
445 for (rowiterator i=A.beforeEnd(); i!=endi; --i)
448 coliterator endj=(*i).end();
449 coliterator j=(*i).begin();
450 for (; j.index()<i.index(); ++j)
451 (*j).mmv(x[j.index()],rhs);
454 (*j).mmv(x[j.index()],rhs);
456 x[i.index()].axpy(w,v);
460 template<
class X,
class Y,
class K>
461 static void dbjac (
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;
470 rowiterator endi=A.end();
471 for (rowiterator i=A.begin(); i!=endi; ++i)
474 coliterator endj=(*i).end();
475 coliterator j=(*i).begin();
476 for (; j.index()<i.index(); ++j)
477 (*j).mmv(x[j.index()],rhs);
480 (*j).mmv(x[j.index()],rhs);
488 struct algmeta_itsteps<0,M> {
489 template<
class X,
class Y,
class K>
490 static void dbgs (
const M& A, X& x,
const Y& b,
const K& )
494 template<
class X,
class Y,
class K>
495 static void bsorf (
const M& A, X& x,
const Y& b,
const K& )
499 template<
class X,
class Y,
class K>
500 static void bsorb (
const M& A, X& x,
const Y& b,
const K& )
504 template<
class X,
class Y,
class K>
505 static void dbjac (
const M& A, X& x,
const Y& b,
const K& )
511 template<
int I,
typename T1,
typename... MultiTypeMatrixArgs>
512 struct algmeta_itsteps<I,MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>> {
514 typename... MultiTypeVectorArgs,
516 static void dbgs (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
517 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
518 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
526 typename... MultiTypeVectorArgs,
528 static void bsorf (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
529 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
530 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
538 typename... MultiTypeVectorArgs,
540 static void bsorb (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
541 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
542 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
550 typename... MultiTypeVectorArgs,
553 static void dbjac (
const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
554 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
555 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
566 template<
class M,
class X,
class Y,
class K>
567 void dbgs (
const M& A, X& x,
const Y& b,
const K& w)
572 template<
class M,
class X,
class Y,
class K,
int l>
573 void dbgs (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
578 template<
class M,
class X,
class Y,
class K>
579 void bsorf (
const M& A, X& x,
const Y& b,
const K& w)
584 template<
class M,
class X,
class Y,
class K,
int l>
585 void bsorf (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
590 template<
class M,
class X,
class Y,
class K>
591 void bsorb (
const M& A, X& x,
const Y& b,
const K& w)
596 template<
class M,
class X,
class Y,
class K,
int l>
597 void bsorb (
const M& A, X& x,
const Y& b,
const K& w,
BL<l> )
599 algmeta_itsteps<l,typename std::remove_cv<M>::type>
::bsorb(A,x,b,w);
602 template<
class M,
class X,
class Y,
class K>
603 void dbjac (
const M& A, X& x,
const Y& b,
const K& w)
608 template<
class M,
class X,
class Y,
class K,
int l>
609 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:170
void bsorb(const M &A, X &x, const Y &b, const K &w)
SSOR step.
Definition: gsetc.hh:591
void ubltsolve(const M &A, X &v, const Y &d)
unit block lower triangular solve
Definition: gsetc.hh:183
void dbjac(const M &A, X &x, const Y &b, const K &w)
Jacobi step.
Definition: gsetc.hh:603
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:567
void bsorf(const M &A, X &x, const Y &b, const K &w, BL< l >)
SOR step.
Definition: gsetc.hh:585
void dbjac(const M &A, X &x, const Y &b, const K &w, BL< l >)
Jacobi step.
Definition: gsetc.hh:609
void dbgs(const M &A, X &x, const Y &b, const K &w, BL< l >)
GS step.
Definition: gsetc.hh:573
void bltsolve(const M &A, X &v, const Y &d, const K &w, BL< l >)
relaxed block lower triangular solve
Definition: gsetc.hh:233
void butsolve(const M &A, X &v, const Y &d, const K &w, BL< l > bl)
relaxed block upper triangular solve
Definition: gsetc.hh:260
void bdsolve(const M &A, X &v, const Y &d)
block diagonal solve, no relaxation
Definition: gsetc.hh:332
void bsorb(const M &A, X &x, const Y &b, const K &w, BL< l >)
Backward SOR step.
Definition: gsetc.hh:597
void bdsolve(const M &A, X &v, const Y &d, const K &w, BL< l >)
block diagonal solve, with relaxation
Definition: gsetc.hh:355
void butsolve(const M &A, X &v, const Y &d)
block upper triangular solve
Definition: gsetc.hh:197
void bsorf(const M &A, X &x, const Y &b, const K &w)
SOR step.
Definition: gsetc.hh:579
void ubutsolve(const M &A, X &v, const Y &d)
unit block upper triangular solve
Definition: gsetc.hh:210
static void dbjac(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:329
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:245
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:301
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:274
static constexpr std::size_t N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:55
Dune namespace.
Definition: alignment.hh:11
compile-time parameter for block recursion depth
Definition: gsetc.hh:40