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

Generated on Tue Oct 20 12:48:34 2009 for Belos by doxygen 1.4.7