dune-istl  2.2.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_GSETC_HH
4 #define DUNE_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 {
47  };
48 
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<>
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<>
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<>
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<>
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>
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<>
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<>
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>
375 
376 #if HAVE_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>
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  typedef typename X::block_type xblock;
403  bblock rhs;
404 
405  X xold(x); // remember old x
406 
407  rowiterator endi=A.end();
408  for (rowiterator i=A.begin(); i!=endi; ++i)
409  {
410  rhs = b[i.index()]; // rhs = b_i
411  coliterator endj=(*i).end();
412  coliterator j=(*i).begin();
413  for (; j.index()<i.index(); ++j) // iterate over a_ij with j < i
414  (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
415  coliterator diag=j++; // *diag = a_ii and increment coliterator j from a_ii to a_i+1,i to skip diagonal
416  for (; j != endj; ++j) // iterate over a_ij with j > i
417  (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j>i} a_ij * xold_j
418  algmeta_itsteps<I-1>::dbgs(*diag,x[i.index()],rhs,w); // if I==1: xnew_i = rhs/a_ii
419  }
420  // 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;
421  x *= w;
422  x.axpy(K(1)-w,xold);
423  }
424 
425 #if HAVE_BOOST
426 #ifdef HAVE_BOOST_FUSION
427 
428  template<typename T11, typename T12, typename T13, typename T14,
429  typename T15, typename T16, typename T17, typename T18,
430  typename T19, typename T21, typename T22, typename T23,
431  typename T24, typename T25, typename T26, typename T27,
432  typename T28, typename T29, class K>
434  MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
435  const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
436  const K& w)
437  {
438  const int rowcount =
439  boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
441  }
442 #endif
443 #endif
444 
445  template<class M, class X, class Y, class K>
446  static void bsorf (const M& A, X& x, const Y& b, const K& w)
447  {
448  typedef typename M::ConstRowIterator rowiterator;
449  typedef typename M::ConstColIterator coliterator;
450  typedef typename Y::block_type bblock;
451  typedef typename X::block_type xblock;
452  bblock rhs;
453  xblock v;
454 
455  // Initialize nested data structure if there are entries
456  if(A.begin()!=A.end())
457  v=x[0];
458 
459  rowiterator endi=A.end();
460  for (rowiterator i=A.begin(); i!=endi; ++i)
461  {
462  rhs = b[i.index()]; // rhs = b_i
463  coliterator endj=(*i).end(); // iterate over a_ij with j < i
464  coliterator j=(*i).begin();
465  for (; j.index()<i.index(); ++j)
466  (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
467  coliterator diag=j; // *diag = a_ii
468  for (; j!=endj; ++j)
469  (*j).mmv(x[j.index()],rhs); // rhs -= sum_{j<i} a_ij * xnew_j
470  algmeta_itsteps<I-1>::bsorf(*diag,v,rhs,w); // if blocksize I==1: v = rhs/a_ii
471  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)
472  }
473  }
474 
475 #if HAVE_BOOST
476 #ifdef HAVE_BOOST_FUSION
477 
478  template<typename T11, typename T12, typename T13, typename T14,
479  typename T15, typename T16, typename T17, typename T18,
480  typename T19, typename T21, typename T22, typename T23,
481  typename T24, typename T25, typename T26, typename T27,
482  typename T28, typename T29, class K>
484  MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
485  const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
486  const K& w)
487  {
488  const int rowcount =
489  mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
491  }
492 #endif
493 #endif
494 
495  template<class M, class X, class Y, class K>
496  static void bsorb (const M& A, X& x, const Y& b, const K& w)
497  {
498  typedef typename M::ConstRowIterator rowiterator;
499  typedef typename M::ConstColIterator coliterator;
500  typedef typename Y::block_type bblock;
501  typedef typename X::block_type xblock;
502  bblock rhs;
503  xblock v;
504 
505  // Initialize nested data structure if there are entries
506  if(A.begin()!=A.end())
507  v=x[0];
508 
509  rowiterator endi=A.beforeBegin();
510  for (rowiterator i=A.beforeEnd(); i!=endi; --i)
511  {
512  rhs = b[i.index()];
513  coliterator endj=(*i).end();
514  coliterator j=(*i).begin();
515  for (; j.index()<i.index(); ++j)
516  (*j).mmv(x[j.index()],rhs);
517  coliterator diag=j;
518  for (; j!=endj; ++j)
519  (*j).mmv(x[j.index()],rhs);
520  algmeta_itsteps<I-1>::bsorb(*diag,v,rhs,w);
521  x[i.index()].axpy(w,v);
522  }
523  }
524 
525 #if HAVE_BOOST
526 #ifdef HAVE_BOOST_FUSION
527 
528  template<typename T11, typename T12, typename T13, typename T14,
529  typename T15, typename T16, typename T17, typename T18,
530  typename T19, typename T21, typename T22, typename T23,
531  typename T24, typename T25, typename T26, typename T27,
532  typename T28, typename T29, class K>
534  MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x,
535  const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b,
536  const K& w)
537  {
538  const int rowcount =
539  boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
541  }
542 #endif
543 #endif
544 
545  template<class M, class X, class Y, class K>
546  static void dbjac (const M& A, X& x, const Y& b, const K& w)
547  {
548  typedef typename M::ConstRowIterator rowiterator;
549  typedef typename M::ConstColIterator coliterator;
550  typedef typename Y::block_type bblock;
551  typedef typename X::block_type xblock;
552  bblock rhs;
553 
554  X v(x); // allocate with same size
555 
556  rowiterator endi=A.end();
557  for (rowiterator i=A.begin(); i!=endi; ++i)
558  {
559  rhs = b[i.index()];
560  coliterator endj=(*i).end();
561  coliterator j=(*i).begin();
562  for (; j.index()<i.index(); ++j)
563  (*j).mmv(x[j.index()],rhs);
564  coliterator diag=j;
565  for (; j!=endj; ++j)
566  (*j).mmv(x[j.index()],rhs);
567  algmeta_itsteps<I-1>::dbjac(*diag,v[i.index()],rhs,w);
568  }
569  x.axpy(w,v);
570  }
571  };
572  // end of recursion
573  template<>
574  struct algmeta_itsteps<0> {
575  template<class M, class X, class Y, class K>
576  static void dbgs (const M& A, X& x, const Y& b, const K& w)
577  {
578  A.solve(x,b);
579  }
580  template<class M, class X, class Y, class K>
581  static void bsorf (const M& A, X& x, const Y& b, const K& w)
582  {
583  A.solve(x,b);
584  }
585  template<class M, class X, class Y, class K>
586  static void bsorb (const M& A, X& x, const Y& b, const K& w)
587  {
588  A.solve(x,b);
589  }
590  template<class M, class X, class Y, class K>
591  static void dbjac (const M& A, X& x, const Y& b, const K& w)
592  {
593  A.solve(x,b);
594  }
595  };
596 
597 
598  // user calls
599 
601  template<class M, class X, class Y, class K>
602  void dbgs (const M& A, X& x, const Y& b, const K& w)
603  {
604  algmeta_itsteps<1>::dbgs(A,x,b,w);
605  }
607  template<class M, class X, class Y, class K, int l>
608  void dbgs (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
609  {
610  algmeta_itsteps<l>::dbgs(A,x,b,w);
611  }
613  template<class M, class X, class Y, class K>
614  void bsorf (const M& A, X& x, const Y& b, const K& w)
615  {
616  algmeta_itsteps<1>::bsorf(A,x,b,w);
617  }
619  template<class M, class X, class Y, class K, int l>
620  void bsorf (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
621  {
622  algmeta_itsteps<l>::bsorf(A,x,b,w);
623  }
625  template<class M, class X, class Y, class K>
626  void bsorb (const M& A, X& x, const Y& b, const K& w)
627  {
628  algmeta_itsteps<1>::bsorb(A,x,b,w);
629  }
631  template<class M, class X, class Y, class K, int l>
632  void bsorb (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
633  {
634  algmeta_itsteps<l>::bsorb(A,x,b,w);
635  }
637  template<class M, class X, class Y, class K>
638  void dbjac (const M& A, X& x, const Y& b, const K& w)
639  {
640  algmeta_itsteps<1>::dbjac(A,x,b,w);
641  }
643  template<class M, class X, class Y, class K, int l>
644  void dbjac (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
645  {
646  algmeta_itsteps<l>::dbjac(A,x,b,w);
647  }
648 
649 
652 } // end namespace
653 
654 #endif