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