Dune Core Modules (2.4.1)

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#include "multitypeblockvector.hh"
12#include "multitypeblockmatrix.hh"
13
14#include "istlexception.hh"
15
16
22namespace Dune {
23
34 //============================================================
35 // parameter types
36 //============================================================
37
39 template<int l>
40 struct BL {
41 enum {recursion_level = l};
42 };
43
44 enum WithDiagType {
45 withdiag=1,
46 nodiag=0
47 };
48
49 enum WithRelaxType {
50 withrelax=1,
51 norelax=0
52 };
53
54 //============================================================
55 // generic triangular solves
56 // consider block decomposition A = L + D + U
57 // we can invert L, L+D, U, U+D
58 // we can apply relaxation or not
59 // we can recurse over a fixed number of levels
60 //============================================================
61
62 // template meta program for triangular solves
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)
67 {
68 // iterator types
69 typedef typename M::ConstRowIterator rowiterator;
70 typedef typename M::ConstColIterator coliterator;
71 typedef typename Y::block_type bblock;
72
73 // local solve at each block and immediate update
74 rowiterator endi=A.end();
75 for (rowiterator i=A.begin(); i!=endi; ++i)
76 {
77 bblock rhs(d[i.index()]);
78 coliterator j;
79 for (j=(*i).begin(); j.index()<i.index(); ++j)
80 (*j).mmv(v[j.index()],rhs);
81 algmeta_btsolve<I-1,diag,relax>::bltsolve(*j,v[i.index()],rhs,w);
82 }
83 }
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)
86 {
87 // iterator types
88 typedef typename M::ConstRowIterator rowiterator;
89 typedef typename M::ConstColIterator coliterator;
90 typedef typename Y::block_type bblock;
91
92 // local solve at each block and immediate update
93 rowiterator rendi=A.beforeBegin();
94 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
95 {
96 bblock rhs(d[i.index()]);
97 coliterator j;
98 for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
99 (*j).mmv(v[j.index()],rhs);
100 algmeta_btsolve<I-1,diag,relax>::butsolve(*j,v[i.index()],rhs,w);
101 }
102 }
103 };
104
105 // recursion end ...
106 template<>
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)
110 {
111 A.solve(v,d);
112 v *= w;
113 }
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)
116 {
117 A.solve(v,d);
118 v *= w;
119 }
120 };
121 template<>
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& /*w*/)
125 {
126 A.solve(v,d);
127 }
128 template<class M, class X, class Y, class K>
129 static void butsolve (const M& A, X& v, const Y& d, const K& /*w*/)
130 {
131 A.solve(v,d);
132 }
133 };
134 template<>
135 struct algmeta_btsolve<0,nodiag,withrelax> {
136 template<class M, class X, class Y, class K>
137 static void bltsolve (const M& /*A*/, X& v, const Y& d, const K& w)
138 {
139 v = d;
140 v *= w;
141 }
142 template<class M, class X, class Y, class K>
143 static void butsolve (const M& /*A*/, X& v, const Y& d, const K& w)
144 {
145 v = d;
146 v *= w;
147 }
148 };
149 template<>
150 struct algmeta_btsolve<0,nodiag,norelax> {
151 template<class M, class X, class Y, class K>
152 static void bltsolve (const M& /*A*/, X& v, const Y& d, const K& /*w*/)
153 {
154 v = d;
155 }
156 template<class M, class X, class Y, class K>
157 static void butsolve (const M& /*A*/, X& v, const Y& d, const K& /*w*/)
158 {
159 v = d;
160 }
161 };
162
163
164 // user calls
165
166 // default block recursion level = 1
167
169 template<class M, class X, class Y>
170 void bltsolve (const M& A, X& v, const Y& d)
171 {
172 typename X::field_type w=1;
174 }
176 template<class M, class X, class Y, class K>
177 void bltsolve (const M& A, X& v, const Y& d, const K& w)
178 {
180 }
182 template<class M, class X, class Y>
183 void ubltsolve (const M& A, X& v, const Y& d)
184 {
185 typename X::field_type w=1;
187 }
189 template<class M, class X, class Y, class K>
190 void ubltsolve (const M& A, X& v, const Y& d, const K& w)
191 {
193 }
194
196 template<class M, class X, class Y>
197 void butsolve (const M& A, X& v, const Y& d)
198 {
199 typename X::field_type w=1;
201 }
203 template<class M, class X, class Y, class K>
204 void butsolve (const M& A, X& v, const Y& d, const K& w)
205 {
207 }
209 template<class M, class X, class Y>
210 void ubutsolve (const M& A, X& v, const Y& d)
211 {
212 typename X::field_type w=1;
214 }
216 template<class M, class X, class Y, class K>
217 void ubutsolve (const M& A, X& v, const Y& d, const K& w)
218 {
220 }
221
222 // general block recursion level >= 0
223
225 template<class M, class X, class Y, int l>
226 void bltsolve (const M& A, X& v, const Y& d, BL<l> /*bl*/)
227 {
228 typename X::field_type w=1;
230 }
232 template<class M, class X, class Y, class K, int l>
233 void bltsolve (const M& A, X& v, const Y& d, const K& w, BL<l> /*bl*/)
234 {
236 }
238 template<class M, class X, class Y, int l>
239 void ubltsolve (const M& A, X& v, const Y& d, BL<l> /*bl*/)
240 {
241 typename X::field_type w=1;
243 }
245 template<class M, class X, class Y, class K, int l>
246 void ubltsolve (const M& A, X& v, const Y& d, const K& w, BL<l> /*bl*/)
247 {
249 }
250
252 template<class M, class X, class Y, int l>
253 void butsolve (const M& A, X& v, const Y& d, BL<l> bl)
254 {
255 typename X::field_type w=1;
257 }
259 template<class M, class X, class Y, class K, int l>
260 void butsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
261 {
263 }
265 template<class M, class X, class Y, int l>
266 void ubutsolve (const M& A, X& v, const Y& d, BL<l> bl)
267 {
268 typename X::field_type w=1;
270 }
272 template<class M, class X, class Y, class K, int l>
273 void ubutsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
274 {
276 }
277
278
279
280 //============================================================
281 // generic block diagonal solves
282 // consider block decomposition A = L + D + U
283 // we can apply relaxation or not
284 // we can recurse over a fixed number of levels
285 //============================================================
286
287 // template meta program for diagonal solves
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)
292 {
293 // iterator types
294 typedef typename M::ConstRowIterator rowiterator;
295 typedef typename M::ConstColIterator coliterator;
296
297 // local solve at each block and immediate update
298 rowiterator rendi=A.beforeBegin();
299 for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
300 {
301 coliterator ii=(*i).find(i.index());
302 algmeta_bdsolve<I-1,relax>::bdsolve(*ii,v[i.index()],d[i.index()],w);
303 }
304 }
305 };
306
307 // recursion end ...
308 template<>
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)
312 {
313 A.solve(v,d);
314 v *= w;
315 }
316 };
317 template<>
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& /*w*/)
321 {
322 A.solve(v,d);
323 }
324 };
325
326 // user calls
327
328 // default block recursion level = 1
329
331 template<class M, class X, class Y>
332 void bdsolve (const M& A, X& v, const Y& d)
333 {
334 typename X::field_type w=1;
336 }
338 template<class M, class X, class Y, class K>
339 void bdsolve (const M& A, X& v, const Y& d, const K& w)
340 {
342 }
343
344 // general block recursion level >= 0
345
347 template<class M, class X, class Y, int l>
348 void bdsolve (const M& A, X& v, const Y& d, BL<l> /*bl*/)
349 {
350 typename X::field_type w=1;
352 }
354 template<class M, class X, class Y, class K, int l>
355 void bdsolve (const M& A, X& v, const Y& d, const K& w, BL<l> /*bl*/)
356 {
358 }
359
360
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
372 // template meta program for iterative solver steps
373 template<int I>
374 struct algmeta_itsteps {
375
376#if HAVE_DUNE_BOOST
377#ifdef HAVE_BOOST_FUSION
378
379 template<typename T11, typename T12, typename T13, typename T14,
380 typename T15, typename T16, typename T17, typename T18,
381 typename T19, typename T21, typename T22, typename T23,
382 typename T24, typename T25, typename T26, typename T27,
383 typename T28, typename T29, class K>
384 static void dbgs (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A,
385 MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
386 const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
387 const K& w)
388 {
389 const int rowcount =
390 boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
392 }
393#endif
394#endif
395
396 template<class M, class X, class Y, class K>
397 static void dbgs (const M& A, X& x, const Y& b, const K& w)
398 {
399 typedef typename M::ConstRowIterator rowiterator;
400 typedef typename M::ConstColIterator coliterator;
401 typedef typename Y::block_type bblock;
402 bblock rhs;
403
404 X xold(x); // remember old x
405
406 rowiterator endi=A.end();
407 for (rowiterator i=A.begin(); i!=endi; ++i)
408 {
409 rhs = b[i.index()]; // rhs = b_i
410 coliterator endj=(*i).end();
411 coliterator j=(*i).begin();
412 for (; j.index()<i.index(); ++j) // iterate over a_ij with j < i
413 (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
414 coliterator diag=j++; // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
415 for (; j != endj; ++j) // iterate over a_ij with j > i
416 (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j>i} a_ij * xold_j
417 algmeta_itsteps<I-1>::dbgs(*diag,x[i.index()],rhs,w); // if I==1: xnew_i = rhs/a_ii
418 }
419 // 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;
420 x *= w;
421 x.axpy(K(1)-w,xold);
422 }
423
424#if HAVE_DUNE_BOOST
425#ifdef HAVE_BOOST_FUSION
426
427 template<typename T11, typename T12, typename T13, typename T14,
428 typename T15, typename T16, typename T17, typename T18,
429 typename T19, typename T21, typename T22, typename T23,
430 typename T24, typename T25, typename T26, typename T27,
431 typename T28, typename T29, class K>
432 static void bsorf (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A,
433 MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
434 const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
435 const K& w)
436 {
437 const int rowcount =
438 boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
440 }
441#endif
442#endif
443
444 template<class M, class X, class Y, class K>
445 static void bsorf (const M& A, X& x, const Y& b, const K& w)
446 {
447 typedef typename M::ConstRowIterator rowiterator;
448 typedef typename M::ConstColIterator coliterator;
449 typedef typename Y::block_type bblock;
450 typedef typename X::block_type xblock;
451 bblock rhs;
452 xblock v;
453
454 // Initialize nested data structure if there are entries
455 if(A.begin()!=A.end())
456 v=x[0];
457
458 rowiterator endi=A.end();
459 for (rowiterator i=A.begin(); i!=endi; ++i)
460 {
461 rhs = b[i.index()]; // rhs = b_i
462 coliterator endj=(*i).end(); // iterate over a_ij with j < i
463 coliterator j=(*i).begin();
464 for (; j.index()<i.index(); ++j)
465 (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
466 coliterator diag=j; // *diag = a_ii
467 for (; j!=endj; ++j)
468 (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
469 algmeta_itsteps<I-1>::bsorf(*diag,v,rhs,w); // if blocksize I==1: v = rhs/a_ii
470 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)
471 }
472 }
473
474#if HAVE_DUNE_BOOST
475#ifdef HAVE_BOOST_FUSION
476
477 template<typename T11, typename T12, typename T13, typename T14,
478 typename T15, typename T16, typename T17, typename T18,
479 typename T19, typename T21, typename T22, typename T23,
480 typename T24, typename T25, typename T26, typename T27,
481 typename T28, typename T29, class K>
482 static void bsorb (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A,
483 MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
484 const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
485 const K& w)
486 {
487 const int rowcount =
488 mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
490 }
491#endif
492#endif
493
494 template<class M, class X, class Y, class K>
495 static void bsorb (const M& A, X& x, const Y& b, const K& w)
496 {
497 typedef typename M::ConstRowIterator rowiterator;
498 typedef typename M::ConstColIterator coliterator;
499 typedef typename Y::block_type bblock;
500 typedef typename X::block_type xblock;
501 bblock rhs;
502 xblock v;
503
504 // Initialize nested data structure if there are entries
505 if(A.begin()!=A.end())
506 v=x[0];
507
508 rowiterator endi=A.beforeBegin();
509 for (rowiterator i=A.beforeEnd(); i!=endi; --i)
510 {
511 rhs = b[i.index()];
512 coliterator endj=(*i).end();
513 coliterator j=(*i).begin();
514 for (; j.index()<i.index(); ++j)
515 (*j).mmv(x[j.index()],rhs);
516 coliterator diag=j;
517 for (; j!=endj; ++j)
518 (*j).mmv(x[j.index()],rhs);
519 algmeta_itsteps<I-1>::bsorb(*diag,v,rhs,w);
520 x[i.index()].axpy(w,v);
521 }
522 }
523
524#if HAVE_DUNE_BOOST
525#ifdef HAVE_BOOST_FUSION
526
527 template<typename T11, typename T12, typename T13, typename T14,
528 typename T15, typename T16, typename T17, typename T18,
529 typename T19, typename T21, typename T22, typename T23,
530 typename T24, typename T25, typename T26, typename T27,
531 typename T28, typename T29, class K>
532 static void dbjac (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A,
533 MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
534 const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
535 const K& w)
536 {
537 const int rowcount =
538 boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
540 }
541#endif
542#endif
543
544 template<class M, class X, class Y, class K>
545 static void dbjac (const M& A, X& x, const Y& b, const K& w)
546 {
547 typedef typename M::ConstRowIterator rowiterator;
548 typedef typename M::ConstColIterator coliterator;
549 typedef typename Y::block_type bblock;
550 bblock rhs;
551
552 X v(x); // allocate with same size
553
554 rowiterator endi=A.end();
555 for (rowiterator i=A.begin(); i!=endi; ++i)
556 {
557 rhs = b[i.index()];
558 coliterator endj=(*i).end();
559 coliterator j=(*i).begin();
560 for (; j.index()<i.index(); ++j)
561 (*j).mmv(x[j.index()],rhs);
562 coliterator diag=j;
563 for (; j!=endj; ++j)
564 (*j).mmv(x[j.index()],rhs);
565 algmeta_itsteps<I-1>::dbjac(*diag,v[i.index()],rhs,w);
566 }
567 x.axpy(w,v);
568 }
569 };
570 // end of recursion
571 template<>
572 struct algmeta_itsteps<0> {
573 template<class M, class X, class Y, class K>
574 static void dbgs (const M& A, X& x, const Y& b, const K& /*w*/)
575 {
576 A.solve(x,b);
577 }
578 template<class M, class X, class Y, class K>
579 static void bsorf (const M& A, X& x, const Y& b, const K& /*w*/)
580 {
581 A.solve(x,b);
582 }
583 template<class M, class X, class Y, class K>
584 static void bsorb (const M& A, X& x, const Y& b, const K& /*w*/)
585 {
586 A.solve(x,b);
587 }
588 template<class M, class X, class Y, class K>
589 static void dbjac (const M& A, X& x, const Y& b, const K& /*w*/)
590 {
591 A.solve(x,b);
592 }
593 };
594
595
596 // user calls
597
599 template<class M, class X, class Y, class K>
600 void dbgs (const M& A, X& x, const Y& b, const K& w)
601 {
603 }
605 template<class M, class X, class Y, class K, int l>
606 void dbgs (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
607 {
609 }
611 template<class M, class X, class Y, class K>
612 void bsorf (const M& A, X& x, const Y& b, const K& w)
613 {
615 }
617 template<class M, class X, class Y, class K, int l>
618 void bsorf (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
619 {
621 }
623 template<class M, class X, class Y, class K>
624 void bsorb (const M& A, X& x, const Y& b, const K& w)
625 {
627 }
629 template<class M, class X, class Y, class K, int l>
630 void bsorb (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
631 {
633 }
635 template<class M, class X, class Y, class K>
636 void dbjac (const M& A, X& x, const Y& b, const K& w)
637 {
639 }
641 template<class M, class X, class Y, class K, int l>
642 void dbjac (const M& A, X& x, const Y& b, const K& w, BL<l> /*bl*/)
643 {
645 }
646
647
650} // end namespace
651
652#endif
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:624
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:636
void dbgs(const M &A, X &x, const Y &b, const K &w)
GS step.
Definition: gsetc.hh:600
void bsorf(const M &A, X &x, const Y &b, const K &w, BL< l >)
SOR step.
Definition: gsetc.hh:618
void dbjac(const M &A, X &x, const Y &b, const K &w, BL< l >)
Jacobi step.
Definition: gsetc.hh:642
void dbgs(const M &A, X &x, const Y &b, const K &w, BL< l >)
GS step.
Definition: gsetc.hh:606
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:630
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:612
void ubutsolve(const M &A, X &v, const Y &d)
unit block upper triangular solve
Definition: gsetc.hh:210
Dune namespace.
Definition: alignment.hh:10
compile-time parameter for block recursion depth
Definition: gsetc.hh:40
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 21, 23:30, 2024)