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 
22 namespace 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  {
602  algmeta_itsteps<1>::dbgs(A,x,b,w);
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  {
608  algmeta_itsteps<l>::dbgs(A,x,b,w);
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  {
614  algmeta_itsteps<1>::bsorf(A,x,b,w);
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  {
620  algmeta_itsteps<l>::bsorf(A,x,b,w);
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  {
626  algmeta_itsteps<1>::bsorb(A,x,b,w);
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  {
632  algmeta_itsteps<l>::bsorb(A,x,b,w);
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  {
638  algmeta_itsteps<1>::dbjac(A,x,b,w);
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  {
644  algmeta_itsteps<l>::dbjac(A,x,b,w);
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.80.0 (May 16, 22:29, 2024)