dune-istl
2.2.0
|
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