DUNE PDELab (2.7)

gsetc.hh
Go to the documentation of this file.
1// -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
2// vi: set et ts=4 sw=2 sts=2:
3#ifndef DUNE_ISTL_GSETC_HH
4#define DUNE_ISTL_GSETC_HH
5
6#include <cmath>
7#include <complex>
8#include <iostream>
9#include <iomanip>
10#include <string>
11
12#include <dune/common/hybridutilities.hh>
13
14#include "multitypeblockvector.hh"
15#include "multitypeblockmatrix.hh"
16
17#include "istlexception.hh"
18
19
25namespace Dune {
26
37 //============================================================
38 // parameter types
39 //============================================================
40
42 template<int l>
43 struct BL {
44 enum {recursion_level = l};
45 };
46
47 enum WithDiagType {
48 withdiag=1,
49 nodiag=0
50 };
51
52 enum WithRelaxType {
53 withrelax=1,
54 norelax=0
55 };
56
57 //============================================================
58 // generic triangular solves
59 // consider block decomposition A = L + D + U
60 // we can invert L, L+D, U, U+D
61 // we can apply relaxation or not
62 // we can recurse over a fixed number of levels
63 //============================================================
64
65 // template meta program for triangular solves
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)
70 {
71 // iterator types
72 typedef typename M::ConstRowIterator rowiterator;
73 typedef typename M::ConstColIterator coliterator;
74 typedef typename Y::block_type bblock;
75
76 // local solve at each block and immediate update
77 rowiterator endi=A.end();
78 for (rowiterator i=A.begin(); i!=endi; ++i)
79 {
80 bblock rhs(d[i.index()]);
81 coliterator j;
82 for (j=(*i).begin(); j.index()<i.index(); ++j)
83 (*j).mmv(v[j.index()],rhs);
84 algmeta_btsolve<I-1,diag,relax>::bltsolve(*j,v[i.index()],rhs,w);
85 }
86 }
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)
89 {
90 // iterator types
91 typedef typename M::ConstRowIterator rowiterator;
92 typedef typename M::ConstColIterator coliterator;
93 typedef typename Y::block_type bblock;
94
95 // local solve at each block and immediate update
96 rowiterator rendi=A.beforeBegin();
97 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
98 {
99 bblock rhs(d[i.index()]);
100 coliterator j;
101 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
102 (*j).mmv(v[j.index()],rhs);
103 algmeta_btsolve<I-1,diag,relax>::butsolve(*j,v[i.index()],rhs,w);
104 }
105 }
106 };
107
108 // recursion end ...
109 template<>
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)
113 {
114 A.solve(v,d);
115 v *= w;
116 }
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)
119 {
120 A.solve(v,d);
121 v *= w;
122 }
123 };
124 template<>
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& /*w*/)
128 {
129 A.solve(v,d);
130 }
131 template<class M, class X, class Y, class K>
132 static void butsolve (const M& A, X& v, const Y& d, const K& /*w*/)
133 {
134 A.solve(v,d);
135 }
136 };
137 template<>
138 struct algmeta_btsolve<0,nodiag,withrelax> {
139 template<class M, class X, class Y, class K>
140 static void bltsolve (const M& /*A*/, X& v, const Y& d, const K& w)
141 {
142 v = d;
143 v *= w;
144 }
145 template<class M, class X, class Y, class K>
146 static void butsolve (const M& /*A*/, X& v, const Y& d, const K& w)
147 {
148 v = d;
149 v *= w;
150 }
151 };
152 template<>
153 struct algmeta_btsolve<0,nodiag,norelax> {
154 template<class M, class X, class Y, class K>
155 static void bltsolve (const M& /*A*/, X& v, const Y& d, const K& /*w*/)
156 {
157 v = d;
158 }
159 template<class M, class X, class Y, class K>
160 static void butsolve (const M& /*A*/, X& v, const Y& d, const K& /*w*/)
161 {
162 v = d;
163 }
164 };
165
166
167 // user calls
168
169 // default block recursion level = 1
170
172 template<class M, class X, class Y>
173 void bltsolve (const M& A, X& v, const Y& d)
174 {
175 typename X::field_type w=1;
177 }
179 template<class M, class X, class Y, class K>
180 void bltsolve (const M& A, X& v, const Y& d, const K& w)
181 {
183 }
185 template<class M, class X, class Y>
186 void ubltsolve (const M& A, X& v, const Y& d)
187 {
188 typename X::field_type w=1;
190 }
192 template<class M, class X, class Y, class K>
193 void ubltsolve (const M& A, X& v, const Y& d, const K& w)
194 {
196 }
197
199 template<class M, class X, class Y>
200 void butsolve (const M& A, X& v, const Y& d)
201 {
202 typename X::field_type w=1;
204 }
206 template<class M, class X, class Y, class K>
207 void butsolve (const M& A, X& v, const Y& d, const K& w)
208 {
210 }
212 template<class M, class X, class Y>
213 void ubutsolve (const M& A, X& v, const Y& d)
214 {
215 typename X::field_type w=1;
217 }
219 template<class M, class X, class Y, class K>
220 void ubutsolve (const M& A, X& v, const Y& d, const K& w)
221 {
223 }
224
225 // general block recursion level >= 0
226
228 template<class M, class X, class Y, int l>
229 void bltsolve (const M& A, X& v, const Y& d, BL<l> /*bl*/)
230 {
231 typename X::field_type w=1;
233 }
235 template<class M, class X, class Y, class K, int l>
236 void bltsolve (const M& A, X& v, const Y& d, const K& w, BL<l> /*bl*/)
237 {
239 }
241 template<class M, class X, class Y, int l>
242 void ubltsolve (const M& A, X& v, const Y& d, BL<l> /*bl*/)
243 {
244 typename X::field_type w=1;
246 }
248 template<class M, class X, class Y, class K, int l>
249 void ubltsolve (const M& A, X& v, const Y& d, const K& w, BL<l> /*bl*/)
250 {
252 }
253
255 template<class M, class X, class Y, int l>
256 void butsolve (const M& A, X& v, const Y& d, BL<l> bl)
257 {
258 typename X::field_type w=1;
260 }
262 template<class M, class X, class Y, class K, int l>
263 void butsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
264 {
266 }
268 template<class M, class X, class Y, int l>
269 void ubutsolve (const M& A, X& v, const Y& d, BL<l> bl)
270 {
271 typename X::field_type w=1;
273 }
275 template<class M, class X, class Y, class K, int l>
276 void ubutsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
277 {
279 }
280
281
282
283 //============================================================
284 // generic block diagonal solves
285 // consider block decomposition A = L + D + U
286 // we can apply relaxation or not
287 // we can recurse over a fixed number of levels
288 //============================================================
289
290 // template meta program for diagonal solves
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)
295 {
296 // iterator types
297 typedef typename M::ConstRowIterator rowiterator;
298 typedef typename M::ConstColIterator coliterator;
299
300 // local solve at each block and immediate update
301 rowiterator rendi=A.beforeBegin();
302 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
303 {
304 coliterator ii=(*i).find(i.index());
305 algmeta_bdsolve<I-1,relax>::bdsolve(*ii,v[i.index()],d[i.index()],w);
306 }
307 }
308 };
309
310 // recursion end ...
311 template<>
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)
315 {
316 A.solve(v,d);
317 v *= w;
318 }
319 };
320 template<>
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& /*w*/)
324 {
325 A.solve(v,d);
326 }
327 };
328
329 // user calls
330
331 // default block recursion level = 1
332
334 template<class M, class X, class Y>
335 void bdsolve (const M& A, X& v, const Y& d)
336 {
337 typename X::field_type w=1;
339 }
341 template<class M, class X, class Y, class K>
342 void bdsolve (const M& A, X& v, const Y& d, const K& w)
343 {
345 }
346
347 // general block recursion level >= 0
348
350 template<class M, class X, class Y, int l>
351 void bdsolve (const M& A, X& v, const Y& d, BL<l> /*bl*/)
352 {
353 typename X::field_type w=1;
355 }
357 template<class M, class X, class Y, class K, int l>
358 void bdsolve (const M& A, X& v, const Y& d, const K& w, BL<l> /*bl*/)
359 {
361 }
362
363
364 //============================================================
365 // generic steps of iteration methods
366 // Jacobi, Gauss-Seidel, SOR, SSOR
367 // work directly on Ax=b, ie solve M(x^{i+1}-x^i) = w (b-Ax^i)
368 // we can recurse over a fixed number of levels
369 //============================================================
370
371 // template meta program for iterative solver steps
372 template<int I, typename M>
373 struct algmeta_itsteps {
374
375 template<class X, class Y, class K>
376 static void dbgs (const M& A, X& x, const Y& b, const K& w)
377 {
378 typedef typename M::ConstRowIterator rowiterator;
379 typedef typename M::ConstColIterator coliterator;
380 typedef typename Y::block_type bblock;
381 bblock rhs;
382
383 X xold(x); // remember old x
384
385 rowiterator endi=A.end();
386 for (rowiterator i=A.begin(); i!=endi; ++i)
387 {
388 rhs = b[i.index()]; // rhs = b_i
389 coliterator endj=(*i).end();
390 coliterator j=(*i).begin();
391 Hybrid::ifElse(IsNumber<typename M::block_type>(),
392 [&](auto id) {
393 for (; j.index()<i.index(); ++j)
394 rhs -= id(*j) * x[j.index()];
395 coliterator diag=j++; // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
396 for (; j != endj; ++j)
397 rhs -= id(*j) * x[j.index()];
398 x[i.index()] = rhs / id(*diag);
399 },
400 [&](auto id) {
401 for (; j.index()<i.index(); ++j) // iterate over a_ij with j < i
402 id(*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
403 coliterator diag=j++; // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
404 for (; j != endj; ++j) // iterate over a_ij with j > i
405 id(*j).mmv(x[j.index()],rhs); // rhs -= sum_{j>i} a_ij * xold_j
406 algmeta_itsteps<I-1,typename M::block_type>::dbgs(*diag,x[i.index()],rhs,w); // if I==1: xnew_i = rhs/a_ii
407 });
408 }
409 // next two lines: xnew_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j) + (1-w)*xold;
410 x *= w;
411 x.axpy(K(1)-w,xold);
412 }
413
414 template<class X, class Y, class K>
415 static void bsorf (const M& A, X& x, const Y& b, const K& w)
416 {
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;
421 bblock rhs;
422 xblock v;
423
424 // Initialize nested data structure if there are entries
425 if(A.begin()!=A.end())
426 v=x[0];
427
428 rowiterator endi=A.end();
429 for (rowiterator i=A.begin(); i!=endi; ++i)
430 {
431 rhs = b[i.index()]; // rhs = b_i
432 coliterator endj=(*i).end(); // iterate over a_ij with j < i
433 coliterator j=(*i).begin();
434 Hybrid::ifElse(IsNumber<typename M::block_type>(),
435 [&](auto id) {
436 for (; j.index()<i.index(); ++j)
437 rhs -= id(*j) * x[j.index()]; // rhs -= sum_{j<i} a_ij * xnew_j
438 coliterator diag=j; // *diag = a_ii
439 for (; j!=endj; ++j)
440 rhs -= id(*j) * x[j.index()]; // rhs -= sum_{j<i} a_ij * xnew_j
441 v = rhs / id(*diag);
442 id(x[i.index()]) += w*id(v); // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
443 },
444 [&](auto id) {
445 for (; j.index()<i.index(); ++j)
446 id(*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
447 coliterator diag=j; // *diag = a_ii
448 for (; j!=endj; ++j)
449 id(*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
450 algmeta_itsteps<I-1,typename M::block_type>::bsorf(*diag,v,rhs,w); // if blocksize I==1: v = rhs/a_ii
451 id(x[i.index()]).axpy(w,v); // x_i = w / a_ii * (b_i - sum_{j<i} a_ij * xnew_j - sum_{j>=i} a_ij * xold_j)
452 });
453 }
454 }
455
456 template<class X, class Y, class K>
457 static void bsorb (const M& A, X& x, const Y& b, const K& w)
458 {
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;
463 bblock rhs;
464 xblock v;
465
466 // Initialize nested data structure if there are entries
467 if(A.begin()!=A.end())
468 v=x[0];
469
470 rowiterator endi=A.beforeBegin();
471 for (rowiterator i=A.beforeEnd(); i!=endi; --i)
472 {
473 rhs = b[i.index()];
474 coliterator endj=(*i).end();
475 coliterator j=(*i).begin();
476 Hybrid::ifElse(IsNumber<typename M::block_type>(),
477 [&](auto id) {
478 for (; j.index()<i.index(); ++j)
479 rhs -= id(*j) * x[j.index()];
480 coliterator diag=j;
481 for (; j!=endj; ++j)
482 rhs -= id(*j) * x[j.index()];
483 v = rhs / id(*diag);
484 x[i.index()] += w*id(v);
485 },
486 [&](auto id) {
487 for (; j.index()<i.index(); ++j)
488 id(*j).mmv(x[j.index()],rhs);
489 coliterator diag=j;
490 for (; j!=endj; ++j)
491 id(*j).mmv(x[j.index()],rhs);
493 id(x[i.index()]).axpy(w,v);
494 });
495 }
496 }
497
498 template<class X, class Y, class K>
499 static void dbjac (const M& A, X& x, const Y& b, const K& w)
500 {
501 typedef typename M::ConstRowIterator rowiterator;
502 typedef typename M::ConstColIterator coliterator;
503 typedef typename Y::block_type bblock;
504 bblock rhs;
505
506 X v(x); // allocate with same size
507
508 rowiterator endi=A.end();
509 for (rowiterator i=A.begin(); i!=endi; ++i)
510 {
511 rhs = b[i.index()];
512 coliterator endj=(*i).end();
513 coliterator j=(*i).begin();
514 Hybrid::ifElse(IsNumber<typename M::block_type>(),
515 [&](auto id) {
516 for (; j.index()<i.index(); ++j)
517 rhs -= id(*j) * x[j.index()];
518 coliterator diag=j;
519 for (; j!=endj; ++j)
520 rhs -= id(*j) * x[j.index()];
521 v[i.index()] = rhs / id(*diag);
522 },
523 [&](auto id) {
524 for (; j.index()<i.index(); ++j)
525 id(*j).mmv(x[j.index()],rhs);
526 coliterator diag=j;
527 for (; j!=endj; ++j)
528 id(*j).mmv(x[j.index()],rhs);
530 });
531 }
532 x.axpy(w,v);
533 }
534 };
535 // end of recursion
536 template<typename M>
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& /*w*/)
540 {
541 A.solve(x,b);
542 }
543 template<class X, class Y, class K>
544 static void bsorf (const M& A, X& x, const Y& b, const K& /*w*/)
545 {
546 A.solve(x,b);
547 }
548 template<class X, class Y, class K>
549 static void bsorb (const M& A, X& x, const Y& b, const K& /*w*/)
550 {
551 A.solve(x,b);
552 }
553 template<class X, class Y, class K>
554 static void dbjac (const M& A, X& x, const Y& b, const K& /*w*/)
555 {
556 A.solve(x,b);
557 }
558 };
559
560 template<int I, typename T1, typename... MultiTypeMatrixArgs>
561 struct algmeta_itsteps<I,MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>> {
562 template<
563 typename... MultiTypeVectorArgs,
564 class K>
565 static void dbgs (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
566 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
567 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
568 const K& w)
569 {
572 }
573
574 template<
575 typename... MultiTypeVectorArgs,
576 class K>
577 static void bsorf (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
578 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
579 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
580 const K& w)
581 {
584 }
585
586 template<
587 typename... MultiTypeVectorArgs,
588 class K>
589 static void bsorb (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
590 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
591 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
592 const K& w)
593 {
596 }
597
598 template<
599 typename... MultiTypeVectorArgs,
600 class K
601 >
602 static void dbjac (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
603 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
604 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
605 const K& w)
606 {
609 }
610 };
611
612 // user calls
613
615 template<class M, class X, class Y, class K>
616 void dbgs (const M& A, X& x, const Y& b, const K& w)
617 {
619 }
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> /*bl*/)
623 {
625 }
627 template<class M, class X, class Y, class K>
628 void bsorf (const M& A, X& x, const Y& b, const K& w)
629 {
631 }
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> /*bl*/)
635 {
637 }
639 template<class M, class X, class Y, class K>
640 void bsorb (const M& A, X& x, const Y& b, const K& w)
641 {
643 }
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> /*bl*/)
647 {
648 algmeta_itsteps<l,typename std::remove_cv<M>::type>::bsorb(A,x,b,w);
649 }
651 template<class M, class X, class Y, class K>
652 void dbjac (const M& A, X& x, const Y& b, const K& w)
653 {
655 }
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> /*bl*/)
659 {
661 }
662
663
666} // end namespace
667
668#endif
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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Jul 15, 22:36, 2024)