Dune Core Modules (2.5.0)

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 // generic steps of iteration methods
363 // Jacobi, Gauss-Seidel, SOR, SSOR
364 // work directly on Ax=b, ie solve M(x^{i+1}-x^i) = w (b-Ax^i)
365 // we can recurse over a fixed number of levels
366 //============================================================
367
368 // template meta program for iterative solver steps
369 template<int I, typename M>
370 struct algmeta_itsteps {
371
372 template<class X, class Y, class K>
373 static void dbgs (const M& A, X& x, const Y& b, const K& w)
374 {
375 typedef typename M::ConstRowIterator rowiterator;
376 typedef typename M::ConstColIterator coliterator;
377 typedef typename Y::block_type bblock;
378 bblock rhs;
379
380 X xold(x); // remember old x
381
382 rowiterator endi=A.end();
383 for (rowiterator i=A.begin(); i!=endi; ++i)
384 {
385 rhs = b[i.index()]; // rhs = b_i
386 coliterator endj=(*i).end();
387 coliterator j=(*i).begin();
388 for (; j.index()<i.index(); ++j) // iterate over a_ij with j < i
389 (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
390 coliterator diag=j++; // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
391 for (; j != endj; ++j) // iterate over a_ij with j > i
392 (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j>i} a_ij * xold_j
393 algmeta_itsteps<I-1,typename M::block_type>::dbgs(*diag,x[i.index()],rhs,w); // if I==1: xnew_i = rhs/a_ii
394 }
395 // 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;
396 x *= w;
397 x.axpy(K(1)-w,xold);
398 }
399
400 template<class X, class Y, class K>
401 static void bsorf (const M& A, X& x, const Y& b, const K& w)
402 {
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;
407 bblock rhs;
408 xblock v;
409
410 // Initialize nested data structure if there are entries
411 if(A.begin()!=A.end())
412 v=x[0];
413
414 rowiterator endi=A.end();
415 for (rowiterator i=A.begin(); i!=endi; ++i)
416 {
417 rhs = b[i.index()]; // rhs = b_i
418 coliterator endj=(*i).end(); // iterate over a_ij with j < i
419 coliterator j=(*i).begin();
420 for (; j.index()<i.index(); ++j)
421 (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
422 coliterator diag=j; // *diag = a_ii
423 for (; j!=endj; ++j)
424 (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
425 algmeta_itsteps<I-1,typename M::block_type>::bsorf(*diag,v,rhs,w); // if blocksize I==1: v = rhs/a_ii
426 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)
427 }
428 }
429
430 template<class X, class Y, class K>
431 static void bsorb (const M& A, X& x, const Y& b, const K& w)
432 {
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;
437 bblock rhs;
438 xblock v;
439
440 // Initialize nested data structure if there are entries
441 if(A.begin()!=A.end())
442 v=x[0];
443
444 rowiterator endi=A.beforeBegin();
445 for (rowiterator i=A.beforeEnd(); i!=endi; --i)
446 {
447 rhs = b[i.index()];
448 coliterator endj=(*i).end();
449 coliterator j=(*i).begin();
450 for (; j.index()<i.index(); ++j)
451 (*j).mmv(x[j.index()],rhs);
452 coliterator diag=j;
453 for (; j!=endj; ++j)
454 (*j).mmv(x[j.index()],rhs);
456 x[i.index()].axpy(w,v);
457 }
458 }
459
460 template<class X, class Y, class K>
461 static void dbjac (const M& A, X& x, const Y& b, const K& w)
462 {
463 typedef typename M::ConstRowIterator rowiterator;
464 typedef typename M::ConstColIterator coliterator;
465 typedef typename Y::block_type bblock;
466 bblock rhs;
467
468 X v(x); // allocate with same size
469
470 rowiterator endi=A.end();
471 for (rowiterator i=A.begin(); i!=endi; ++i)
472 {
473 rhs = b[i.index()];
474 coliterator endj=(*i).end();
475 coliterator j=(*i).begin();
476 for (; j.index()<i.index(); ++j)
477 (*j).mmv(x[j.index()],rhs);
478 coliterator diag=j;
479 for (; j!=endj; ++j)
480 (*j).mmv(x[j.index()],rhs);
482 }
483 x.axpy(w,v);
484 }
485 };
486 // end of recursion
487 template<typename M>
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& /*w*/)
491 {
492 A.solve(x,b);
493 }
494 template<class X, class Y, class K>
495 static void bsorf (const M& A, X& x, const Y& b, const K& /*w*/)
496 {
497 A.solve(x,b);
498 }
499 template<class X, class Y, class K>
500 static void bsorb (const M& A, X& x, const Y& b, const K& /*w*/)
501 {
502 A.solve(x,b);
503 }
504 template<class X, class Y, class K>
505 static void dbjac (const M& A, X& x, const Y& b, const K& /*w*/)
506 {
507 A.solve(x,b);
508 }
509 };
510
511 template<int I, typename T1, typename... MultiTypeMatrixArgs>
512 struct algmeta_itsteps<I,MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>> {
513 template<
514 typename... MultiTypeVectorArgs,
515 class K>
516 static void dbgs (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
517 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
518 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
519 const K& w)
520 {
523 }
524
525 template<
526 typename... MultiTypeVectorArgs,
527 class K>
528 static void bsorf (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
529 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
530 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
531 const K& w)
532 {
535 }
536
537 template<
538 typename... MultiTypeVectorArgs,
539 class K>
540 static void bsorb (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
541 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
542 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
543 const K& w)
544 {
547 }
548
549 template<
550 typename... MultiTypeVectorArgs,
551 class K
552 >
553 static void dbjac (const MultiTypeBlockMatrix<T1, MultiTypeMatrixArgs...>& A,
554 MultiTypeBlockVector<MultiTypeVectorArgs...>& x,
555 const MultiTypeBlockVector<MultiTypeVectorArgs...>& b,
556 const K& w)
557 {
560 }
561 };
562
563 // user calls
564
566 template<class M, class X, class Y, class K>
567 void dbgs (const M& A, X& x, const Y& b, const K& w)
568 {
570 }
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> /*bl*/)
574 {
576 }
578 template<class M, class X, class Y, class K>
579 void bsorf (const M& A, X& x, const Y& b, const K& w)
580 {
582 }
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> /*bl*/)
586 {
588 }
590 template<class M, class X, class Y, class K>
591 void bsorb (const M& A, X& x, const Y& b, const K& w)
592 {
594 }
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> /*bl*/)
598 {
599 algmeta_itsteps<l,typename std::remove_cv<M>::type>::bsorb(A,x,b,w);
600 }
602 template<class M, class X, class Y, class K>
603 void dbjac (const M& A, X& x, const Y& b, const K& w)
604 {
606 }
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> /*bl*/)
610 {
612 }
613
614
617} // end namespace
618
619#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: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
Creative Commons License   |  Legal Statements / Impressum  |  Hosted by TU Dresden  |  generated with Hugo v0.111.3 (Nov 12, 23:30, 2024)