dune-istl  2.2.0
gsetc.hh
Go to the documentation of this file.
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*-
00002 // vi: set et ts=4 sw=2 sts=2:
00003 #ifndef DUNE_GSETC_HH
00004 #define DUNE_GSETC_HH
00005 
00006 #include<cmath>
00007 #include<complex>
00008 #include<iostream>
00009 #include<iomanip>
00010 #include<string>
00011 #include "multitypeblockvector.hh"
00012 #include "multitypeblockmatrix.hh"
00013 
00014 #include "istlexception.hh"
00015 
00016 
00022 namespace Dune {
00023    
00034   //============================================================
00035   // parameter types
00036   //============================================================
00037 
00039   template<int l>
00040   struct BL {
00041         enum {recursion_level = l};
00042   };
00043 
00044   enum WithDiagType {
00045         withdiag=1,
00046         nodiag=0
00047   };
00048 
00049   enum WithRelaxType {
00050         withrelax=1,
00051         norelax=0
00052   };
00053 
00054   //============================================================
00055   // generic triangular solves
00056   // consider block decomposition A = L + D + U
00057   // we can invert L, L+D, U, U+D
00058   // we can apply relaxation or not
00059   // we can recurse over a fixed number of levels
00060   //============================================================
00061 
00062   // template meta program for triangular solves
00063   template<int I, WithDiagType diag, WithRelaxType relax>
00064   struct algmeta_btsolve {
00065         template<class M, class X, class Y, class K>
00066         static void bltsolve (const M& A, X& v, const Y& d, const K& w)
00067         {
00068           // iterator types
00069           typedef typename M::ConstRowIterator rowiterator;
00070           typedef typename M::ConstColIterator coliterator;
00071           typedef typename Y::block_type bblock;
00072 
00073           // local solve at each block and immediate update
00074           rowiterator endi=A.end();
00075           for (rowiterator i=A.begin(); i!=endi; ++i)
00076                 {
00077                   bblock rhs(d[i.index()]);
00078                   coliterator j;
00079                   for (j=(*i).begin(); j.index()<i.index(); ++j)
00080                         (*j).mmv(v[j.index()],rhs);
00081                   algmeta_btsolve<I-1,diag,relax>::bltsolve(*j,v[i.index()],rhs,w);
00082                 }
00083         }
00084         template<class M, class X, class Y, class K>
00085         static void butsolve (const M& A, X& v, const Y& d, const K& w)
00086         {
00087           // iterator types
00088           typedef typename M::ConstRowIterator rowiterator;
00089           typedef typename M::ConstColIterator coliterator;
00090           typedef typename Y::block_type bblock;
00091 
00092           // local solve at each block and immediate update
00093           rowiterator rendi=A.beforeBegin();
00094           for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
00095                 {
00096                   bblock rhs(d[i.index()]);
00097                   coliterator j;
00098                   for (j=(*i).beforeEnd(); j.index()>i.index(); --j)
00099                         (*j).mmv(v[j.index()],rhs);
00100                   algmeta_btsolve<I-1,diag,relax>::butsolve(*j,v[i.index()],rhs,w);
00101                 }
00102         }
00103   };
00104 
00105   // recursion end ...
00106   template<>
00107   struct algmeta_btsolve<0,withdiag,withrelax> {
00108         template<class M, class X, class Y, class K>
00109         static void bltsolve (const M& A, X& v, const Y& d, const K& w)
00110         {
00111           A.solve(v,d);
00112           v *= w;
00113         }
00114         template<class M, class X, class Y, class K>
00115         static void butsolve (const M& A, X& v, const Y& d, const K& w)
00116         {
00117           A.solve(v,d);
00118           v *= w;
00119         }
00120   };
00121   template<>
00122   struct algmeta_btsolve<0,withdiag,norelax> {
00123         template<class M, class X, class Y, class K>
00124         static void bltsolve (const M& A, X& v, const Y& d, const K& w)
00125         {
00126           A.solve(v,d);
00127         }
00128         template<class M, class X, class Y, class K>
00129         static void butsolve (const M& A, X& v, const Y& d, const K& w)
00130         {
00131           A.solve(v,d);
00132         }
00133   };
00134   template<>
00135   struct algmeta_btsolve<0,nodiag,withrelax> {
00136         template<class M, class X, class Y, class K>
00137         static void bltsolve (const M& A, X& v, const Y& d, const K& w)
00138         {
00139           v = d;
00140           v *= w;
00141         }
00142         template<class M, class X, class Y, class K>
00143         static void butsolve (const M& A, X& v, const Y& d, const K& w)
00144         {
00145           v = d;
00146           v *= w;
00147         }
00148   };
00149   template<>
00150   struct algmeta_btsolve<0,nodiag,norelax> {
00151         template<class M, class X, class Y, class K>
00152         static void bltsolve (const M& A, X& v, const Y& d, const K& w)
00153         {
00154           v = d;
00155         }
00156         template<class M, class X, class Y, class K>
00157         static void butsolve (const M& A, X& v, const Y& d, const K& w)
00158         {
00159           v = d;
00160         }
00161   };
00162 
00163 
00164   // user calls
00165 
00166   // default block recursion level = 1
00167 
00169   template<class M, class X, class Y>
00170   void bltsolve (const M& A, X& v, const Y& d)
00171   {
00172         typename X::field_type w=1;
00173         algmeta_btsolve<1,withdiag,norelax>::bltsolve(A,v,d,w);
00174   }
00176   template<class M, class X, class Y, class K>
00177   void bltsolve (const M& A, X& v, const Y& d, const K& w)
00178   {
00179         algmeta_btsolve<1,withdiag,withrelax>::bltsolve(A,v,d,w);
00180   }
00182   template<class M, class X, class Y>
00183   void ubltsolve (const M& A, X& v, const Y& d)
00184   {
00185         typename X::field_type w=1;
00186         algmeta_btsolve<1,nodiag,norelax>::bltsolve(A,v,d,w);
00187   }
00189   template<class M, class X, class Y, class K>
00190   void ubltsolve (const M& A, X& v, const Y& d, const K& w)
00191   {
00192         algmeta_btsolve<1,nodiag,withrelax>::bltsolve(A,v,d,w);
00193   }
00194 
00196   template<class M, class X, class Y>
00197   void butsolve (const M& A, X& v, const Y& d)
00198   {
00199         typename X::field_type w=1;
00200         algmeta_btsolve<1,withdiag,norelax>::butsolve(A,v,d,w);
00201   }
00203   template<class M, class X, class Y, class K>
00204   void butsolve (const M& A, X& v, const Y& d, const K& w)
00205   {
00206         algmeta_btsolve<1,withdiag,withrelax>::butsolve(A,v,d,w);
00207   }
00209   template<class M, class X, class Y>
00210   void ubutsolve (const M& A, X& v, const Y& d)
00211   {
00212         typename X::field_type w=1;
00213         algmeta_btsolve<1,nodiag,norelax>::butsolve(A,v,d,w);
00214   }
00216   template<class M, class X, class Y, class K>
00217   void ubutsolve (const M& A, X& v, const Y& d, const K& w)
00218   {
00219         algmeta_btsolve<1,nodiag,withrelax>::butsolve(A,v,d,w);
00220   }
00221 
00222   // general block recursion level >= 0
00223 
00225   template<class M, class X, class Y, int l>
00226   void bltsolve (const M& A, X& v, const Y& d, BL<l> bl)
00227   {
00228         typename X::field_type w=1;
00229         algmeta_btsolve<l,withdiag,norelax>::bltsolve(A,v,d,w);
00230   }
00232   template<class M, class X, class Y, class K, int l>
00233   void bltsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
00234   {
00235         algmeta_btsolve<l,withdiag,withrelax>::bltsolve(A,v,d,w);
00236   }
00238   template<class M, class X, class Y, int l>
00239   void ubltsolve (const M& A, X& v, const Y& d, BL<l> bl)
00240   {
00241         typename X::field_type w=1;
00242         algmeta_btsolve<l,nodiag,norelax>::bltsolve(A,v,d,w);
00243   }
00245   template<class M, class X, class Y, class K, int l>
00246   void ubltsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
00247   {
00248         algmeta_btsolve<l,nodiag,withrelax>::bltsolve(A,v,d,w);
00249   }
00250 
00252   template<class M, class X, class Y, int l>
00253   void butsolve (const M& A, X& v, const Y& d, BL<l> bl)
00254   {
00255         typename X::field_type w=1;
00256         algmeta_btsolve<l,withdiag,norelax>::butsolve(A,v,d,w);
00257   }
00259   template<class M, class X, class Y, class K, int l>
00260   void butsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
00261   {
00262         algmeta_btsolve<l,withdiag,withrelax>::butsolve(A,v,d,w);
00263   }
00265   template<class M, class X, class Y, int l>
00266   void ubutsolve (const M& A, X& v, const Y& d, BL<l> bl)
00267   {
00268         typename X::field_type w=1;
00269         algmeta_btsolve<l,nodiag,norelax>::butsolve(A,v,d,w);
00270   }
00272   template<class M, class X, class Y, class K, int l>
00273   void ubutsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
00274   {
00275         algmeta_btsolve<l,nodiag,withrelax>::butsolve(A,v,d,w);
00276   }
00277 
00278 
00279 
00280   //============================================================
00281   // generic block diagonal solves
00282   // consider block decomposition A = L + D + U
00283   // we can apply relaxation or not
00284   // we can recurse over a fixed number of levels
00285   //============================================================
00286 
00287   // template meta program for diagonal solves
00288   template<int I, WithRelaxType relax>
00289   struct algmeta_bdsolve {
00290         template<class M, class X, class Y, class K>
00291         static void bdsolve (const M& A, X& v, const Y& d, const K& w)
00292         {
00293           // iterator types
00294           typedef typename M::ConstRowIterator rowiterator;
00295           typedef typename M::ConstColIterator coliterator;
00296 
00297           // local solve at each block and immediate update
00298           rowiterator rendi=A.beforeBegin();
00299           for (rowiterator i=A.beforeEnd(); i!=rendi; --i)
00300                 {
00301                   coliterator ii=(*i).find(i.index());
00302                   algmeta_bdsolve<I-1,relax>::bdsolve(*ii,v[i.index()],d[i.index()],w);
00303                 }
00304         }
00305   };
00306 
00307   // recursion end ...
00308   template<>
00309   struct algmeta_bdsolve<0,withrelax> {
00310         template<class M, class X, class Y, class K>
00311         static void bdsolve (const M& A, X& v, const Y& d, const K& w)
00312         {
00313           A.solve(v,d);
00314           v *= w;
00315         }
00316   };
00317   template<>
00318   struct algmeta_bdsolve<0,norelax> {
00319         template<class M, class X, class Y, class K>
00320         static void bdsolve (const M& A, X& v, const Y& d, const K& w)
00321         {
00322           A.solve(v,d);
00323         }
00324   };
00325 
00326   // user calls
00327 
00328   // default block recursion level = 1
00329 
00331   template<class M, class X, class Y>
00332   void bdsolve (const M& A, X& v, const Y& d)
00333   {
00334         typename X::field_type w=1;
00335         algmeta_bdsolve<1,norelax>::bdsolve(A,v,d,w);
00336   }
00338   template<class M, class X, class Y, class K>
00339   void bdsolve (const M& A, X& v, const Y& d, const K& w)
00340   {
00341         algmeta_bdsolve<1,withrelax>::bdsolve(A,v,d,w);
00342   }
00343 
00344   // general block recursion level >= 0
00345 
00347   template<class M, class X, class Y, int l>
00348   void bdsolve (const M& A, X& v, const Y& d, BL<l> bl)
00349   {
00350         typename X::field_type w=1;
00351         algmeta_bdsolve<l,norelax>::bdsolve(A,v,d,w);
00352   }
00354   template<class M, class X, class Y, class K, int l>
00355   void bdsolve (const M& A, X& v, const Y& d, const K& w, BL<l> bl)
00356   {
00357         algmeta_bdsolve<l,withrelax>::bdsolve(A,v,d,w);
00358   }
00359 
00360 
00361 
00362 
00363 
00364   //============================================================
00365   // generic steps of iteration methods
00366   // Jacobi, Gauss-Seidel, SOR, SSOR
00367   // work directly on Ax=b, ie solve M(x^{i+1}-x^i) = w (b-Ax^i)
00368   // we can recurse over a fixed number of levels
00369   //============================================================
00370 
00371 
00372   // template meta program for iterative solver steps
00373   template<int I>
00374   struct algmeta_itsteps {
00375 
00376 #if HAVE_BOOST
00377 #ifdef HAVE_BOOST_FUSION
00378 
00379     template<typename T11, typename T12, typename T13, typename T14,
00380              typename T15, typename T16, typename T17, typename T18, 
00381              typename T19, typename T21, typename T22, typename T23, 
00382              typename T24, typename T25, typename T26, typename T27, 
00383              typename T28, typename T29, class K>
00384     static void dbgs (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A, 
00385                       MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x, 
00386                       const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b, 
00387                       const K& w)
00388     {   
00389       const int rowcount = 
00390         boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
00391       Dune::MultiTypeBlockMatrix_Solver<I,0,rowcount>::dbgs(A, x, b, w);
00392     }
00393 #endif
00394 #endif
00395 
00396         template<class M, class X, class Y, class K>
00397         static void dbgs (const M& A, X& x, const Y& b, const K& w)
00398         {
00399           typedef typename M::ConstRowIterator rowiterator;
00400           typedef typename M::ConstColIterator coliterator;
00401           typedef typename Y::block_type bblock;
00402           typedef typename X::block_type xblock;
00403           bblock rhs;
00404 
00405           X xold(x); // remember old x
00406 
00407           rowiterator endi=A.end();
00408           for (rowiterator i=A.begin(); i!=endi; ++i)
00409                 {
00410                   rhs = b[i.index()];
00411                   coliterator endj=(*i).end();
00412                   coliterator j=(*i).begin();
00413                   for (; j.index()<i.index(); ++j)
00414                         (*j).mmv(x[j.index()],rhs);
00415                   coliterator diag=j; 
00416                   for (; j != endj; ++j)
00417                         (*j).mmv(x[j.index()],rhs);
00418                   algmeta_itsteps<I-1>::dbgs(*diag,x[i.index()],rhs,w);
00419                 }
00420           x *= w;
00421           x.axpy(1-w,xold);
00422         }
00423 
00424 #if HAVE_BOOST
00425 #ifdef HAVE_BOOST_FUSION
00426 
00427     template<typename T11, typename T12, typename T13, typename T14, 
00428              typename T15, typename T16, typename T17, typename T18, 
00429              typename T19, typename T21, typename T22, typename T23, 
00430              typename T24, typename T25, typename T26, typename T27, 
00431              typename T28, typename T29, class K>
00432     static void bsorf (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A, 
00433                        MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x, 
00434                        const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b, 
00435                        const K& w)
00436     {   
00437       const int rowcount = 
00438         boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
00439       Dune::MultiTypeBlockMatrix_Solver<I,0,rowcount>::bsorf(A, x, b, w);
00440     }
00441 #endif
00442 #endif
00443 
00444         template<class M, class X, class Y, class K>
00445         static void bsorf (const M& A, X& x, const Y& b, const K& w)
00446         {
00447           typedef typename M::ConstRowIterator rowiterator;
00448           typedef typename M::ConstColIterator coliterator;
00449           typedef typename Y::block_type bblock;
00450           typedef typename X::block_type xblock;
00451           bblock rhs;
00452           xblock v;
00453 
00454           // Initialize nested data structure if there are entries
00455           if(A.begin()!=A.end())
00456             v=x[0];
00457 
00458           rowiterator endi=A.end();
00459           for (rowiterator i=A.begin(); i!=endi; ++i)
00460                 {
00461                   rhs = b[i.index()];
00462                   coliterator endj=(*i).end();
00463                   coliterator j=(*i).begin();
00464                   for (; j.index()<i.index(); ++j)
00465                         (*j).mmv(x[j.index()],rhs);
00466                   coliterator diag=j; 
00467                   for (; j!=endj; ++j)
00468                         (*j).mmv(x[j.index()],rhs);
00469                   algmeta_itsteps<I-1>::bsorf(*diag,v,rhs,w);
00470                   x[i.index()].axpy(w,v);
00471                 }
00472         }
00473 
00474 #if HAVE_BOOST
00475 #ifdef HAVE_BOOST_FUSION
00476 
00477     template<typename T11, typename T12, typename T13, typename T14, 
00478              typename T15, typename T16, typename T17, typename T18, 
00479              typename T19, typename T21, typename T22, typename T23, 
00480              typename T24, typename T25, typename T26, typename T27, 
00481              typename T28, typename T29, class K>
00482     static void bsorb (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A, 
00483                        MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x, 
00484                        const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b, 
00485                        const K& w)
00486     {   
00487       const int rowcount = 
00488         mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
00489       Dune::MultiTypeBlockMatrix_Solver<I,rowcount-1,rowcount>::bsorb(A, x, b, w);
00490     }
00491 #endif
00492 #endif
00493 
00494         template<class M, class X, class Y, class K>
00495         static void bsorb (const M& A, X& x, const Y& b, const K& w)
00496         {
00497           typedef typename M::ConstRowIterator rowiterator;
00498           typedef typename M::ConstColIterator coliterator;
00499           typedef typename Y::block_type bblock;
00500           typedef typename X::block_type xblock;
00501           bblock rhs;
00502           xblock v;
00503 
00504           // Initialize nested data structure if there are entries
00505           if(A.begin()!=A.end())
00506             v=x[0];
00507           
00508           rowiterator endi=A.beforeBegin();
00509           for (rowiterator i=A.beforeEnd(); i!=endi; --i)
00510                 {
00511                   rhs = b[i.index()];
00512                   coliterator endj=(*i).end();
00513                   coliterator j=(*i).begin();
00514                   for (; j.index()<i.index(); ++j)
00515                         (*j).mmv(x[j.index()],rhs);
00516                   coliterator diag=j; 
00517                   for (; j!=endj; ++j)
00518                         (*j).mmv(x[j.index()],rhs);
00519                   algmeta_itsteps<I-1>::bsorb(*diag,v,rhs,w);
00520                   x[i.index()].axpy(w,v);
00521                 }
00522         }
00523 
00524 #if HAVE_BOOST
00525 #ifdef HAVE_BOOST_FUSION
00526 
00527     template<typename T11, typename T12, typename T13, typename T14, 
00528              typename T15, typename T16, typename T17, typename T18, 
00529              typename T19, typename T21, typename T22, typename T23, 
00530              typename T24, typename T25, typename T26, typename T27, 
00531              typename T28, typename T29, class K>
00532     static void dbjac (const MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19>& A,
00533                        MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& x, 
00534                        const MultiTypeBlockVector<T21,T22,T23,T24,T25,T26,T27,T28,T29>& b, 
00535                        const K& w)
00536     {   
00537       const int rowcount = 
00538         boost::mpl::size<MultiTypeBlockMatrix<T11,T12,T13,T14,T15,T16,T17,T18,T19> >::value;
00539       Dune::MultiTypeBlockMatrix_Solver<I,0,rowcount >::dbjac(A, x, b, w);
00540     }
00541 #endif
00542 #endif
00543 
00544         template<class M, class X, class Y, class K>
00545         static void dbjac (const M& A, X& x, const Y& b, const K& w)
00546         {
00547           typedef typename M::ConstRowIterator rowiterator;
00548           typedef typename M::ConstColIterator coliterator;
00549           typedef typename Y::block_type bblock;
00550           typedef typename X::block_type xblock;
00551           bblock rhs;
00552 
00553           X v(x); // allocate with same size
00554 
00555           rowiterator endi=A.end();
00556           for (rowiterator i=A.begin(); i!=endi; ++i)
00557                 {
00558                   rhs = b[i.index()];
00559                   coliterator endj=(*i).end();
00560                   coliterator j=(*i).begin();
00561                   for (; j.index()<i.index(); ++j)
00562                         (*j).mmv(x[j.index()],rhs);
00563                   coliterator diag=j; 
00564                   for (; j!=endj; ++j)
00565                         (*j).mmv(x[j.index()],rhs);
00566                   algmeta_itsteps<I-1>::dbjac(*diag,v[i.index()],rhs,w);
00567                 }
00568           x.axpy(w,v);
00569         }
00570   };
00571   // end of recursion
00572   template<>
00573   struct algmeta_itsteps<0> {
00574         template<class M, class X, class Y, class K>
00575         static void dbgs (const M& A, X& x, const Y& b, const K& w)
00576         {
00577           A.solve(x,b);
00578         }
00579         template<class M, class X, class Y, class K>
00580         static void bsorf (const M& A, X& x, const Y& b, const K& w)
00581         {
00582           A.solve(x,b);
00583         }
00584         template<class M, class X, class Y, class K>
00585         static void bsorb (const M& A, X& x, const Y& b, const K& w)
00586         {
00587           A.solve(x,b);
00588         }
00589         template<class M, class X, class Y, class K>
00590         static void dbjac (const M& A, X& x, const Y& b, const K& w)
00591         {
00592           A.solve(x,b);
00593         }
00594   };
00595 
00596 
00597   // user calls
00598 
00600   template<class M, class X, class Y, class K>
00601   void dbgs (const M& A, X& x, const Y& b, const K& w)
00602   {
00603         algmeta_itsteps<1>::dbgs(A,x,b,w);
00604   }
00606   template<class M, class X, class Y, class K, int l>
00607   void dbgs (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
00608   {
00609         algmeta_itsteps<l>::dbgs(A,x,b,w);
00610   }
00612   template<class M, class X, class Y, class K>
00613   void bsorf (const M& A, X& x, const Y& b, const K& w)
00614   {
00615         algmeta_itsteps<1>::bsorf(A,x,b,w);
00616   }
00618   template<class M, class X, class Y, class K, int l>
00619   void bsorf (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
00620   {
00621         algmeta_itsteps<l>::bsorf(A,x,b,w);
00622   }
00624   template<class M, class X, class Y, class K>
00625   void bsorb (const M& A, X& x, const Y& b, const K& w)
00626   {
00627         algmeta_itsteps<1>::bsorb(A,x,b,w);
00628   }
00630   template<class M, class X, class Y, class K, int l>
00631   void bsorb (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
00632   {
00633         algmeta_itsteps<l>::bsorb(A,x,b,w);
00634   }
00636   template<class M, class X, class Y, class K>
00637   void dbjac (const M& A, X& x, const Y& b, const K& w)
00638   {
00639         algmeta_itsteps<1>::dbjac(A,x,b,w);
00640   }
00642   template<class M, class X, class Y, class K, int l>
00643   void dbjac (const M& A, X& x, const Y& b, const K& w, BL<l> bl)
00644   {
00645         algmeta_itsteps<l>::dbjac(A,x,b,w);
00646   }
00647 
00648 
00651 } // end namespace
00652 
00653 #endif