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 MultiVecTraitsExt<ScalarType,MV> MVText;
00087   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00088   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00089   typedef typename SCT::magnitudeType MagnitudeType;
00090 
00092 
00093 
00099   CGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00100       const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00101       const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00102       Teuchos::ParameterList &params );
00103 
00105   virtual ~CGIter() {};
00107 
00108 
00110 
00111   
00124   void iterate();
00125 
00140   void initializeCG(CGIterationState<ScalarType,MV> newstate);
00141 
00145   void initialize()
00146   {
00147     CGIterationState<ScalarType,MV> empty;
00148     initializeCG(empty);
00149   }
00150   
00157   CGIterationState<ScalarType,MV> getState() const {
00158     CGIterationState<ScalarType,MV> state;
00159     state.R = R_;
00160     state.P = P_;
00161     state.AP = AP_;
00162     state.Z = Z_;
00163     return state;
00164   }
00165 
00167 
00168   
00170 
00171 
00173   int getNumIters() const { return iter_; }
00174   
00176   void resetNumIters( int iter = 0 ) { iter_ = iter; }
00177 
00180   Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00181 
00183 
00185   Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00186 
00188   
00190 
00191 
00193   const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00194 
00196   int getBlockSize() const { return 1; }
00197 
00199   void setBlockSize(int blockSize) {
00200     TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00201            "Belos::CGIter::setBlockSize(): Cannot use a block size that is not one.");
00202   }
00203 
00205   bool isInitialized() { return initialized_; }
00206 
00208 
00209   private:
00210 
00211   //
00212   // Internal methods
00213   //
00215   void setStateSize();
00216   
00217   //
00218   // Classes inputed through constructor that define the linear problem to be solved.
00219   //
00220   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00221   const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00222   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00223 
00224   //  
00225   // Current solver state
00226   //
00227   // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00228   // is capable of running; _initialize is controlled  by the initialize() member method
00229   // For the implications of the state of initialized_, please see documentation for initialize()
00230   bool initialized_;
00231 
00232   // stateStorageInitialized_ specifies that the state storage has been initialized.
00233   // This initialization may be postponed if the linear problem was generated without 
00234   // the right-hand side or solution vectors.
00235   bool stateStorageInitialized_;
00236 
00237   // Current number of iterations performed.
00238   int iter_;
00239   
00240   // 
00241   // State Storage
00242   // 
00243   // Residual
00244   Teuchos::RCP<MV> R_;
00245   //
00246   // Preconditioned residual
00247   Teuchos::RCP<MV> Z_;
00248   //
00249   // Direction vector
00250   Teuchos::RCP<MV> P_;
00251   //
00252   // Operator applied to direction vector
00253   Teuchos::RCP<MV> AP_;
00254 
00255 };
00256 
00258   // Constructor.
00259   template<class ScalarType, class MV, class OP>
00260   CGIter<ScalarType,MV,OP>::CGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00261                const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00262                const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00263                Teuchos::ParameterList &params ):
00264     lp_(problem),
00265     om_(printer),
00266     stest_(tester),
00267     initialized_(false),
00268     stateStorageInitialized_(false),
00269     iter_(0)
00270   {
00271   }
00272 
00274   // Setup the state storage.
00275   template <class ScalarType, class MV, class OP>
00276   void CGIter<ScalarType,MV,OP>::setStateSize ()
00277   {
00278     if (!stateStorageInitialized_) {
00279 
00280       // Check if there is any multivector to clone from.
00281       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00282       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00283       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00284   stateStorageInitialized_ = false;
00285   return;
00286       }
00287       else {
00288   
00289   // Initialize the state storage
00290   // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00291   if (R_ == Teuchos::null) {
00292     // Get the multivector that is not null.
00293     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00294     TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00295            "Belos::CGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00296     R_ = MVT::Clone( *tmp, 1 );
00297     Z_ = MVT::Clone( *tmp, 1 );
00298     P_ = MVT::Clone( *tmp, 1 );
00299     AP_ = MVT::Clone( *tmp, 1 );
00300   }
00301   
00302   // State storage has now been initialized.
00303   stateStorageInitialized_ = true;
00304       }
00305     }
00306   }
00307 
00308 
00310   // Initialize this iteration object
00311   template <class ScalarType, class MV, class OP>
00312   void CGIter<ScalarType,MV,OP>::initializeCG(CGIterationState<ScalarType,MV> newstate)
00313   {
00314     // Initialize the state storage if it isn't already.
00315     if (!stateStorageInitialized_) 
00316       setStateSize();
00317 
00318     TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00319            "Belos::CGIter::initialize(): Cannot initialize state storage!");
00320     
00321     // NOTE:  In CGIter R_, the initial residual, is required!!!  
00322     //
00323     std::string errstr("Belos::CGIter::initialize(): Specified multivectors must have a consistent length and width.");
00324 
00325     // Create convenience variables for zero and one.
00326     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00327     const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
00328 
00329     if (newstate.R != Teuchos::null) {
00330 
00331       TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*newstate.R) != MVText::GetGlobalLength(*R_),
00332                           std::invalid_argument, errstr );
00333       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != 1,
00334                           std::invalid_argument, errstr );
00335 
00336       // Copy basis vectors from newstate into V
00337       if (newstate.R != R_) {
00338         // copy over the initial residual (unpreconditioned).
00339   MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
00340       }
00341 
00342       // Compute initial direction vectors
00343       // Initially, they are set to the preconditioned residuals
00344       //
00345       if ( lp_->getLeftPrec() != Teuchos::null ) {
00346         lp_->applyLeftPrec( *R_, *Z_ );
00347         if ( lp_->getRightPrec() != Teuchos::null ) {
00348           Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
00349           lp_->applyRightPrec( *Z_, *tmp );
00350           Z_ = tmp;
00351         }
00352       }
00353       else if ( lp_->getRightPrec() != Teuchos::null ) {
00354         lp_->applyRightPrec( *R_, *Z_ );
00355       } 
00356       else {
00357         Z_ = R_;
00358       }
00359       MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00360     }
00361     else {
00362 
00363       TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00364                          "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
00365     }
00366 
00367     // The solver is initialized
00368     initialized_ = true;
00369   }
00370 
00371 
00373   // Iterate until the status test informs us we should stop.
00374   template <class ScalarType, class MV, class OP>
00375   void CGIter<ScalarType,MV,OP>::iterate()
00376   {
00377     //
00378     // Allocate/initialize data structures
00379     //
00380     if (initialized_ == false) {
00381       initialize();
00382     }
00383 
00384     // Allocate memory for scalars.
00385     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
00386     Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
00387     Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 ), pAp( 1, 1 );
00388 
00389     // Create convenience variables for zero and one.
00390     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00391     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00392     
00393     // Get the current solution vector.
00394     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00395 
00396     // Check that the current solution vector only has one column. 
00397     TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
00398                         "Belos::CGIter::iterate(): current linear system has more than one vector!" );
00399 
00400     // Compute first <r,z> a.k.a. rHz
00401     MVT::MvTransMv( one, *R_, *Z_, rHz );
00402     
00404     // Iterate until the status test tells us to stop.
00405     //
00406     while (stest_->checkStatus(this) != Passed) {
00407       
00408       // Increment the iteration
00409       iter_++;
00410 
00411       // Multiply the current direction vector by A and store in AP_
00412       lp_->applyOp( *P_, *AP_ );
00413       
00414       // Compute alpha := <R_,Z_> / <P_,AP_>
00415       MVT::MvTransMv( one, *P_, *AP_, pAp );
00416       alpha(0,0) = rHz(0,0) / pAp(0,0);
00417       
00418       // Check that alpha is a positive number!
00419       TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(alpha(0,0)) <= zero, CGIterateFailure,
00420         "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00421       //
00422       // Update the solution vector x := x + alpha * P_
00423       //
00424       MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
00425       lp_->updateSolution();
00426       //
00427       // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
00428       //
00429       rHz_old(0,0) = rHz(0,0);
00430       //
00431       // Compute the new residual R_ := R_ - alpha * AP_
00432       //
00433       MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
00434       //
00435       // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ], 
00436       // and the new direction vector p.
00437       //
00438       if ( lp_->getLeftPrec() != Teuchos::null ) {
00439         lp_->applyLeftPrec( *R_, *Z_ );
00440         if ( lp_->getRightPrec() != Teuchos::null ) {
00441           Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, 1 );
00442           lp_->applyRightPrec( *Z_, *tmp );
00443           Z_ = tmp;
00444         }
00445       }
00446       else if ( lp_->getRightPrec() != Teuchos::null ) {
00447         lp_->applyRightPrec( *R_, *Z_ );
00448       } 
00449       else {
00450         Z_ = R_;
00451       }
00452       //
00453       MVT::MvTransMv( one, *R_, *Z_, rHz );
00454       //
00455       beta(0,0) = rHz(0,0) / rHz_old(0,0);
00456       //
00457       MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
00458       
00459     } // end while (sTest_->checkStatus(this) != Passed)
00460   }
00461 
00462 } // end Belos namespace
00463 
00464 #endif /* BELOS_CG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines