BelosCG.hpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 //
00004 //                 Belos: Block Linear Solvers Package
00005 //                 Copyright (2004) Sandia Corporation
00006 //
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 //
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00025 //
00026 // ***********************************************************************
00027 // @HEADER
00028 //
00029 // This file contains an implementation of the CG algorithm
00030 // for solving real symmetric positive definite linear systems of 
00031 // equations Ax = b, where b is a single-vector and x is the corresponding solution.
00032 // This implementation uses the preconditioned inner-product to retain symmetry of
00033 // the linear system.
00034 //
00035 #ifndef BELOS_CG_HPP
00036 #define BELOS_CG_HPP
00037 
00045 #include "BelosConfigDefs.hpp"
00046 #include "BelosIterativeSolver.hpp"
00047 #include "BelosLinearProblem.hpp"
00048 #include "BelosOutputManager.hpp"
00049 #include "BelosOperator.hpp"
00050 #include "BelosStatusTest.hpp"
00051 #include "BelosMultiVecTraits.hpp"
00052 
00053 #include "Teuchos_LAPACK.hpp"
00054 #include "Teuchos_SerialDenseMatrix.hpp"
00055 #include "Teuchos_ScalarTraits.hpp"
00056 #include "Teuchos_TimeMonitor.hpp"
00057 
00067 namespace Belos {
00068   
00069   template <class ScalarType, class MV, class OP>
00070   class CG : public IterativeSolver<ScalarType,MV,OP> { 
00071   public:
00072     //
00073     // Convenience typedefs
00074     //
00075     typedef MultiVecTraits<ScalarType,MV> MVT;
00076     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00077     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00078     typedef typename SCT::magnitudeType MagnitudeType;
00079     
00081 
00082 
00084     CG(const RefCountPtr<LinearProblem<ScalarType,MV,OP> > &lp, 
00085        const RefCountPtr<StatusTest<ScalarType,MV,OP> > &stest, 
00086        const RefCountPtr<OutputManager<ScalarType> > &om,
00087        const RefCountPtr<ParameterList> &pl = Teuchos::null
00088        );
00089     
00091     virtual ~CG() {};
00093     
00095 
00096     
00098     int GetNumIters() const { return( _iter ); }
00099     
00101     int GetNumRestarts() const { return(0); }
00102     
00104 
00107     RefCountPtr<const MV> GetNativeResiduals( std::vector<MagnitudeType> *normvec ) const;
00108     
00110 
00114     RefCountPtr<MV> GetCurrentSoln() { return MVT::CloneCopy( *_cur_block_sol ); };
00115   
00117 
00119     RefCountPtr<LinearProblem<ScalarType,MV,OP> > GetLinearProblem() const { return( _lp ); }
00120     
00121     RefCountPtr<StatusTest<ScalarType,MV,OP> > GetStatusTest() const { return( _stest ); }
00122 
00124     
00126 
00127     
00132     void Solve();
00134 
00137     
00139     std::string description() const;
00140     
00142     
00143   private:
00144     
00146     RefCountPtr<LinearProblem<ScalarType,MV,OP> > _lp; 
00147     
00149     RefCountPtr<StatusTest<ScalarType,MV,OP> > _stest; 
00150     
00152     RefCountPtr<OutputManager<ScalarType> > _om;
00153     
00155     RefCountPtr<ParameterList> _pl;     
00156 
00158     RefCountPtr<MV> _cur_block_sol;
00159     
00161     RefCountPtr<MV> _cur_block_rhs; 
00162     
00164     RefCountPtr<MV> _residvec;
00165     
00167     RefCountPtr<ostream> _os;
00168     
00170     int _iter;
00171     
00173     bool _restartTimers;
00174 
00176     Teuchos::RefCountPtr<Teuchos::Time> _timerOp, _timerPrec, _timerTotal;
00177   };
00178   
00179   //
00180   // Implementation
00181   //
00182   
00183   template <class ScalarType, class MV, class OP>
00184   CG<ScalarType,MV,OP>::CG(const RefCountPtr<LinearProblem<ScalarType,MV,OP> > &lp,
00185          const RefCountPtr<StatusTest<ScalarType,MV,OP> > &stest,
00186          const RefCountPtr<OutputManager<ScalarType> > &om,
00187          const RefCountPtr<ParameterList> &pl 
00188          ) : 
00189     _lp(lp), 
00190     _stest(stest),
00191     _om(om),
00192     _pl(pl),
00193     _os(om->GetOStream()),
00194     _iter(0),
00195     _restartTimers(true),
00196     _timerOp(Teuchos::TimeMonitor::getNewTimer("Operation Op*x")),
00197     _timerPrec(Teuchos::TimeMonitor::getNewTimer("Operation Prec*x")),
00198     _timerTotal(Teuchos::TimeMonitor::getNewTimer("Total time"))
00199   { 
00200   }
00201   
00202   template <class ScalarType, class MV, class OP>
00203   RefCountPtr<const MV> 
00204   CG<ScalarType,MV,OP>::GetNativeResiduals( std::vector<MagnitudeType> *normvec ) const 
00205   {
00206     std::vector<int> index( 1 );
00207     index[0] = 0;
00208     RefCountPtr<MV> ResidMV = MVT::CloneView( *_residvec, index );
00209     return ResidMV;
00210   }
00211   
00212   template <class ScalarType, class MV, class OP>
00213   void 
00214   CG<ScalarType,MV,OP>::Solve () 
00215   {
00216     Teuchos::TimeMonitor LocalTimer(*_timerTotal,_restartTimers);
00217 
00218     if ( _restartTimers ) {
00219       _timerOp->reset();
00220       _timerPrec->reset();
00221     }
00222 
00223     // Get the output stream from the OutputManager
00224     _os = _om->GetOStream();
00225     //
00226     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00227     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00228     bool exit_flg = false;
00229     bool isPrec = ( _lp->GetLeftPrec().get()!=NULL );
00230     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
00231     Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
00232     Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 ), pAp( 1, 1 );
00233     RefCountPtr<MV> _p, _Ap, _z;
00234     //
00235     // Retrieve the first linear system to be solved.
00236     //
00237     _cur_block_sol = _lp->GetCurrLHSVec();
00238     _cur_block_rhs = _lp->GetCurrRHSVec();
00239     //
00240     //  Start executable statements. 
00241     //
00242     while (_cur_block_sol.get() && _cur_block_rhs.get() ) {
00243       //
00244       // Only continue if the linear system is single-vector.
00245       //
00246       if ( _lp->GetBlockSize() > 1 ) return;
00247       //
00248       if (_om->isVerbosityAndPrint( IterationDetails )) {
00249   *_os << endl;
00250   *_os << "===================================================" << endl;
00251   *_os << "Solving linear system(s):  " << _lp->GetRHSIndex() << " through " << _lp->GetRHSIndex()+_lp->GetNumToSolve() << endl;
00252   *_os << endl;
00253       } 
00254       //
00255       _residvec = MVT::Clone( *_cur_block_sol, 1 );
00256       _p = MVT::CloneCopy( *_cur_block_sol );
00257       _Ap = MVT::Clone( *_cur_block_sol, 1 ); 
00258       _z = MVT::Clone( *_cur_block_sol, 1 ); 
00259       //
00260       // ************ Compute the initial residual ********************************
00261       //
00262       // p0 = r0 = cur_block_rhs - A * cur_block_sol 
00263       //
00264       // Multiply the current solution by A and store in _Ap
00265       //       _Ap = A*_p 
00266       //
00267       {
00268   Teuchos::TimeMonitor OpTimer(*_timerOp);
00269   _lp->ApplyOp( *_p, *_Ap );
00270       }
00271       //
00272       // Compute initial residual and store in _residvec
00273       //     _residvec = cur_block_rhs - _Ap
00274       //
00275       MVT::MvAddMv(one, *_cur_block_rhs, -one, *_Ap, *_residvec);
00276       //
00277       //----------------Compute initial direction vectors--------------------------
00278       // Initially, they are set to the preconditioned residuals
00279       //
00280       if ( isPrec ) {
00281   {
00282     Teuchos::TimeMonitor PrecTimer(*_timerPrec);
00283     _lp->ApplyLeftPrec( *_residvec, *_z ); 
00284   }
00285   MVT::MvAddMv( one, *_z, zero, *_z, *_p );
00286       } else {
00287   MVT::MvAddMv( one, *_residvec, zero, *_residvec, *_z ); 
00288   MVT::MvAddMv( one, *_residvec, zero, *_residvec, *_p );
00289       }
00290       //
00291       // Compute first <r,z> a.k.a. rHz
00292       // 
00293       MVT::MvTransMv( one, *_residvec, *_z, rHz );
00294       //
00295       // ***************************************************************************
00296       // ************************Main CG Loop***************************************
00297       // ***************************************************************************
00298       // 
00299       for (_iter=0; _stest->CheckStatus(this) == Unconverged && !exit_flg; _iter++) 
00300   {
00301     //
00302     // Multiply the current direction vector by A and store in _Ap
00303     //       _Ap = A*_p 
00304     //
00305     {
00306       Teuchos::TimeMonitor OpTimer(*_timerOp);
00307       _lp->ApplyOp( *_p, *_Ap );
00308     }
00309     //
00310     // Compute alpha := <_residvec, _z> / <_p, _Ap >
00311     //
00312     MVT::MvTransMv( one, *_p, *_Ap, pAp );
00313     //
00314     alpha(0,0) = rHz(0,0) / pAp(0,0);
00315     //
00316     // Check that alpha is a positive number!
00317     //
00318     if ( SCT::real(alpha(0,0)) <= zero ) {
00319       if (_om->isVerbosityAndPrint( Errors )) {
00320         *_os << " Exiting CG iteration " << endl;
00321         *_os << " ERROR: Non-positive value for p^H*A*p ("<< SCT::real(alpha(0,0)) <<") !!! "<< endl;
00322       }
00323       break; // Get out from this solve.
00324     }
00325     //
00326     // Update the solution vector x := x + alpha * _p
00327     //
00328     MVT::MvAddMv( one, *_cur_block_sol, alpha(0,0), *_p, *_cur_block_sol );
00329     _lp->SolutionUpdated();
00330     //
00331     // Save the denominator of beta before residual is updated [ old <_residvec, _z> ]
00332     //
00333     rHz_old(0,0) = rHz(0,0);
00334     //
00335     // Compute the new residual _residvec := _residvec - alpha * _Ap
00336     //
00337     MVT::MvAddMv( one, *_residvec, -alpha(0,0), *_Ap, *_residvec );
00338     //
00339     // Compute beta := [ new <_residvec, _z> ] / [ old <_residvec, _z> ], 
00340     // and the new direction vector p.
00341     //
00342     if ( isPrec ) {
00343       {
00344         Teuchos::TimeMonitor PrecTimer(*_timerPrec);
00345         _lp->ApplyLeftPrec( *_residvec, *_z );
00346       }
00347     } else {
00348       MVT::MvAddMv( one, *_residvec, zero, *_residvec, *_z ); 
00349     }
00350     //
00351     MVT::MvTransMv( one, *_residvec, *_z, rHz );
00352     //
00353     beta(0,0) = rHz(0,0) / rHz_old(0,0);
00354     //
00355     MVT::MvAddMv( one, *_z, beta(0,0), *_p, *_p );
00356     //
00357   } // end of the main CG loop -- for(_iter = 0;...)
00358       // *******************************************************************************
00359       //
00360       // Inform the linear problem manager that we are done with the current block of linear systems.
00361       //
00362       _lp->SetCurrLSVec();
00363       //
00364       // Get the next block of linear systems, if it returns the null pointer we are done.
00365       //
00366       _cur_block_sol = _lp->GetCurrLHSVec();
00367       _cur_block_rhs = _lp->GetCurrRHSVec();
00368       //
00369       // Print out solver status.
00370       //
00371       if (_om->isVerbosityAndPrint( FinalSummary )) {
00372   _stest->Print(*_os);
00373       }
00374       //
00375     } // end while ( _cur_block_sol && _cur_block_rhs )
00376     // **********************************************************************************
00377     //
00378     // Print timing details 
00379 
00380     // Stop timer.
00381     _timerTotal->stop();
00382 
00383     // Reset format that will be used to print the summary
00384     Teuchos::TimeMonitor::format().setPageWidth(54);
00385 
00386     if (_om->isVerbosity( Belos::TimingDetails )) {
00387       if (_om->doPrint())
00388         *_os <<"********************TIMING DETAILS********************"<<endl;
00389       Teuchos::TimeMonitor::summarize( *_os );
00390       if (_om->doPrint())
00391         *_os <<"******************************************************"<<endl;
00392     }
00393   } // end CGSolve()
00394 
00395 
00396   // Overridden from Teuchos::Describable
00397   template<class ScalarType, class MV, class OP>
00398   std::string 
00399   CG<ScalarType,MV,OP>::description() const
00400   {
00401     std::ostringstream oss;
00402     oss << "Belos::CG<...,"<<Teuchos::ScalarTraits<ScalarType>::name()<<">";
00403     oss << "{";
00404     oss << "}";
00405     return oss.str();
00406   }
00407 
00408   //
00409 } // namespace Belos
00410 //
00411 #endif // BELOS_CG_HPP
00412 //
00413 // End of file BelosCG.hpp
00414 
00415 

Generated on Thu Sep 18 12:30:12 2008 for Belos by doxygen 1.3.9.1