Belos Version of the Day
BelosPseudoBlockCGIter.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_PSEUDO_BLOCK_CG_ITER_HPP
00043 #define BELOS_PSEUDO_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_SerialDenseMatrix.hpp"
00062 #include "Teuchos_SerialDenseVector.hpp"
00063 #include "Teuchos_ScalarTraits.hpp"
00064 #include "Teuchos_ParameterList.hpp"
00065 #include "Teuchos_TimeMonitor.hpp"
00066 
00078 namespace Belos {
00079   
00080   template<class ScalarType, class MV, class OP>
00081   class PseudoBlockCGIter : virtual public CGIteration<ScalarType,MV,OP> {
00082     
00083   public:
00084     
00085     //
00086     // Convenience typedefs
00087     //
00088     typedef MultiVecTraits<ScalarType,MV> MVT;
00089     typedef OperatorTraits<ScalarType,MV,OP> OPT;
00090     typedef Teuchos::ScalarTraits<ScalarType> SCT;
00091     typedef typename SCT::magnitudeType MagnitudeType;
00092     
00094 
00095     
00101     PseudoBlockCGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00102         const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00103         const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00104         Teuchos::ParameterList &params );
00105     
00107     virtual ~PseudoBlockCGIter() {};
00109     
00110     
00112 
00113     
00127     void iterate();
00128     
00149     void initializeCG(CGIterationState<ScalarType,MV> newstate);
00150     
00154     void initialize()
00155     {
00156       CGIterationState<ScalarType,MV> empty;
00157       initializeCG(empty);
00158     }
00159     
00167     CGIterationState<ScalarType,MV> getState() const {
00168       CGIterationState<ScalarType,MV> state;
00169       state.R = R_;
00170       state.P = P_;
00171       state.AP = AP_;
00172       state.Z = Z_;
00173       return state;
00174     }
00175     
00177     
00178     
00180 
00181     
00183     int getNumIters() const { return iter_; }
00184     
00186     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00187     
00190     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00191     
00193 
00195     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00196     
00198     
00200 
00201     
00203     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00204     
00206     int getBlockSize() const { return 1; }
00207     
00209     void setBlockSize(int blockSize) { 
00210       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00211        "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
00212     }
00213     
00215     bool isInitialized() { return initialized_; }
00216     
00218     
00219   private:
00220     
00221     //
00222     // Classes inputed through constructor that define the linear problem to be solved.
00223     //
00224     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00225     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00226     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00227     
00228     //
00229     // Algorithmic parameters
00230     //  
00231     // numRHS_ is the current number of linear systems being solved.
00232     int numRHS_;
00233     
00234     // 
00235     // Current solver state
00236     //
00237     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00238     // is capable of running; _initialize is controlled  by the initialize() member method
00239     // For the implications of the state of initialized_, please see documentation for initialize()
00240     bool initialized_;
00241     
00242     // Current number of iterations performed.
00243     int iter_;
00244 
00245     // Current number of iterations performed.
00246     bool assertPositiveDefiniteness_;
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 vector
00258     Teuchos::RCP<MV> P_;
00259     //
00260     // Operator applied to direction vector
00261     Teuchos::RCP<MV> AP_;
00262 
00263   };
00264   
00266   // Constructor.
00267   template<class ScalarType, class MV, class OP>
00268   PseudoBlockCGIter<ScalarType,MV,OP>::PseudoBlockCGIter(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                      Teuchos::ParameterList &params ):
00272     lp_(problem),
00273     om_(printer),
00274     stest_(tester),
00275     numRHS_(0),
00276     initialized_(false),
00277     iter_(0),
00278     assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) )
00279   {
00280   }
00281   
00282 
00284   // Initialize this iteration object
00285   template <class ScalarType, class MV, class OP>
00286   void PseudoBlockCGIter<ScalarType,MV,OP>::initializeCG(CGIterationState<ScalarType,MV> newstate)
00287   {
00288     // Check if there is any multivector to clone from.
00289     Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
00290     Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
00291     TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
00292            "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
00293 
00294     // Get the multivector that is not null.
00295     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00296 
00297     // Get the number of right-hand sides we're solving for now.
00298     int numRHS = MVT::GetNumberVecs(*tmp);
00299     numRHS_ = numRHS;
00300   
00301     // Initialize the state storage
00302     // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
00303     if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
00304       R_ = MVT::Clone( *tmp, numRHS_ );
00305       Z_ = MVT::Clone( *tmp, numRHS_ );
00306       P_ = MVT::Clone( *tmp, numRHS_ );
00307       AP_ = MVT::Clone( *tmp, numRHS_ );
00308     }
00309   
00310     // NOTE:  In CGIter R_, the initial residual, is required!!!  
00311     //
00312     std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
00313 
00314     // Create convenience variables for zero and one.
00315     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00316     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00317 
00318     if (!Teuchos::is_null(newstate.R)) {
00319 
00320       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
00321                           std::invalid_argument, errstr );
00322       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
00323                           std::invalid_argument, errstr );
00324 
00325       // Copy basis vectors from newstate into V
00326       if (newstate.R != R_) {
00327         // copy over the initial residual (unpreconditioned).
00328   MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
00329       }
00330 
00331       // Compute initial direction vectors
00332       // Initially, they are set to the preconditioned residuals
00333       //
00334       if ( lp_->getLeftPrec() != Teuchos::null ) {
00335         lp_->applyLeftPrec( *R_, *Z_ );
00336         if ( lp_->getRightPrec() != Teuchos::null ) {
00337           Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
00338           lp_->applyRightPrec( *Z_, *tmp );
00339           Z_ = tmp;
00340         }
00341       }
00342       else if ( lp_->getRightPrec() != Teuchos::null ) {
00343   lp_->applyRightPrec( *R_, *Z_ ); 
00344       } 
00345       else {
00346   Z_ = R_;
00347       }
00348       MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00349     }
00350     else {
00351 
00352       TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
00353                          "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
00354     }
00355 
00356     // The solver is initialized
00357     initialized_ = true;
00358   }
00359 
00360 
00362   // Iterate until the status test informs us we should stop.
00363   template <class ScalarType, class MV, class OP>
00364   void PseudoBlockCGIter<ScalarType,MV,OP>::iterate()
00365   {
00366     //
00367     // Allocate/initialize data structures
00368     //
00369     if (initialized_ == false) {
00370       initialize();
00371     }
00372 
00373     // Allocate memory for scalars.
00374     int i=0;
00375     std::vector<int> index(1);
00376     std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
00377     Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ );
00378 
00379     // Create convenience variables for zero and one.
00380     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00381     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00382     
00383     // Get the current solution std::vector.
00384     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00385 
00386     // Compute first <r,z> a.k.a. rHz
00387     MVT::MvDot( *R_, *Z_, rHz );
00388 
00389     if ( assertPositiveDefiniteness_ )
00390         for (i=0; i<numRHS_; ++i)
00391             TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
00392                                 CGIterateFailure,
00393                                 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
00394 
00396     // Iterate until the status test tells us to stop.
00397     //
00398     while (stest_->checkStatus(this) != Passed) {
00399       
00400       // Increment the iteration
00401       iter_++;
00402 
00403       // Multiply the current direction std::vector by A and store in AP_
00404       lp_->applyOp( *P_, *AP_ );
00405       
00406       // Compute alpha := <R_,Z_> / <P_,AP_>
00407       MVT::MvDot( *P_, *AP_, pAp );
00408 
00409       for (i=0; i<numRHS_; ++i) {
00410         if ( assertPositiveDefiniteness_ )
00411             // Check that pAp[i] is a positive number!
00412             TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
00413                                 CGIterateFailure,
00414                                 "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00415 
00416         alpha(i,i) = rHz[i] / pAp[i];
00417       }
00418 
00419       //
00420       // Update the solution std::vector x := x + alpha * P_
00421       //
00422       MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
00423       lp_->updateSolution();
00424       //
00425       // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
00426       //
00427       for (i=0; i<numRHS_; ++i) {
00428         rHz_old[i] = rHz[i];
00429       }
00430       //
00431       // Compute the new residual R_ := R_ - alpha * AP_
00432       //
00433       MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
00434       //
00435       // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ], 
00436       // and the new direction std::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_, numRHS_ );
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::MvDot( *R_, *Z_, rHz );
00454       if ( assertPositiveDefiniteness_ )
00455           for (i=0; i<numRHS_; ++i)
00456               TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
00457                                   CGIterateFailure,
00458                                   "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
00459       //
00460       // Update the search directions.
00461       for (i=0; i<numRHS_; ++i) {
00462         beta(i,i) = rHz[i] / rHz_old[i];
00463   index[0] = i;
00464   Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
00465   Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
00466         MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
00467       }
00468       //      
00469     } // end while (sTest_->checkStatus(this) != Passed)
00470   }
00471 
00472 } // end Belos namespace
00473 
00474 #endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines