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