dune-istl  2.2.0
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>
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  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()];
411  coliterator endj=(*i).end();
412  coliterator j=(*i).begin();
413  for (; j.index()<i.index(); ++j)
414  (*j).mmv(x[j.index()],rhs);
415  coliterator diag=j;
416  for (; j != endj; ++j)
417  (*j).mmv(x[j.index()],rhs);
418  algmeta_itsteps<I-1>::dbgs(*diag,x[i.index()],rhs,w);
419  }
420  x *= w;
421  x.axpy(1-w,xold);
422  }
423 
424 #if HAVE_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()];
462  coliterator endj=(*i).end();
463  coliterator j=(*i).begin();
464  for (; j.index()<i.index(); ++j)
465  (*j).mmv(x[j.index()],rhs);
466  coliterator diag=j;
467  for (; j!=endj; ++j)
468  (*j).mmv(x[j.index()],rhs);
469  algmeta_itsteps<I-1>::bsorf(*diag,v,rhs,w);
470  x[i.index()].axpy(w,v);
471  }
472  }
473 
474 #if HAVE_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_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  typedef typename X::block_type xblock;
551  bblock rhs;
552 
553  X v(x); // allocate with same size
554 
555  rowiterator endi=A.end();
556  for (rowiterator i=A.begin(); i!=endi; ++i)
557  {
558  rhs = b[i.index()];
559  coliterator endj=(*i).end();
560  coliterator j=(*i).begin();
561  for (; j.index()<i.index(); ++j)
562  (*j).mmv(x[j.index()],rhs);
563  coliterator diag=j;
564  for (; j!=endj; ++j)
565  (*j).mmv(x[j.index()],rhs);
566  algmeta_itsteps<I-1>::dbjac(*diag,v[i.index()],rhs,w);
567  }
568  x.axpy(w,v);
569  }
570  };
571  // end of recursion
572  template<>
573  struct algmeta_itsteps<0> {
574  template<class M, class X, class Y, class K>
575  static void dbgs (const M& A, X& x, const Y& b, const K& w)
576  {
577  A.solve(x,b);
578  }
579  template<class M, class X, class Y, class K>
580  static void bsorf (const M& A, X& x, const Y& b, const K& w)
581  {
582  A.solve(x,b);
583  }
584  template<class M, class X, class Y, class K>
585  static void bsorb (const M& A, X& x, const Y& b, const K& w)
586  {
587  A.solve(x,b);
588  }
589  template<class M, class X, class Y, class K>
590  static void dbjac (const M& A, X& x, const Y& b, const K& w)
591  {
592  A.solve(x,b);
593  }
594  };
595 
596 
597  // user calls
598 
600  template<class M, class X, class Y, class K>
601  void dbgs (const M& A, X& x, const Y& b, const K& w)
602  {
603  algmeta_itsteps<1>::dbgs(A,x,b,w);
604  }
606  template<class M, class X, class Y, class K, int l>
607  void dbgs (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
608  {
609  algmeta_itsteps<l>::dbgs(A,x,b,w);
610  }
612  template<class M, class X, class Y, class K>
613  void bsorf (const M& A, X& x, const Y& b, const K& w)
614  {
615  algmeta_itsteps<1>::bsorf(A,x,b,w);
616  }
618  template<class M, class X, class Y, class K, int l>
619  void bsorf (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
620  {
621  algmeta_itsteps<l>::bsorf(A,x,b,w);
622  }
624  template<class M, class X, class Y, class K>
625  void bsorb (const M& A, X& x, const Y& b, const K& w)
626  {
627  algmeta_itsteps<1>::bsorb(A,x,b,w);
628  }
630  template<class M, class X, class Y, class K, int l>
631  void bsorb (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
632  {
633  algmeta_itsteps<l>::bsorb(A,x,b,w);
634  }
636  template<class M, class X, class Y, class K>
637  void dbjac (const M& A, X& x, const Y& b, const K& w)
638  {
639  algmeta_itsteps<1>::dbjac(A,x,b,w);
640  }
642  template<class M, class X, class Y, class K, int l>
643  void dbjac (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
644  {
645  algmeta_itsteps<l>::dbjac(A,x,b,w);
646  }
647 
648 
651 } // end namespace
652 
653 #endif