DUNE PDELab (git)

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