Dune Core Modules (2.7.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 
12 #include <dune/common/hybridutilities.hh>
13 
14 #include "multitypeblockvector.hh"
15 #include "multitypeblockmatrix.hh"
16 
17 #include "istlexception.hh"
18 
19 
25 namespace 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);
529  algmeta_itsteps<I-1,typename M::block_type>::dbjac(*diag,v[i.index()],rhs,w);
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.80.0 (May 9, 22:29, 2024)