BelosBlockCGIter.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_BLOCK_CG_ITER_HPP
00030 #define BELOS_BLOCK_CG_ITER_HPP
00031 
00036 #include "BelosConfigDefs.hpp"
00037 #include "BelosTypes.hpp"
00038 #include "BelosCGIteration.hpp"
00039 
00040 #include "BelosLinearProblem.hpp"
00041 #include "BelosMatOrthoManager.hpp"
00042 #include "BelosOutputManager.hpp"
00043 #include "BelosStatusTest.hpp"
00044 #include "BelosOperatorTraits.hpp"
00045 #include "BelosMultiVecTraits.hpp"
00046 
00047 #include "Teuchos_BLAS.hpp"
00048 #include "Teuchos_LAPACK.hpp"
00049 #include "Teuchos_SerialDenseMatrix.hpp"
00050 #include "Teuchos_SerialDenseVector.hpp"
00051 #include "Teuchos_ScalarTraits.hpp"
00052 #include "Teuchos_ParameterList.hpp"
00053 #include "Teuchos_TimeMonitor.hpp"
00054 
00063 namespace Belos {
00064   
00065 template<class ScalarType, class MV, class OP>
00066 class BlockCGIter : virtual public CGIteration<ScalarType,MV,OP> {
00067 
00068   public:
00069     
00070   //
00071   // Convenience typedefs
00072   //
00073   typedef MultiVecTraits<ScalarType,MV> MVT;
00074   typedef OperatorTraits<ScalarType,MV,OP> OPT;
00075   typedef Teuchos::ScalarTraits<ScalarType> SCT;
00076   typedef typename SCT::magnitudeType MagnitudeType;
00077 
00079 
00080 
00086   BlockCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00087          const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00088          const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00089          const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00090          Teuchos::ParameterList &params );
00091 
00093   virtual ~BlockCGIter() {};
00095 
00096 
00098 
00099   
00112   void iterate();
00113 
00128   void initializeCG(CGIterationState<ScalarType,MV> newstate);
00129 
00133   void initialize()
00134   {
00135     CGIterationState<ScalarType,MV> empty;
00136     initializeCG(empty);
00137   }
00138   
00145   CGIterationState<ScalarType,MV> getState() const {
00146     CGIterationState<ScalarType,MV> state;
00147     state.R = R_;
00148     state.P = P_;
00149     state.AP = AP_;
00150     state.Z = Z_;
00151     return state;
00152   }
00153 
00155 
00156   
00158 
00159 
00161   int getNumIters() const { return iter_; }
00162   
00164   void resetNumIters( int iter=0 ) { iter_ = iter; }
00165 
00168   Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00169 
00171 
00173   Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00174 
00176   
00178 
00179 
00181   const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00182 
00184   int getBlockSize() const { return blockSize_; }
00185 
00187   void setBlockSize(int blockSize);
00188 
00190   bool isInitialized() { return initialized_; }
00191 
00193 
00194   private:
00195 
00196   //
00197   // Internal methods
00198   //
00200   void setStateSize();
00201   
00202   //
00203   // Classes inputed through constructor that define the linear problem to be solved.
00204   //
00205   const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00206   const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00207   const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00208   const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;
00209 
00210   //
00211   // Algorithmic parameters
00212   //
00213   // blockSize_ is the solver block size.
00214   int blockSize_;
00215 
00216   //  
00217   // Current solver state
00218   //
00219   // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00220   // is capable of running; _initialize is controlled  by the initialize() member method
00221   // For the implications of the state of initialized_, please see documentation for initialize()
00222   bool initialized_;
00223 
00224   // stateStorageInitialized_ specified that the state storage has be initialized.
00225   // This initialization may be postponed if the linear problem was generated without 
00226   // the right-hand side or solution vectors.
00227   bool stateStorageInitialized_;
00228 
00229   // Current subspace dimension, and number of iterations performed.
00230   int iter_;
00231   
00232   // 
00233   // State Storage
00234   // 
00235   // Residual
00236   Teuchos::RCP<MV> R_;
00237   //
00238   // Preconditioned residual
00239   Teuchos::RCP<MV> Z_;
00240   //
00241   // Direction std::vector
00242   Teuchos::RCP<MV> P_;
00243   //
00244   // Operator applied to direction std::vector
00245   Teuchos::RCP<MV> AP_;
00246 
00247 };
00248 
00250   // Constructor.
00251   template<class ScalarType, class MV, class OP>
00252   BlockCGIter<ScalarType,MV,OP>::BlockCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00253                const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00254                const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00255                const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
00256                Teuchos::ParameterList &params ):
00257     lp_(problem),
00258     om_(printer),
00259     stest_(tester),
00260     ortho_(ortho),
00261     blockSize_(0),
00262     initialized_(false),
00263     stateStorageInitialized_(false),
00264     iter_(0)
00265   {
00266     // Set the block size and allocate data
00267     int bs = params.get("Block Size", 1);
00268     setBlockSize( bs );
00269   }
00270 
00272   // Setup the state storage.
00273   template <class ScalarType, class MV, class OP>
00274   void BlockCGIter<ScalarType,MV,OP>::setStateSize ()
00275   {
00276     if (!stateStorageInitialized_) {
00277 
00278       // Check if there is any multivector to clone from.
00279       Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
00280       Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
00281       if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
00282   stateStorageInitialized_ = false;
00283   return;
00284       }
00285       else {
00286   
00287   // Initialize the state storage
00288   // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
00289   if (R_ == Teuchos::null || MVT::GetNumberVecs(*R_)!=blockSize_) {
00290     // Get the multivector that is not null.
00291     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00292     TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
00293            "Belos::BlockCGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
00294     R_ = MVT::Clone( *tmp, blockSize_ );
00295     Z_ = MVT::Clone( *tmp, blockSize_ );
00296     P_ = MVT::Clone( *tmp, blockSize_ );
00297     AP_ = MVT::Clone( *tmp, blockSize_ );
00298   }
00299   
00300   // State storage has now been initialized.
00301   stateStorageInitialized_ = true;
00302       }
00303     }
00304   }
00305 
00307   // Set the block size and make necessary adjustments.
00308   template <class ScalarType, class MV, class OP>
00309   void BlockCGIter<ScalarType,MV,OP>::setBlockSize (int blockSize)
00310   {
00311     // This routine only allocates space; it doesn't not perform any computation
00312     // any change in size will invalidate the state of the solver.
00313 
00314     TEST_FOR_EXCEPTION(blockSize <= 0, std::invalid_argument, "Belos::BlockGmresIter::setBlockSize was passed a non-positive argument.");
00315     if (blockSize == blockSize_) {
00316       // do nothing
00317       return;
00318     }
00319 
00320     if (blockSize!=blockSize_)
00321       stateStorageInitialized_ = false;
00322 
00323     blockSize_ = blockSize;
00324 
00325     initialized_ = false;
00326 
00327     // Use the current blockSize_ to initialize the state storage.
00328     setStateSize();
00329 
00330   }
00331 
00333   // Initialize this iteration object
00334   template <class ScalarType, class MV, class OP>
00335   void BlockCGIter<ScalarType,MV,OP>::initializeCG(CGIterationState<ScalarType,MV> newstate)
00336   {
00337     // Initialize the state storage if it isn't already.
00338     if (!stateStorageInitialized_) 
00339       setStateSize();
00340 
00341     TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
00342            "Belos::BlockCGIter::initialize(): Cannot initialize state storage!");
00343     
00344     // NOTE:  In BlockCGIter R_, the initial residual, is required!!!  
00345     //
00346     std::string errstr("Belos::BlockCGIter::initialize(): Specified multivectors must have a consistent length and width.");
00347 
00348     // Create convenience variables for zero and one.
00349     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00350     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00351 
00352     if (newstate.R != Teuchos::null) {
00353 
00354       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
00355                           std::invalid_argument, errstr );
00356       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != blockSize_,
00357                           std::invalid_argument, errstr );
00358 
00359       // Copy basis vectors from newstate into V
00360       if (newstate.R != R_) {
00361         // copy over the initial residual (unpreconditioned).
00362   MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
00363       }
00364       // Compute initial direction vectors
00365       // Initially, they are set to the preconditioned residuals
00366       //
00367       if ( lp_->getLeftPrec() != Teuchos::null ) {
00368   lp_->applyLeftPrec( *R_, *Z_ ); 
00369   MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00370       } else {
00371   Z_ = R_;
00372   MVT::MvAddMv( one, *R_, zero, *R_, *P_ );
00373       }
00374     }
00375     else {
00376 
00377       TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
00378                          "Belos::BlockCGIter::initialize(): BlockCGStateIterState does not have initial residual.");
00379     }
00380 
00381     // The solver is initialized
00382     initialized_ = true;
00383   }
00384 
00385 
00387   // Iterate until the status test informs us we should stop.
00388   template <class ScalarType, class MV, class OP>
00389   void BlockCGIter<ScalarType,MV,OP>::iterate()
00390   {
00391     //
00392     // Allocate/initialize data structures
00393     //
00394     if (initialized_ == false) {
00395       initialize();
00396     }
00397     // Allocate data needed for LAPACK work.
00398     int info = 0;
00399     char UPLO = 'U';
00400     Teuchos::LAPACK<int,ScalarType> lapack;
00401 
00402     // Allocate memory for scalars.
00403     Teuchos::SerialDenseMatrix<int,ScalarType> alpha( blockSize_, blockSize_ );
00404     Teuchos::SerialDenseMatrix<int,ScalarType> beta( blockSize_, blockSize_ );
00405     Teuchos::SerialDenseMatrix<int,ScalarType> rHz( blockSize_, blockSize_ ), 
00406       rHz_old( blockSize_, blockSize_ ), pAp( blockSize_, blockSize_ );
00407 
00408     // Create convenience variables for zero and one.
00409     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00410     
00411     // Get the current solution std::vector.
00412     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00413 
00414     // Check that the current solution std::vector has blockSize_ columns. 
00415     TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != blockSize_, CGIterateFailure,
00416                         "Belos::BlockCGIter::iterate(): current linear system does not have the right number of vectors!" );
00417     int rank = ortho_->normalize( *P_, Teuchos::null );
00418     TEST_FOR_EXCEPTION(rank != blockSize_,CGIterationOrthoFailure,
00419                          "Belos::BlockCGIter::iterate(): Failed to compute initial block of orthonormal direction vectors.");
00420 
00421 
00423     // Iterate until the status test tells us to stop.
00424     //
00425     while (stest_->checkStatus(this) != Passed) {
00426         
00427       // Increment the iteration
00428       iter_++;
00429     
00430       // Multiply the current direction std::vector by A and store in Ap_
00431       lp_->applyOp( *P_, *AP_ );
00432       
00433       // Compute alpha := <P_,R_> / <P_,AP_>
00434       // 1) Compute P^T * A * P = pAp and P^T * R 
00435       // 2) Compute the Cholesky Factorization of pAp
00436       // 3) Back and forward solves to compute alpha
00437       //
00438       MVT::MvTransMv( one, *P_, *R_, alpha );
00439       MVT::MvTransMv( one, *P_, *AP_, pAp );      
00440      
00441       // Compute Cholesky factorization of pAp
00442       lapack.POTRF(UPLO, blockSize_, pAp.values(), blockSize_, &info);
00443       TEST_FOR_EXCEPTION(info != 0,CGIterationLAPACKFailure,
00444                          "Belos::BlockCGIter::iterate(): Failed to compute Cholesky factorization using LAPACK routine POTRF.");
00445 
00446       // Compute alpha by performing a back and forward solve with the Cholesky factorization in pAp.
00447       lapack.POTRS(UPLO, blockSize_, blockSize_, pAp.values(), blockSize_, 
00448        alpha.values(), blockSize_, &info);
00449       TEST_FOR_EXCEPTION(info != 0,CGIterationLAPACKFailure,
00450                          "Belos::BlockCGIter::iterate(): Failed to compute alpha using Cholesky factorization (POTRS).");
00451       
00452       //
00453       // Update the solution std::vector X := X + alpha * P_
00454       //
00455       MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
00456       lp_->updateSolution();
00457       //
00458       // Compute the new residual R_ := R_ - alpha * AP_
00459       //
00460       MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
00461       //
00462       // Compute the new preconditioned residual, Z_.
00463       if ( lp_->getLeftPrec() != Teuchos::null ) {
00464   lp_->applyLeftPrec( *R_, *Z_ );
00465       } else {
00466   Z_ = R_;
00467       }
00468       //
00469       // Compute beta := <AP_,Z_> / <P_,AP_> 
00470       // 1) Compute AP_^T * Z_ 
00471       // 2) Compute the Cholesky Factorization of pAp (already have)
00472       // 3) Back and forward solves to compute beta
00473 
00474       // Compute <AP_,Z>
00475       MVT::MvTransMv( -one, *AP_, *Z_, beta );
00476       //
00477       lapack.POTRS(UPLO, blockSize_, blockSize_, pAp.values(), blockSize_, 
00478        beta.values(), blockSize_, &info);
00479       TEST_FOR_EXCEPTION(info != 0,CGIterationLAPACKFailure,
00480                          "Belos::BlockCGIter::iterate(): Failed to compute beta using Cholesky factorization (POTRS).");
00481       //
00482       // Compute the new direction vectors P_ = Z_ + P_ * beta 
00483       //
00484       Teuchos::RCP<MV> Pnew = MVT::CloneCopy( *Z_ );
00485       MVT::MvTimesMatAddMv(one, *P_, beta, one, *Pnew);
00486       P_ = Pnew;
00487 
00488       // Compute orthonormal block of new direction vectors.
00489       rank = ortho_->normalize( *P_, Teuchos::null );
00490       TEST_FOR_EXCEPTION(rank != blockSize_,CGIterationOrthoFailure,
00491                          "Belos::BlockCGIter::iterate(): Failed to compute block of orthonormal direction vectors.");
00492       
00493     } // end while (sTest_->checkStatus(this) != Passed)
00494   }
00495 
00496 } // end Belos namespace
00497 
00498 #endif /* BELOS_BLOCK_CG_ITER_HPP */

Generated on Wed May 12 21:45:50 2010 for Belos by  doxygen 1.4.7