BelosCGIter.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 #ifndef BELOS_CG_ITER_HPP
00030 #define BELOS_CG_ITER_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 #include "BelosCGIteration.hpp"
00039 
00040 #include "BelosLinearProblem.hpp"
00041 #include "BelosOutputManager.hpp"
00042 #include "BelosStatusTest.hpp"
00043 #include "BelosOperatorTraits.hpp"
00044 #include "BelosMultiVecTraits.hpp"
00045 
00046 #include "Teuchos_SerialDenseMatrix.hpp"
00047 #include "Teuchos_SerialDenseVector.hpp"
00048 #include "Teuchos_ScalarTraits.hpp"
00049 #include "Teuchos_ParameterList.hpp"
00050 #include "Teuchos_TimeMonitor.hpp"
00051 
00062 namespace Belos {
00063   
00064 template<class ScalarType, class MV, class OP>
00065 class CGIter : virtual public CGIteration<ScalarType,MV,OP> {
00066 
00067   public:
00068     
00069   //
00070   // Convenience typedefs
00071   //
00072   typedef MultiVecTraits<ScalarType,MV> MVT;
00073   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00074   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00075   typedef typename SCT::magnitudeType MagnitudeType;
00076 
00078 
00079 
00085   CGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00086       const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00087       const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00088       Teuchos::ParameterList &params );
00089 
00091   virtual ~CGIter() {};
00093 
00094 
00096 
00097   
00110   void iterate();
00111 
00126   void initializeCG(CGIterationState<ScalarType,MV> newstate);
00127 
00131   void initialize()
00132   {
00133     CGIterationState<ScalarType,MV> empty;
00134     initializeCG(empty);
00135   }
00136   
00143   CGIterationState<ScalarType,MV> getState() const {
00144     CGIterationState<ScalarType,MV> state;
00145     state.R = R_;
00146     state.P = P_;
00147     state.AP = AP_;
00148     state.Z = Z_;
00149     return state;
00150   }
00151 
00153 
00154   
00156 
00157 
00159   int getNumIters() const { return iter_; }
00160   
00162   void resetNumIters( int iter = 0 ) { iter_ = iter; }
00163 
00166   Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00167 
00169 
00171   Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00172 
00174   
00176 
00177 
00179   const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00180 
00182   int getBlockSize() const { return 1; }
00183 
00185   void setBlockSize(int blockSize) {
00186     TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00187            "Belos::CGIter::setBlockSize(): Cannot use a block size that is not one.");
00188   }
00189 
00191   bool isInitialized() { return initialized_; }
00192 
00194 
00195   private:
00196 
00197   //
00198   // Internal methods
00199   //
00201   void setStateSize();
00202   
00203   //
00204   // Classes inputed through constructor that define the linear problem to be solved.
00205   //
00206   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00207   const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00208   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00209 
00210   //  
00211   // Current solver state
00212   //
00213   // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00214   // is capable of running; _initialize is controlled  by the initialize() member method
00215   // For the implications of the state of initialized_, please see documentation for initialize()
00216   bool initialized_;
00217 
00218   // stateStorageInitialized_ specifies that the state storage has been initialized.
00219   // This initialization may be postponed if the linear problem was generated without 
00220   // the right-hand side or solution vectors.
00221   bool stateStorageInitialized_;
00222 
00223   // Current number of iterations performed.
00224   int iter_;
00225   
00226   // 
00227   // State Storage
00228   // 
00229   // Residual
00230   Teuchos::RCP<MV> R_;
00231   //
00232   // Preconditioned residual
00233   Teuchos::RCP<MV> Z_;
00234   //
00235   // Direction vector
00236   Teuchos::RCP<MV> P_;
00237   //
00238   // Operator applied to direction vector
00239   Teuchos::RCP<MV> AP_;
00240 
00241 };
00242 
00244   // Constructor.
00245   template<class ScalarType, class MV, class OP>
00246   CGIter<ScalarType,MV,OP>::CGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00247                const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00248                const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00249                Teuchos::ParameterList &params ):
00250     lp_(problem),
00251     om_(printer),
00252     stest_(tester),
00253     initialized_(false),
00254     stateStorageInitialized_(false),
00255     iter_(0)
00256   {
00257   }
00258 
00260   // Setup the state storage.
00261   template <class ScalarType, class MV, class OP>
00262   void CGIter<ScalarType,MV,OP>::setStateSize ()
00263   {
00264     if (!stateStorageInitialized_) {
00265 
00266       // Check if there is any multivector to clone from.
00267       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00268       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00269       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00270   stateStorageInitialized_ = false;
00271   return;
00272       }
00273       else {
00274   
00275   // Initialize the state storage
00276   // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00277   if (R_ == Teuchos::null) {
00278     // Get the multivector that is not null.
00279     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00280     TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00281            "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00282     R_ = MVT::Clone( *tmp, 1 );
00283     Z_ = MVT::Clone( *tmp, 1 );
00284     P_ = MVT::Clone( *tmp, 1 );
00285     AP_ = MVT::Clone( *tmp, 1 );
00286   }
00287   
00288   // State storage has now been initialized.
00289   stateStorageInitialized_ = true;
00290       }
00291     }
00292   }
00293 
00294 
00296   // Initialize this iteration object
00297   template <class ScalarType, class MV, class OP>
00298   void CGIter<ScalarType,MV,OP>::initializeCG(CGIterationState<ScalarType,MV> newstate)
00299   {
00300     // Initialize the state storage if it isn't already.
00301     if (!stateStorageInitialized_) 
00302       setStateSize();
00303 
00304     TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00305            "Belos::CGIter::initialize(): Cannot initialize state storage!");
00306     
00307     // NOTE:  In CGIter R_, the initial residual, is required!!!  
00308     //
00309     std::string errstr("Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
00310 
00311     // Create convenience variables for zero and one.
00312     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00313     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00314 
00315     if (newstate.R != Teuchos::null) {
00316 
00317       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
00318                           std::invalid_argument, errstr );
00319       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
00320                           std::invalid_argument, errstr );
00321 
00322       // Copy basis vectors from newstate into V
00323       if (newstate.R != R_) {
00324         // copy over the initial residual (unpreconditioned).
00325   MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
00326       }
00327 
00328       // Compute initial direction vectors
00329       // Initially, they are set to the preconditioned residuals
00330       //
00331       if ( lp_->getLeftPrec() != Teuchos::null ) {
00332   lp_->applyLeftPrec( *R_, *Z_ ); 
00333   MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00334       } else {
00335   Z_ = R_;
00336   MVT::MvAddMv( one, *R_, zero, *R_, *P_ );
00337       }
00338     }
00339     else {
00340 
00341       TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00342                          "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
00343     }
00344 
00345     // The solver is initialized
00346     initialized_ = true;
00347   }
00348 
00349 
00351   // Iterate until the status test informs us we should stop.
00352   template <class ScalarType, class MV, class OP>
00353   void CGIter<ScalarType,MV,OP>::iterate()
00354   {
00355     //
00356     // Allocate/initialize data structures
00357     //
00358     if (initialized_ == false) {
00359       initialize();
00360     }
00361 
00362     // Allocate memory for scalars.
00363     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
00364     Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
00365     Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 ), pAp( 1, 1 );
00366 
00367     // Create convenience variables for zero and one.
00368     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00369     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00370     
00371     // Get the current solution vector.
00372     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00373 
00374     // Check that the current solution vector only has one column. 
00375     TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
00376                         "Belos::CGIter::iterate(): current linear system has more than one vector!" );
00377 
00378     // Compute first <r,z> a.k.a. rHz
00379     MVT::MvTransMv( one, *R_, *Z_, rHz );
00380     
00382     // Iterate until the status test tells us to stop.
00383     //
00384     while (stest_->checkStatus(this) != Passed) {
00385       
00386       // Increment the iteration
00387       iter_++;
00388 
00389       // Multiply the current direction vector by A and store in AP_
00390       lp_->applyOp( *P_, *AP_ );
00391       
00392       // Compute alpha := <R_,Z_> / <P_,AP_>
00393       MVT::MvTransMv( one, *P_, *AP_, pAp );
00394       alpha(0,0) = rHz(0,0) / pAp(0,0);
00395       
00396       // Check that alpha is a positive number!
00397       TEST_FOR_EXCEPTION( SCT::real(alpha(0,0)) <= zero, CGIterateFailure,
00398         "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00399       //
00400       // Update the solution vector x := x + alpha * P_
00401       //
00402       MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
00403       lp_->updateSolution();
00404       //
00405       // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
00406       //
00407       rHz_old(0,0) = rHz(0,0);
00408       //
00409       // Compute the new residual R_ := R_ - alpha * AP_
00410       //
00411       MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
00412       //
00413       // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ], 
00414       // and the new direction vector p.
00415       //
00416       if ( lp_->getLeftPrec() != Teuchos::null ) {
00417   lp_->applyLeftPrec( *R_, *Z_ );
00418       } else {
00419   Z_ = R_;
00420       }
00421       //
00422       MVT::MvTransMv( one, *R_, *Z_, rHz );
00423       //
00424       beta(0,0) = rHz(0,0) / rHz_old(0,0);
00425       //
00426       MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
00427       
00428     } // end while (sTest_->checkStatus(this) != Passed)
00429   }
00430 
00431 } // end Belos namespace
00432 
00433 #endif /* BELOS_CG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:05:14 2011 for Belos Package Browser (Single Doxygen Collection) by  doxygen 1.6.3