Belos Package Browser (Single Doxygen Collection) Development
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   MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00347       } else {
00348   Z_ = R_;
00349   MVT::MvAddMv( one, *R_, zero, *R_, *P_ );
00350       }
00351     }
00352     else {
00353 
00354       TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00355                          "Belos::CGIter::initialize(): CGIterationState does not have initial residual.");
00356     }
00357 
00358     // The solver is initialized
00359     initialized_ = true;
00360   }
00361 
00362 
00364   // Iterate until the status test informs us we should stop.
00365   template <class ScalarType, class MV, class OP>
00366   void CGIter<ScalarType,MV,OP>::iterate()
00367   {
00368     //
00369     // Allocate/initialize data structures
00370     //
00371     if (initialized_ == false) {
00372       initialize();
00373     }
00374 
00375     // Allocate memory for scalars.
00376     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
00377     Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
00378     Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 ), pAp( 1, 1 );
00379 
00380     // Create convenience variables for zero and one.
00381     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00382     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00383     
00384     // Get the current solution vector.
00385     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00386 
00387     // Check that the current solution vector only has one column. 
00388     TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, CGIterateFailure,
00389                         "Belos::CGIter::iterate(): current linear system has more than one vector!" );
00390 
00391     // Compute first <r,z> a.k.a. rHz
00392     MVT::MvTransMv( one, *R_, *Z_, rHz );
00393     
00395     // Iterate until the status test tells us to stop.
00396     //
00397     while (stest_->checkStatus(this) != Passed) {
00398       
00399       // Increment the iteration
00400       iter_++;
00401 
00402       // Multiply the current direction vector by A and store in AP_
00403       lp_->applyOp( *P_, *AP_ );
00404       
00405       // Compute alpha := <R_,Z_> / <P_,AP_>
00406       MVT::MvTransMv( one, *P_, *AP_, pAp );
00407       alpha(0,0) = rHz(0,0) / pAp(0,0);
00408       
00409       // Check that alpha is a positive number!
00410       TEST_FOR_EXCEPTION( SCT::real(alpha(0,0)) <= zero, CGIterateFailure,
00411         "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00412       //
00413       // Update the solution vector x := x + alpha * P_
00414       //
00415       MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
00416       lp_->updateSolution();
00417       //
00418       // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
00419       //
00420       rHz_old(0,0) = rHz(0,0);
00421       //
00422       // Compute the new residual R_ := R_ - alpha * AP_
00423       //
00424       MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
00425       //
00426       // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ], 
00427       // and the new direction vector p.
00428       //
00429       if ( lp_->getLeftPrec() != Teuchos::null ) {
00430   lp_->applyLeftPrec( *R_, *Z_ );
00431       } else {
00432   Z_ = R_;
00433       }
00434       //
00435       MVT::MvTransMv( one, *R_, *Z_, rHz );
00436       //
00437       beta(0,0) = rHz(0,0) / rHz_old(0,0);
00438       //
00439       MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
00440       
00441     } // end while (sTest_->checkStatus(this) != Passed)
00442   }
00443 
00444 } // end Belos namespace
00445 
00446 #endif /* BELOS_CG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines