dune-istl
2.2.0
|
00001 // -*- tab-width: 4; indent-tabs-mode: nil; c-basic-offset: 2 -*- 00002 // vi: set ts=4 sw=2 et sts=2: 00003 00004 #ifndef DUNE_SOLVERS_HH 00005 #define DUNE_SOLVERS_HH 00006 00007 #include<cmath> 00008 #include<complex> 00009 #include<iostream> 00010 #include<iomanip> 00011 #include<string> 00012 00013 #include "istlexception.hh" 00014 #include "operators.hh" 00015 #include "preconditioners.hh" 00016 #include "scalarproducts.hh" 00017 #include <dune/common/timer.hh> 00018 #include <dune/common/ftraits.hh> 00019 #include <dune/common/static_assert.hh> 00020 00021 namespace Dune { 00046 struct InverseOperatorResult 00047 { 00049 InverseOperatorResult () 00050 { 00051 clear(); 00052 } 00053 00055 void clear () 00056 { 00057 iterations = 0; 00058 reduction = 0; 00059 converged = false; 00060 conv_rate = 1; 00061 elapsed = 0; 00062 } 00063 00065 int iterations; 00066 00068 double reduction; 00069 00071 bool converged; 00072 00074 double conv_rate; 00075 00077 double elapsed; 00078 }; 00079 00080 00081 //===================================================================== 00093 template<class X, class Y> 00094 class InverseOperator { 00095 public: 00097 typedef X domain_type; 00098 00100 typedef Y range_type; 00101 00103 typedef typename X::field_type field_type; 00104 00114 virtual void apply (X& x, Y& b, InverseOperatorResult& res) = 0; 00115 00126 virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) = 0; 00127 00129 virtual ~InverseOperator () {} 00130 00131 protected: 00132 // spacing values 00133 enum { iterationSpacing = 5 , normSpacing = 16 }; 00134 00136 void printHeader(std::ostream& s) const 00137 { 00138 s << std::setw(iterationSpacing) << " Iter"; 00139 s << std::setw(normSpacing) << "Defect"; 00140 s << std::setw(normSpacing) << "Rate" << std::endl; 00141 } 00142 00144 template <class DataType> 00145 void printOutput(std::ostream& s, 00146 const double iter, 00147 const DataType& norm, 00148 const DataType& norm_old) const 00149 { 00150 const DataType rate = norm/norm_old; 00151 s << std::setw(iterationSpacing) << iter << " "; 00152 s << std::setw(normSpacing) << norm << " "; 00153 s << std::setw(normSpacing) << rate << std::endl; 00154 } 00155 00157 template <class DataType> 00158 void printOutput(std::ostream& s, 00159 const double iter, 00160 const DataType& norm) const 00161 { 00162 s << std::setw(iterationSpacing) << iter << " "; 00163 s << std::setw(normSpacing) << norm << std::endl; 00164 } 00165 }; 00166 00167 00168 //===================================================================== 00169 // Implementation of this interface 00170 //===================================================================== 00171 00180 template<class X> 00181 class LoopSolver : public InverseOperator<X,X> { 00182 public: 00184 typedef X domain_type; 00186 typedef X range_type; 00188 typedef typename X::field_type field_type; 00189 00209 template<class L, class P> 00210 LoopSolver (L& op, P& prec, 00211 double reduction, int maxit, int verbose) : 00212 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) 00213 { 00214 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category), 00215 "L and P have to have the same category!"); 00216 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential), 00217 "L has to be sequential!"); 00218 } 00219 00240 template<class L, class S, class P> 00241 LoopSolver (L& op, S& sp, P& prec, 00242 double reduction, int maxit, int verbose) : 00243 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) 00244 { 00245 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category), 00246 "L and P must have the same category!"); 00247 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category), 00248 "L and S must have the same category!"); 00249 } 00250 00251 00253 virtual void apply (X& x, X& b, InverseOperatorResult& res) 00254 { 00255 // clear solver statistics 00256 res.clear(); 00257 00258 // start a timer 00259 Timer watch; 00260 00261 // prepare preconditioner 00262 _prec.pre(x,b); 00263 00264 // overwrite b with defect 00265 _op.applyscaleadd(-1,x,b); 00266 00267 // compute norm, \todo parallelization 00268 double def0 = _sp.norm(b); 00269 00270 // printing 00271 if (_verbose>0) 00272 { 00273 std::cout << "=== LoopSolver" << std::endl; 00274 if (_verbose>1) 00275 { 00276 this->printHeader(std::cout); 00277 this->printOutput(std::cout,0,def0); 00278 } 00279 } 00280 00281 // allocate correction vector 00282 X v(x); 00283 00284 // iteration loop 00285 int i=1; double def=def0; 00286 for ( ; i<=_maxit; i++ ) 00287 { 00288 v = 0; // clear correction 00289 _prec.apply(v,b); // apply preconditioner 00290 x += v; // update solution 00291 _op.applyscaleadd(-1,v,b); // update defect 00292 double defnew=_sp.norm(b); // comp defect norm 00293 if (_verbose>1) // print 00294 this->printOutput(std::cout,i,defnew,def); 00295 //std::cout << i << " " << defnew << " " << defnew/def << std::endl; 00296 def = defnew; // update norm 00297 if (def<def0*_reduction || def<1E-30) // convergence check 00298 { 00299 res.converged = true; 00300 break; 00301 } 00302 } 00303 00304 // print 00305 if (_verbose==1) 00306 this->printOutput(std::cout,i,def); 00307 00308 // postprocess preconditioner 00309 _prec.post(x); 00310 00311 // fill statistics 00312 res.iterations = i; 00313 res.reduction = def/def0; 00314 res.conv_rate = pow(res.reduction,1.0/i); 00315 res.elapsed = watch.elapsed(); 00316 00317 // final print 00318 if (_verbose>0) 00319 { 00320 std::cout << "=== rate=" << res.conv_rate 00321 << ", T=" << res.elapsed 00322 << ", TIT=" << res.elapsed/i 00323 << ", IT=" << i << std::endl; 00324 } 00325 } 00326 00328 virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) 00329 { 00330 std::swap(_reduction,reduction); 00331 (*this).apply(x,b,res); 00332 std::swap(_reduction,reduction); 00333 } 00334 00335 private: 00336 SeqScalarProduct<X> ssp; 00337 LinearOperator<X,X>& _op; 00338 Preconditioner<X,X>& _prec; 00339 ScalarProduct<X>& _sp; 00340 double _reduction; 00341 int _maxit; 00342 int _verbose; 00343 }; 00344 00345 00346 // all these solvers are taken from the SUMO library 00348 template<class X> 00349 class GradientSolver : public InverseOperator<X,X> { 00350 public: 00352 typedef X domain_type; 00354 typedef X range_type; 00356 typedef typename X::field_type field_type; 00357 00363 template<class L, class P> 00364 GradientSolver (L& op, P& prec, 00365 double reduction, int maxit, int verbose) : 00366 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) 00367 { 00368 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category), 00369 "L and P have to have the same category!"); 00370 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential), 00371 "L has to be sequential!"); 00372 } 00378 template<class L, class S, class P> 00379 GradientSolver (L& op, S& sp, P& prec, 00380 double reduction, int maxit, int verbose) : 00381 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) 00382 { 00383 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category), 00384 "L and P have to have the same category!"); 00385 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(S::category), 00386 "L and S have to have the same category!"); 00387 } 00388 00394 virtual void apply (X& x, X& b, InverseOperatorResult& res) 00395 { 00396 res.clear(); // clear solver statistics 00397 Timer watch; // start a timer 00398 _prec.pre(x,b); // prepare preconditioner 00399 _op.applyscaleadd(-1,x,b); // overwrite b with defect 00400 00401 X p(x); // create local vectors 00402 X q(b); 00403 00404 double def0 = _sp.norm(b);// compute norm 00405 00406 if (_verbose>0) // printing 00407 { 00408 std::cout << "=== GradientSolver" << std::endl; 00409 if (_verbose>1) 00410 { 00411 this->printHeader(std::cout); 00412 this->printOutput(std::cout,0,def0); 00413 } 00414 } 00415 00416 int i=1; double def=def0; // loop variables 00417 field_type lambda; 00418 for ( ; i<=_maxit; i++ ) 00419 { 00420 p = 0; // clear correction 00421 _prec.apply(p,b); // apply preconditioner 00422 _op.apply(p,q); // q=Ap 00423 lambda = _sp.dot(p,b)/_sp.dot(q,p);// minimization 00424 x.axpy(lambda,p); // update solution 00425 b.axpy(-lambda,q); // update defect 00426 00427 double defnew=_sp.norm(b);// comp defect norm 00428 if (_verbose>1) // print 00429 this->printOutput(std::cout,i,defnew,def); 00430 00431 def = defnew; // update norm 00432 if (def<def0*_reduction || def<1E-30) // convergence check 00433 { 00434 res.converged = true; 00435 break; 00436 } 00437 } 00438 00439 if (_verbose==1) // printing for non verbose 00440 this->printOutput(std::cout,i,def); 00441 00442 _prec.post(x); // postprocess preconditioner 00443 res.iterations = i; // fill statistics 00444 res.reduction = def/def0; 00445 res.conv_rate = pow(res.reduction,1.0/i); 00446 res.elapsed = watch.elapsed(); 00447 if (_verbose>0) // final print 00448 std::cout << "=== rate=" << res.conv_rate 00449 << ", T=" << res.elapsed 00450 << ", TIT=" << res.elapsed/i 00451 << ", IT=" << i << std::endl; 00452 } 00453 00459 virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) 00460 { 00461 std::swap(_reduction,reduction); 00462 (*this).apply(x,b,res); 00463 std::swap(_reduction,reduction); 00464 } 00465 00466 private: 00467 SeqScalarProduct<X> ssp; 00468 LinearOperator<X,X>& _op; 00469 Preconditioner<X,X>& _prec; 00470 ScalarProduct<X>& _sp; 00471 double _reduction; 00472 int _maxit; 00473 int _verbose; 00474 }; 00475 00476 00477 00479 template<class X> 00480 class CGSolver : public InverseOperator<X,X> { 00481 public: 00483 typedef X domain_type; 00485 typedef X range_type; 00487 typedef typename X::field_type field_type; 00488 00494 template<class L, class P> 00495 CGSolver (L& op, P& prec, double reduction, int maxit, int verbose) : 00496 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) 00497 { 00498 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category), 00499 "L and P must have the same category!"); 00500 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential), 00501 "L must be sequential!"); 00502 } 00508 template<class L, class S, class P> 00509 CGSolver (L& op, S& sp, P& prec, double reduction, int maxit, int verbose) : 00510 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) 00511 { 00512 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category), 00513 "L and P must have the same category!"); 00514 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category), 00515 "L and S must have the same category!"); 00516 } 00517 00523 virtual void apply (X& x, X& b, InverseOperatorResult& res) 00524 { 00525 res.clear(); // clear solver statistics 00526 Timer watch; // start a timer 00527 _prec.pre(x,b); // prepare preconditioner 00528 _op.applyscaleadd(-1,x,b); // overwrite b with defect 00529 00530 X p(x); // the search direction 00531 X q(x); // a temporary vector 00532 00533 double def0 = _sp.norm(b);// compute norm 00534 if (def0<1E-30) // convergence check 00535 { 00536 res.converged = true; 00537 res.iterations = 0; // fill statistics 00538 res.reduction = 0; 00539 res.conv_rate = 0; 00540 res.elapsed=0; 00541 if (_verbose>0) // final print 00542 std::cout << "=== rate=" << res.conv_rate 00543 << ", T=" << res.elapsed << ", TIT=" << res.elapsed 00544 << ", IT=0" << std::endl; 00545 return; 00546 } 00547 00548 if (_verbose>0) // printing 00549 { 00550 std::cout << "=== CGSolver" << std::endl; 00551 if (_verbose>1) { 00552 this->printHeader(std::cout); 00553 this->printOutput(std::cout,0,def0); 00554 } 00555 } 00556 00557 // some local variables 00558 double def=def0; // loop variables 00559 field_type rho,rholast,lambda,alpha,beta; 00560 00561 // determine initial search direction 00562 p = 0; // clear correction 00563 _prec.apply(p,b); // apply preconditioner 00564 rholast = _sp.dot(p,b); // orthogonalization 00565 00566 // the loop 00567 int i=1; 00568 for ( ; i<=_maxit; i++ ) 00569 { 00570 // minimize in given search direction p 00571 _op.apply(p,q); // q=Ap 00572 alpha = _sp.dot(p,q); // scalar product 00573 lambda = rholast/alpha; // minimization 00574 x.axpy(lambda,p); // update solution 00575 b.axpy(-lambda,q); // update defect 00576 00577 // convergence test 00578 double defnew=_sp.norm(b);// comp defect norm 00579 00580 if (_verbose>1) // print 00581 this->printOutput(std::cout,i,defnew,def); 00582 00583 def = defnew; // update norm 00584 if (def<def0*_reduction || def<1E-30) // convergence check 00585 { 00586 res.converged = true; 00587 break; 00588 } 00589 00590 // determine new search direction 00591 q = 0; // clear correction 00592 _prec.apply(q,b); // apply preconditioner 00593 rho = _sp.dot(q,b); // orthogonalization 00594 beta = rho/rholast; // scaling factor 00595 p *= beta; // scale old search direction 00596 p += q; // orthogonalization with correction 00597 rholast = rho; // remember rho for recurrence 00598 } 00599 00600 if (_verbose==1) // printing for non verbose 00601 this->printOutput(std::cout,i,def); 00602 00603 _prec.post(x); // postprocess preconditioner 00604 res.iterations = i; // fill statistics 00605 res.reduction = def/def0; 00606 res.conv_rate = pow(res.reduction,1.0/i); 00607 res.elapsed = watch.elapsed(); 00608 00609 if (_verbose>0) // final print 00610 { 00611 std::cout << "=== rate=" << res.conv_rate 00612 << ", T=" << res.elapsed 00613 << ", TIT=" << res.elapsed/i 00614 << ", IT=" << i << std::endl; 00615 } 00616 } 00617 00623 virtual void apply (X& x, X& b, double reduction, 00624 InverseOperatorResult& res) 00625 { 00626 std::swap(_reduction,reduction); 00627 (*this).apply(x,b,res); 00628 std::swap(_reduction,reduction); 00629 } 00630 00631 private: 00632 SeqScalarProduct<X> ssp; 00633 LinearOperator<X,X>& _op; 00634 Preconditioner<X,X>& _prec; 00635 ScalarProduct<X>& _sp; 00636 double _reduction; 00637 int _maxit; 00638 int _verbose; 00639 }; 00640 00641 00642 // Ronald Kriemanns BiCG-STAB implementation from Sumo 00644 template<class X> 00645 class BiCGSTABSolver : public InverseOperator<X,X> { 00646 public: 00648 typedef X domain_type; 00650 typedef X range_type; 00652 typedef typename X::field_type field_type; 00653 00659 template<class L, class P> 00660 BiCGSTABSolver (L& op, P& prec, 00661 double reduction, int maxit, int verbose) : 00662 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) 00663 { 00664 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(P::category), "L and P must be of the same category!"); 00665 dune_static_assert(static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential), "L must be sequential!"); 00666 } 00672 template<class L, class S, class P> 00673 BiCGSTABSolver (L& op, S& sp, P& prec, 00674 double reduction, int maxit, int verbose) : 00675 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) 00676 { 00677 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category), 00678 "L and P must have the same category!"); 00679 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category), 00680 "L and S must have the same category!"); 00681 } 00682 00688 virtual void apply (X& x, X& b, InverseOperatorResult& res) 00689 { 00690 const double EPSILON=1e-80; 00691 00692 double it; 00693 field_type rho, rho_new, alpha, beta, h, omega; 00694 field_type norm, norm_old, norm_0; 00695 00696 // 00697 // get vectors and matrix 00698 // 00699 X& r=b; 00700 X p(x); 00701 X v(x); 00702 X t(x); 00703 X y(x); 00704 X rt(x); 00705 00706 // 00707 // begin iteration 00708 // 00709 00710 // r = r - Ax; rt = r 00711 res.clear(); // clear solver statistics 00712 Timer watch; // start a timer 00713 _prec.pre(x,r); // prepare preconditioner 00714 _op.applyscaleadd(-1,x,r); // overwrite b with defect 00715 00716 rt=r; 00717 00718 norm = norm_old = norm_0 = _sp.norm(r); 00719 00720 p=0; 00721 v=0; 00722 00723 rho = 1; 00724 alpha = 1; 00725 omega = 1; 00726 00727 if (_verbose>0) // printing 00728 { 00729 std::cout << "=== BiCGSTABSolver" << std::endl; 00730 if (_verbose>1) 00731 { 00732 this->printHeader(std::cout); 00733 this->printOutput(std::cout,0,norm_0); 00734 //std::cout << " Iter Defect Rate" << std::endl; 00735 //std::cout << " 0" << std::setw(14) << norm_0 << std::endl; 00736 } 00737 } 00738 00739 if ( norm < (_reduction * norm_0) || norm<1E-30) 00740 { 00741 res.converged = 1; 00742 _prec.post(x); // postprocess preconditioner 00743 res.iterations = 0; // fill statistics 00744 res.reduction = 0; 00745 res.conv_rate = 0; 00746 res.elapsed = watch.elapsed(); 00747 return; 00748 } 00749 00750 // 00751 // iteration 00752 // 00753 00754 for (it = 0.5; it < _maxit; it+=.5) 00755 { 00756 // 00757 // preprocess, set vecsizes etc. 00758 // 00759 00760 // rho_new = < rt , r > 00761 rho_new = _sp.dot(rt,r); 00762 00763 // look if breakdown occured 00764 if (std::abs(rho) <= EPSILON) 00765 DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - rho " 00766 << rho << " <= EPSILON " << EPSILON 00767 << " after " << it << " iterations"); 00768 if (std::abs(omega) <= EPSILON) 00769 DUNE_THROW(ISTLError,"breakdown in BiCGSTAB - omega " 00770 << omega << " <= EPSILON " << EPSILON 00771 << " after " << it << " iterations"); 00772 00773 00774 if (it<1) 00775 p = r; 00776 else 00777 { 00778 beta = ( rho_new / rho ) * ( alpha / omega ); 00779 p.axpy(-omega,v); // p = r + beta (p - omega*v) 00780 p *= beta; 00781 p += r; 00782 } 00783 00784 // y = W^-1 * p 00785 y = 0; 00786 _prec.apply(y,p); // apply preconditioner 00787 00788 // v = A * y 00789 _op.apply(y,v); 00790 00791 // alpha = rho_new / < rt, v > 00792 h = _sp.dot(rt,v); 00793 00794 if ( std::abs(h) < EPSILON ) 00795 DUNE_THROW(ISTLError,"h=0 in BiCGSTAB"); 00796 00797 alpha = rho_new / h; 00798 00799 // apply first correction to x 00800 // x <- x + alpha y 00801 x.axpy(alpha,y); 00802 00803 // r = r - alpha*v 00804 r.axpy(-alpha,v); 00805 00806 // 00807 // test stop criteria 00808 // 00809 00810 norm = _sp.norm(r); 00811 00812 if (_verbose>1) // print 00813 { 00814 this->printOutput(std::cout,it,norm,norm_old); 00815 } 00816 00817 if ( norm < (_reduction * norm_0) ) 00818 { 00819 res.converged = 1; 00820 break; 00821 } 00822 it+=.5; 00823 00824 norm_old = norm; 00825 00826 // y = W^-1 * r 00827 y = 0; 00828 _prec.apply(y,r); 00829 00830 // t = A * y 00831 _op.apply(y,t); 00832 00833 // omega = < t, r > / < t, t > 00834 omega = _sp.dot(t,r)/_sp.dot(t,t); 00835 00836 // apply second correction to x 00837 // x <- x + omega y 00838 x.axpy(omega,y); 00839 00840 // r = s - omega*t (remember : r = s) 00841 r.axpy(-omega,t); 00842 00843 rho = rho_new; 00844 00845 // 00846 // test stop criteria 00847 // 00848 00849 norm = _sp.norm(r); 00850 00851 if (_verbose > 1) // print 00852 { 00853 this->printOutput(std::cout,it,norm,norm_old); 00854 } 00855 00856 if ( norm < (_reduction * norm_0) || norm<1E-30) 00857 { 00858 res.converged = 1; 00859 break; 00860 } 00861 00862 norm_old = norm; 00863 } // end for 00864 00865 if (_verbose==1) // printing for non verbose 00866 this->printOutput(std::cout,it,norm); 00867 00868 _prec.post(x); // postprocess preconditioner 00869 res.iterations = static_cast<int>(std::ceil(it)); // fill statistics 00870 res.reduction = norm/norm_0; 00871 res.conv_rate = pow(res.reduction,1.0/it); 00872 res.elapsed = watch.elapsed(); 00873 if (_verbose>0) // final print 00874 std::cout << "=== rate=" << res.conv_rate 00875 << ", T=" << res.elapsed 00876 << ", TIT=" << res.elapsed/it 00877 << ", IT=" << it << std::endl; 00878 } 00879 00885 virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) 00886 { 00887 std::swap(_reduction,reduction); 00888 (*this).apply(x,b,res); 00889 std::swap(_reduction,reduction); 00890 } 00891 00892 private: 00893 SeqScalarProduct<X> ssp; 00894 LinearOperator<X,X>& _op; 00895 Preconditioner<X,X>& _prec; 00896 ScalarProduct<X>& _sp; 00897 double _reduction; 00898 int _maxit; 00899 int _verbose; 00900 }; 00901 00908 template<class X> 00909 class MINRESSolver : public InverseOperator<X,X> { 00910 public: 00912 typedef X domain_type; 00914 typedef X range_type; 00916 typedef typename X::field_type field_type; 00918 typedef typename FieldTraits<field_type>::real_type real_type; 00919 00925 template<class L, class P> 00926 MINRESSolver (L& op, P& prec, double reduction, int maxit, int verbose) : 00927 ssp(), _op(op), _prec(prec), _sp(ssp), _reduction(reduction), _maxit(maxit), _verbose(verbose) 00928 { 00929 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category), 00930 "L and P must have the same category!"); 00931 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential), 00932 "L must be sequential!"); 00933 } 00939 template<class L, class S, class P> 00940 MINRESSolver (L& op, S& sp, P& prec, double reduction, int maxit, int verbose) : 00941 _op(op), _prec(prec), _sp(sp), _reduction(reduction), _maxit(maxit), _verbose(verbose) 00942 { 00943 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(P::category), 00944 "L and P must have the same category!"); 00945 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(S::category), 00946 "L and S must have the same category!"); 00947 } 00948 00954 virtual void apply (X& x, X& b, InverseOperatorResult& res) 00955 { 00956 res.clear(); // clear solver statistics 00957 Timer watch; // start a timer 00958 _prec.pre(x,b); // prepare preconditioner 00959 _op.applyscaleadd(-1,x,b); // overwrite b with defect/residual 00960 00961 real_type def0 = _sp.norm(b); // compute residual norm 00962 00963 if (def0<1E-30) // convergence check 00964 { 00965 res.converged = true; 00966 res.iterations = 0; // fill statistics 00967 res.reduction = 0; 00968 res.conv_rate = 0; 00969 res.elapsed=0; 00970 if (_verbose>0) // final print 00971 std::cout << "=== rate=" << res.conv_rate << ", T=" << res.elapsed << ", TIT=" << res.elapsed << ", IT=0" << std::endl; 00972 return; 00973 } 00974 00975 if (_verbose>0) // printing 00976 { 00977 std::cout << "=== MINRESSolver" << std::endl; 00978 if (_verbose>1) { 00979 this->printHeader(std::cout); 00980 this->printOutput(std::cout,0,def0); 00981 } 00982 } 00983 00984 // some local variables 00985 real_type def=def0; // the defect/residual norm 00986 field_type alpha, // recurrence coefficients as computed in the Lanczos alg making up the matrix T 00987 c[2]={0.0, 0.0}, // diagonal entry of Givens rotation 00988 s[2]={0.0, 0.0}; // off-diagonal entries of Givens rotation 00989 real_type beta; 00990 00991 field_type T[3]={0.0, 0.0, 0.0}; // recurrence coefficients (column k of Matrix T) 00992 00993 X z(b.size()), // some temporary vectors 00994 dummy(b.size()); 00995 00996 field_type xi[2]={1.0, 0.0}; 00997 00998 // initialize 00999 z = 0.0; // clear correction 01000 01001 _prec.apply(z,b); // apply preconditioner z=M^-1*b 01002 01003 beta = std::sqrt(std::abs(_sp.dot(z,b))); 01004 real_type beta0 = beta; 01005 01006 X p[3]; // the search directions 01007 X q[3]; // Orthonormal basis vectors (in unpreconditioned case) 01008 01009 q[0].resize(b.size()); 01010 q[1].resize(b.size()); 01011 q[2].resize(b.size()); 01012 q[0] = 0.0; 01013 q[1] = b; 01014 q[1] /= beta; 01015 q[2] = 0.0; 01016 01017 p[0].resize(b.size()); 01018 p[1].resize(b.size()); 01019 p[2].resize(b.size()); 01020 p[0] = 0.0; 01021 p[1] = 0.0; 01022 p[2] = 0.0; 01023 01024 01025 z /= beta; // this is w_current 01026 01027 // the loop 01028 int i=1; 01029 for ( ; i<=_maxit; i++) 01030 { 01031 dummy = z; // remember z_old for the computation of the search direction p in the next iteration 01032 01033 int i1 = i%3, 01034 i0 = (i1+2)%3, 01035 i2 = (i1+1)%3; 01036 01037 // Symmetrically Preconditioned Lanczos (Greenbaum p.121) 01038 _op.apply(z,q[i2]); // q[i2] = Az 01039 q[i2].axpy(-beta, q[i0]); 01040 alpha = _sp.dot(q[i2],z); 01041 q[i2].axpy(-alpha, q[i1]); 01042 01043 z=0.0; 01044 _prec.apply(z,q[i2]); 01045 01046 beta = std::sqrt(std::abs(_sp.dot(q[i2],z))); 01047 01048 q[i2] /= beta; 01049 z /= beta; 01050 01051 // QR Factorization of recurrence coefficient matrix 01052 // apply previous Givens rotations to last column of T 01053 T[1] = T[2]; 01054 if (i>2) 01055 { 01056 T[0] = s[i%2]*T[1]; 01057 T[1] = c[i%2]*T[1]; 01058 } 01059 if (i>1) 01060 { 01061 T[2] = c[(i+1)%2]*alpha - s[(i+1)%2]*T[1]; 01062 T[1] = c[(i+1)%2]*T[1] + s[(i+1)%2]*alpha; 01063 } 01064 else 01065 T[2] = alpha; 01066 01067 // recompute c, s -> current Givens rotation \TODO use BLAS-routine drotg instead for greater robustness 01068 // cblas_drotg (a, b, c, s); 01069 c[i%2] = 1.0/std::sqrt(T[2]*T[2] + beta*beta); 01070 s[i%2] = beta*c[i%2]; 01071 c[i%2] *= T[2]; 01072 01073 // apply current Givens rotation to T eliminating the last entry... 01074 T[2] = c[i%2]*T[2] + s[i%2]*beta; 01075 01076 // ...and to xi, the right hand side of the least squares problem min_y||beta*xi-T*y|| 01077 xi[i%2] = -s[i%2]*xi[(i+1)%2]; 01078 xi[(i+1)%2] *= c[i%2]; 01079 01080 // compute correction direction 01081 p[i2] = dummy; 01082 p[i2].axpy(-T[1],p[i1]); 01083 p[i2].axpy(-T[0],p[i0]); 01084 p[i2] /= T[2]; 01085 01086 // apply correction/update solution 01087 x.axpy(beta0*xi[(i+1)%2], p[i2]); 01088 01089 // remember beta_old 01090 T[2] = beta; 01091 01092 // update residual - not necessary if in the preconditioned case we are content with the residual norm of the 01093 // preconditioned system as convergence test 01094 // _op.apply(p[i2],dummy); 01095 // b.axpy(-beta0*xi[(i+1)%2],dummy); 01096 01097 // convergence test 01098 real_type defnew = std::abs(beta0*xi[i%2]); // the last entry the QR-transformed least squares RHS is the new residual norm 01099 01100 if (_verbose>1) // print 01101 this->printOutput(std::cout,i,defnew,def); 01102 01103 def = defnew; // update norm 01104 if (def<def0*_reduction || def<1E-30 || i==_maxit) // convergence check 01105 { 01106 res.converged = true; 01107 break; 01108 } 01109 } 01110 01111 if (_verbose==1) // printing for non verbose 01112 this->printOutput(std::cout,i,def); 01113 01114 _prec.post(x); // postprocess preconditioner 01115 res.iterations = i; // fill statistics 01116 res.reduction = def/def0; 01117 res.conv_rate = pow(res.reduction,1.0/i); 01118 res.elapsed = watch.elapsed(); 01119 01120 if (_verbose>0) // final print 01121 { 01122 std::cout << "=== rate=" << res.conv_rate 01123 << ", T=" << res.elapsed 01124 << ", TIT=" << res.elapsed/i 01125 << ", IT=" << i << std::endl; 01126 } 01127 01128 } 01129 01135 virtual void apply (X& x, X& b, double reduction, InverseOperatorResult& res) 01136 { 01137 std::swap(_reduction,reduction); 01138 (*this).apply(x,b,res); 01139 std::swap(_reduction,reduction); 01140 } 01141 01142 private: 01143 SeqScalarProduct<X> ssp; 01144 LinearOperator<X,X>& _op; 01145 Preconditioner<X,X>& _prec; 01146 ScalarProduct<X>& _sp; 01147 double _reduction; 01148 int _maxit; 01149 int _verbose; 01150 }; 01151 01163 template<class X, class Y=X, class F = Y> 01164 class RestartedGMResSolver : public InverseOperator<X,Y> 01165 { 01166 public: 01168 typedef X domain_type; 01170 typedef Y range_type; 01172 typedef typename X::field_type field_type; 01174 typedef F basis_type; 01175 01183 template<class L, class P> 01184 RestartedGMResSolver (L& op, P& prec, double reduction, int restart, int maxit, int verbose, bool recalc_defect = false) : 01185 _A_(op), _M(prec), 01186 ssp(), _sp(ssp), _restart(restart), 01187 _reduction(reduction), _maxit(maxit), _verbose(verbose), 01188 _recalc_defect(recalc_defect) 01189 { 01190 dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category), 01191 "P and L must be the same category!"); 01192 dune_static_assert( static_cast<int>(L::category) == static_cast<int>(SolverCategory::sequential), 01193 "L must be sequential!"); 01194 } 01195 01203 template<class L, class S, class P> 01204 RestartedGMResSolver (L& op, S& sp, P& prec, double reduction, int restart, int maxit, int verbose, bool recalc_defect = false) : 01205 _A_(op), _M(prec), 01206 _sp(sp), _restart(restart), 01207 _reduction(reduction), _maxit(maxit), _verbose(verbose), 01208 _recalc_defect(recalc_defect) 01209 { 01210 dune_static_assert(static_cast<int>(P::category) == static_cast<int>(L::category), 01211 "P and L must have the same category!"); 01212 dune_static_assert(static_cast<int>(P::category) == static_cast<int>(S::category), 01213 "P and S must have the same category!"); 01214 } 01215 01217 virtual void apply (X& x, X& b, InverseOperatorResult& res) 01218 { 01219 apply(x,b,_reduction,res); 01220 }; 01221 01227 virtual void apply (X& x, Y& b, double reduction, InverseOperatorResult& res) 01228 { 01229 int m = _restart; 01230 field_type norm; 01231 field_type norm_old = 0.0; 01232 field_type norm_0; 01233 field_type beta; 01234 int i, j = 1, k; 01235 std::vector<field_type> s(m+1), cs(m), sn(m); 01236 // helper vector 01237 X w(b); 01238 std::vector< std::vector<field_type> > H(m+1,s); 01239 std::vector<F> v(m+1,b); 01240 01241 // start timer 01242 Timer watch; // start a timer 01243 01244 // clear solver statistics 01245 res.clear(); 01246 _M.pre(x,b); 01247 if (_recalc_defect) 01248 { 01249 // norm_0 = norm(M^-1 b) 01250 w = 0.0; _M.apply(w,b); // w = M^-1 b 01251 norm_0 = _sp.norm(w); 01252 // r = _M.solve(b - A * x); 01253 w = b; 01254 _A_.applyscaleadd(-1,x, /* => */ w); // w = b - Ax; 01255 v[0] = 0.0; _M.apply(v[0],w); // r = M^-1 w 01256 beta = _sp.norm(v[0]); 01257 } 01258 else 01259 { 01260 // norm_0 = norm(M^-1 b) 01261 w = 0.0; _M.apply(w,b); // w = M^-1 b 01262 // r = _M.solve(b - A * x); 01263 _A_.applyscaleadd(-1,x, /* => */ b); // b = b - Ax; 01264 norm_0 = _sp.norm(b); 01265 v[0] = 0.0; _M.apply(v[0],b); // r = M^-1 b 01266 beta = _sp.norm(v[0]); 01267 } 01268 01269 // avoid division by zero 01270 if (norm_0 == 0.0) 01271 norm_0 = 1.0; 01272 norm = norm_old = _sp.norm(v[0]); 01273 01274 // print header 01275 if (_verbose > 0) 01276 { 01277 std::cout << "=== RestartedGMResSolver" << std::endl; 01278 if (_verbose > 1) 01279 { 01280 this->printHeader(std::cout); 01281 this->printOutput(std::cout,0,norm_0); 01282 this->printOutput(std::cout,0,norm, norm_0); 01283 } 01284 } 01285 01286 // check convergence 01287 if (norm <= reduction * norm_0) { 01288 _M.post(x); // postprocess preconditioner 01289 res.converged = true; 01290 if (_verbose > 0) // final print 01291 print_result(res); 01292 return; 01293 } 01294 01295 while (j <= _maxit && res.converged != true) { 01296 v[0] *= (1.0 / beta); 01297 for (i=1; i<=m; i++) s[i] = 0.0; 01298 s[0] = beta; 01299 01300 for (i = 0; i < m && j <= _maxit && res.converged != true; i++, j++) { 01301 w = 0.0; 01302 v[i+1] = 0.0; // use v[i+1] as temporary vector 01303 _A_.apply(v[i], /* => */ v[i+1]); 01304 _M.apply(w, v[i+1]); 01305 for (k = 0; k <= i; k++) { 01306 H[k][i] = _sp.dot(w, v[k]); 01307 // w -= H[k][i] * v[k]; 01308 w.axpy(-H[k][i], v[k]); 01309 } 01310 H[i+1][i] = _sp.norm(w); 01311 if (H[i+1][i] == 0.0) 01312 DUNE_THROW(ISTLError,"breakdown in GMRes - |w| " 01313 << w << " == 0.0 after " << j << " iterations"); 01314 // v[i+1] = w * (1.0 / H[i+1][i]); 01315 v[i+1] = w; v[i+1] *= (1.0 / H[i+1][i]); 01316 01317 for (k = 0; k < i; k++) 01318 applyPlaneRotation(H[k][i], H[k+1][i], cs[k], sn[k]); 01319 01320 generatePlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]); 01321 applyPlaneRotation(H[i][i], H[i+1][i], cs[i], sn[i]); 01322 applyPlaneRotation(s[i], s[i+1], cs[i], sn[i]); 01323 01324 norm = std::abs(s[i+1]); 01325 01326 if (_verbose > 1) // print 01327 { 01328 this->printOutput(std::cout,j,norm,norm_old); 01329 } 01330 01331 norm_old = norm; 01332 01333 if (norm < reduction * norm_0) { 01334 res.converged = true; 01335 } 01336 } 01337 01338 if (_recalc_defect) 01339 { 01340 // update x 01341 update(x, i - 1, H, s, v); 01342 01343 // update residuum 01344 // r = M^-1 (b - A * x); 01345 w = b; _A_.applyscaleadd(-1,x, /* => */ w); 01346 _M.apply(v[0], w); 01347 beta = _sp.norm(v[0]); 01348 norm = beta; 01349 } 01350 else 01351 { 01352 // calc update vector 01353 w = 0; 01354 update(w, i - 1, H, s, v); 01355 01356 // update x 01357 x += w; 01358 01359 // r = M^-1 (b - A * x); 01360 // update defect 01361 _A_.applyscaleadd(-1,w, /* => */ b); 01362 // r = M^-1 (b - A * x); 01363 v[0] = 0.0; _M.apply(v[0],b); // r = M^-1 b 01364 beta = _sp.norm(v[0]); 01365 norm = beta; 01366 01367 res.converged = false; 01368 } 01369 01370 if (_verbose > 1) // print 01371 { 01372 this->printOutput(std::cout,j,norm,norm_old); 01373 } 01374 01375 norm_old = norm; 01376 01377 if (norm < reduction * norm_0) { 01378 // fill statistics 01379 res.converged = true; 01380 } 01381 01382 if (res.converged != true && _verbose > 0) 01383 std::cout << "=== GMRes::restart\n"; 01384 } 01385 01386 _M.post(x); // postprocess preconditioner 01387 01388 res.iterations = j; 01389 res.reduction = norm / norm_0; 01390 res.conv_rate = pow(res.reduction,1.0/j); 01391 res.elapsed = watch.elapsed(); 01392 01393 if (_verbose>0) 01394 print_result(res); 01395 } 01396 private: 01397 01398 void 01399 print_result (const InverseOperatorResult & res) const 01400 { 01401 int j = res.iterations>0?res.iterations:1; 01402 std::cout << "=== rate=" << res.conv_rate 01403 << ", T=" << res.elapsed 01404 << ", TIT=" << res.elapsed/j 01405 << ", IT=" << res.iterations 01406 << std::endl; 01407 } 01408 01409 static void 01410 update(X &x, int k, 01411 std::vector< std::vector<field_type> > & h, 01412 std::vector<field_type> & s, std::vector<F> v) 01413 { 01414 std::vector<field_type> y(s); 01415 01416 // Backsolve: 01417 for (int i = k; i >= 0; i--) { 01418 y[i] /= h[i][i]; 01419 for (int j = i - 1; j >= 0; j--) 01420 y[j] -= h[j][i] * y[i]; 01421 } 01422 01423 for (int j = 0; j <= k; j++) 01424 // x += v[j] * y[j]; 01425 x.axpy(y[j],v[j]); 01426 } 01427 01428 void 01429 generatePlaneRotation(field_type &dx, field_type &dy, field_type &cs, field_type &sn) 01430 { 01431 if (dy == 0.0) { 01432 cs = 1.0; 01433 sn = 0.0; 01434 } else if (std::abs(dy) > std::abs(dx)) { 01435 field_type temp = dx / dy; 01436 sn = 1.0 / std::sqrt( 1.0 + temp*temp ); 01437 cs = temp * sn; 01438 } else { 01439 field_type temp = dy / dx; 01440 cs = 1.0 / std::sqrt( 1.0 + temp*temp ); 01441 sn = temp * cs; 01442 } 01443 } 01444 01445 01446 void 01447 applyPlaneRotation(field_type &dx, field_type &dy, field_type &cs, field_type &sn) 01448 { 01449 field_type temp = cs * dx + sn * dy; 01450 dy = -sn * dx + cs * dy; 01451 dx = temp; 01452 } 01453 01454 LinearOperator<X,X>& _A_; 01455 Preconditioner<X,X>& _M; 01456 SeqScalarProduct<X> ssp; 01457 ScalarProduct<X>& _sp; 01458 int _restart; 01459 double _reduction; 01460 int _maxit; 01461 int _verbose; 01462 bool _recalc_defect; 01463 }; 01464 01467 } // end namespace 01468 01469 #endif