Dune Core Modules (unstable)

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 
27 namespace 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);
535  algmeta_itsteps<I-1,typename M::block_type>::dbjac(*diag,v[i.index()],rhs,w);
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:569
static void dbgs(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:485
static void bsorb(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:541
static constexpr size_type N()
Return the number of matrix rows.
Definition: multitypeblockmatrix.hh:64
static void bsorf(const TMatrix &A, TVector &x, const TVector &b, const K &w)
Definition: multitypeblockmatrix.hh:514
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.80.0 (May 1, 22:29, 2024)