Belos Version of the Day
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 the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
00008 // the U.S. Government retains certain rights in this software.
00009 //
00010 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Michael A. Heroux (maherou@sandia.gov)
00038 //
00039 // ************************************************************************
00040 //@HEADER
00041 
00042 #ifndef BELOS_CG_ITER_HPP
00043 #define BELOS_CG_ITER_HPP
00044 
00049 #include "BelosConfigDefs.hpp"
00050 #include "BelosTypes.hpp"
00051 #include "BelosCGIteration.hpp"
00052 
00053 #include "BelosLinearProblem.hpp"
00054 #include "BelosOutputManager.hpp"
00055 #include "BelosStatusTest.hpp"
00056 #include "BelosOperatorTraits.hpp"
00057 #include "BelosMultiVecTraits.hpp"
00058 
00059 #include "Teuchos_SerialDenseMatrix.hpp"
00060 #include "Teuchos_SerialDenseVector.hpp"
00061 #include "Teuchos_ScalarTraits.hpp"
00062 #include "Teuchos_ParameterList.hpp"
00063 #include "Teuchos_TimeMonitor.hpp"
00064 
00075 namespace Belos {
00076   
00077 template<class ScalarType, class MV, class OP>
00078 class CGIter : virtual public CGIteration<ScalarType,MV,OP> {
00079 
00080   public:
00081     
00082   //
00083   // Convenience typedefs
00084   //
00085   typedef MultiVecTraits<ScalarType,MV> MVT;
00086   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00087   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00088   typedef typename SCT::magnitudeType MagnitudeType;
00089 
00091 
00092 
00098   CGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00099       const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00100       const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00101       Teuchos::ParameterList &params );
00102 
00104   virtual ~CGIter() {};
00106 
00107 
00109 
00110   
00123   void iterate();
00124 
00139   void initializeCG(CGIterationState<ScalarType,MV> newstate);
00140 
00144   void initialize()
00145   {
00146     CGIterationState<ScalarType,MV> empty;
00147     initializeCG(empty);
00148   }
00149   
00156   CGIterationState<ScalarType,MV> getState() const {
00157     CGIterationState<ScalarType,MV> state;
00158     state.R = R_;
00159     state.P = P_;
00160     state.AP = AP_;
00161     state.Z = Z_;
00162     return state;
00163   }
00164 
00166 
00167   
00169 
00170 
00172   int getNumIters() const { return iter_; }
00173   
00175   void resetNumIters( int iter = 0 ) { iter_ = iter; }
00176 
00179   Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00180 
00182 
00184   Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00185 
00187   
00189 
00190 
00192   const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00193 
00195   int getBlockSize() const { return 1; }
00196 
00198   void setBlockSize(int blockSize) {
00199     TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00200            "Belos::CGIter::setBlockSize(): Cannot use a block size that is not one.");
00201   }
00202 
00204   bool isInitialized() { return initialized_; }
00205 
00207 
00208   private:
00209 
00210   //
00211   // Internal methods
00212   //
00214   void setStateSize();
00215   
00216   //
00217   // Classes inputed through constructor that define the linear problem to be solved.
00218   //
00219   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00220   const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00221   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00222 
00223   //  
00224   // Current solver state
00225   //
00226   // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00227   // is capable of running; _initialize is controlled  by the initialize() member method
00228   // For the implications of the state of initialized_, please see documentation for initialize()
00229   bool initialized_;
00230 
00231   // stateStorageInitialized_ specifies that the state storage has been initialized.
00232   // This initialization may be postponed if the linear problem was generated without 
00233   // the right-hand side or solution vectors.
00234   bool stateStorageInitialized_;
00235 
00236   // Current number of iterations performed.
00237   int iter_;
00238   
00239   // 
00240   // State Storage
00241   // 
00242   // Residual
00243   Teuchos::RCP<MV> R_;
00244   //
00245   // Preconditioned residual
00246   Teuchos::RCP<MV> Z_;
00247   //
00248   // Direction vector
00249   Teuchos::RCP<MV> P_;
00250   //
00251   // Operator applied to direction vector
00252   Teuchos::RCP<MV> AP_;
00253 
00254 };
00255 
00257   // Constructor.
00258   template<class ScalarType, class MV, class OP>
00259   CGIter<ScalarType,MV,OP>::CGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00260                const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00261                const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00262                Teuchos::ParameterList &params ):
00263     lp_(problem),
00264     om_(printer),
00265     stest_(tester),
00266     initialized_(false),
00267     stateStorageInitialized_(false),
00268     iter_(0)
00269   {
00270   }
00271 
00273   // Setup the state storage.
00274   template <class ScalarType, class MV, class OP>
00275   void CGIter<ScalarType,MV,OP>::setStateSize ()
00276   {
00277     if (!stateStorageInitialized_) {
00278 
00279       // Check if there is any multivector to clone from.
00280       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00281       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00282       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00283   stateStorageInitialized_ = false;
00284   return;
00285       }
00286       else {
00287   
00288   // Initialize the state storage
00289   // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00290   if (R_ == Teuchos::null) {
00291     // Get the multivector that is not null.
00292     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00293     TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00294            "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00295     R_ = MVT::Clone( *tmp, 1 );
00296     Z_ = MVT::Clone( *tmp, 1 );
00297     P_ = MVT::Clone( *tmp, 1 );
00298     AP_ = MVT::Clone( *tmp, 1 );
00299   }
00300   
00301   // State storage has now been initialized.
00302   stateStorageInitialized_ = true;
00303       }
00304     }
00305   }
00306 
00307 
00309   // Initialize this iteration object
00310   template <class ScalarType, class MV, class OP>
00311   void CGIter<ScalarType,MV,OP>::initializeCG(CGIterationState<ScalarType,MV> newstate)
00312   {
00313     // Initialize the state storage if it isn't already.
00314     if (!stateStorageInitialized_) 
00315       setStateSize();
00316 
00317     TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00318            "Belos::CGIter::initialize(): Cannot initialize state storage!");
00319     
00320     // NOTE:  In CGIter R_, the initial residual, is required!!!  
00321     //
00322     std::string errstr("Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
00323 
00324     // Create convenience variables for zero and one.
00325     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00326     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00327 
00328     if (newstate.R != Teuchos::null) {
00329 
00330       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
00331                           std::invalid_argument, errstr );
00332       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
00333                           std::invalid_argument, errstr );
00334 
00335       // Copy basis vectors from newstate into V
00336       if (newstate.R != R_) {
00337         // copy over the initial residual (unpreconditioned).
00338   MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
00339       }
00340 
00341       // Compute initial direction vectors
00342       // Initially, they are set to the preconditioned residuals
00343       //
00344       if ( lp_->getLeftPrec() != Teuchos::null ) {
00345         lp_->applyLeftPrec( *R_, *Z_ );
00346         if ( lp_->getRightPrec() != Teuchos::null ) {
00347           Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
00348           lp_->applyRightPrec( *Z_, *tmp );
00349           Z_ = tmp;
00350         }
00351       }
00352       else if ( lp_->getRightPrec() != Teuchos::null ) {
00353         lp_->applyRightPrec( *R_, *Z_ );
00354       } 
00355       else {
00356         Z_ = R_;
00357       }
00358       MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00359     }
00360     else {
00361 
00362       TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00363                          "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
00364     }
00365 
00366     // The solver is initialized
00367     initialized_ = true;
00368   }
00369 
00370 
00372   // Iterate until the status test informs us we should stop.
00373   template <class ScalarType, class MV, class OP>
00374   void CGIter<ScalarType,MV,OP>::iterate()
00375   {
00376     //
00377     // Allocate/initialize data structures
00378     //
00379     if (initialized_ == false) {
00380       initialize();
00381     }
00382 
00383     // Allocate memory for scalars.
00384     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
00385     Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
00386     Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 ), pAp( 1, 1 );
00387 
00388     // Create convenience variables for zero and one.
00389     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00390     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00391     
00392     // Get the current solution vector.
00393     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00394 
00395     // Check that the current solution vector only has one column. 
00396     TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
00397                         "Belos::CGIter::iterate(): current linear system has more than one vector!" );
00398 
00399     // Compute first <r,z> a.k.a. rHz
00400     MVT::MvTransMv( one, *R_, *Z_, rHz );
00401     
00403     // Iterate until the status test tells us to stop.
00404     //
00405     while (stest_->checkStatus(this) != Passed) {
00406       
00407       // Increment the iteration
00408       iter_++;
00409 
00410       // Multiply the current direction vector by A and store in AP_
00411       lp_->applyOp( *P_, *AP_ );
00412       
00413       // Compute alpha := <R_,Z_> / <P_,AP_>
00414       MVT::MvTransMv( one, *P_, *AP_, pAp );
00415       alpha(0,0) = rHz(0,0) / pAp(0,0);
00416       
00417       // Check that alpha is a positive number!
00418       TEST_FOR_EXCEPTION( SCT::real(alpha(0,0)) <= zero, CGIterateFailure,
00419         "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00420       //
00421       // Update the solution vector x := x + alpha * P_
00422       //
00423       MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
00424       lp_->updateSolution();
00425       //
00426       // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
00427       //
00428       rHz_old(0,0) = rHz(0,0);
00429       //
00430       // Compute the new residual R_ := R_ - alpha * AP_
00431       //
00432       MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
00433       //
00434       // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ], 
00435       // and the new direction vector p.
00436       //
00437       if ( lp_->getLeftPrec() != Teuchos::null ) {
00438         lp_->applyLeftPrec( *R_, *Z_ );
00439         if ( lp_->getRightPrec() != Teuchos::null ) {
00440           Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
00441           lp_->applyRightPrec( *Z_, *tmp );
00442           Z_ = tmp;
00443         }
00444       }
00445       else if ( lp_->getRightPrec() != Teuchos::null ) {
00446         lp_->applyRightPrec( *R_, *Z_ );
00447       } 
00448       else {
00449         Z_ = R_;
00450       }
00451       //
00452       MVT::MvTransMv( one, *R_, *Z_, rHz );
00453       //
00454       beta(0,0) = rHz(0,0) / rHz_old(0,0);
00455       //
00456       MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
00457       
00458     } // end while (sTest_->checkStatus(this) != Passed)
00459   }
00460 
00461 } // end Belos namespace
00462 
00463 #endif /* BELOS_CG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines