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 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     PseudoBlockCGIter( 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         Teuchos::ParameterList &params );
00106     
00108     virtual ~PseudoBlockCGIter() {};
00110     
00111     
00113 
00114     
00128     void iterate();
00129     
00150     void initializeCG(CGIterationState<ScalarType,MV> newstate);
00151     
00155     void initialize()
00156     {
00157       CGIterationState<ScalarType,MV> empty;
00158       initializeCG(empty);
00159     }
00160     
00168     CGIterationState<ScalarType,MV> getState() const {
00169       CGIterationState<ScalarType,MV> state;
00170       state.R = R_;
00171       state.P = P_;
00172       state.AP = AP_;
00173       state.Z = Z_;
00174       return state;
00175     }
00176     
00178     
00179     
00181 
00182     
00184     int getNumIters() const { return iter_; }
00185     
00187     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00188     
00191     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00192     
00194 
00196     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00197     
00199     
00201 
00202     
00204     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00205     
00207     int getBlockSize() const { return 1; }
00208     
00210     void setBlockSize(int blockSize) { 
00211       TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00212        "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
00213     }
00214     
00216     bool isInitialized() { return initialized_; }
00217     
00219     
00220   private:
00221     
00222     //
00223     // Classes inputed through constructor that define the linear problem to be solved.
00224     //
00225     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00226     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00227     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00228     
00229     //
00230     // Algorithmic parameters
00231     //  
00232     // numRHS_ is the current number of linear systems being solved.
00233     int numRHS_;
00234     
00235     // 
00236     // Current solver state
00237     //
00238     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00239     // is capable of running; _initialize is controlled  by the initialize() member method
00240     // For the implications of the state of initialized_, please see documentation for initialize()
00241     bool initialized_;
00242     
00243     // Current number of iterations performed.
00244     int iter_;
00245 
00246     // Current number of iterations performed.
00247     bool assertPositiveDefiniteness_;
00248 
00249     // 
00250     // State Storage
00251     //
00252     // Residual
00253     Teuchos::RCP<MV> R_;
00254     //
00255     // Preconditioned residual
00256     Teuchos::RCP<MV> Z_;
00257     //
00258     // Direction vector
00259     Teuchos::RCP<MV> P_;
00260     //
00261     // Operator applied to direction vector
00262     Teuchos::RCP<MV> AP_;
00263 
00264   };
00265   
00267   // Constructor.
00268   template<class ScalarType, class MV, class OP>
00269   PseudoBlockCGIter<ScalarType,MV,OP>::PseudoBlockCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00270                      const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00271                      const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00272                      Teuchos::ParameterList &params ):
00273     lp_(problem),
00274     om_(printer),
00275     stest_(tester),
00276     numRHS_(0),
00277     initialized_(false),
00278     iter_(0),
00279     assertPositiveDefiniteness_( params.get("Assert Positive Definiteness", true) )
00280   {
00281   }
00282   
00283 
00285   // Initialize this iteration object
00286   template <class ScalarType, class MV, class OP>
00287   void PseudoBlockCGIter<ScalarType,MV,OP>::initializeCG(CGIterationState<ScalarType,MV> newstate)
00288   {
00289     // Check if there is any multivector to clone from.
00290     Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
00291     Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
00292     TEUCHOS_TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
00293            "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
00294 
00295     // Get the multivector that is not null.
00296     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00297 
00298     // Get the number of right-hand sides we're solving for now.
00299     int numRHS = MVT::GetNumberVecs(*tmp);
00300     numRHS_ = numRHS;
00301   
00302     // Initialize the state storage
00303     // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
00304     if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
00305       R_ = MVT::Clone( *tmp, numRHS_ );
00306       Z_ = MVT::Clone( *tmp, numRHS_ );
00307       P_ = MVT::Clone( *tmp, numRHS_ );
00308       AP_ = MVT::Clone( *tmp, numRHS_ );
00309     }
00310   
00311     // NOTE:  In CGIter R_, the initial residual, is required!!!  
00312     //
00313     std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
00314 
00315     // Create convenience variables for zero and one.
00316     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00317     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00318 
00319     if (!Teuchos::is_null(newstate.R)) {
00320 
00321       TEUCHOS_TEST_FOR_EXCEPTION( MVText::GetGlobalLength(*newstate.R) != MVText::GetGlobalLength(*R_),
00322                           std::invalid_argument, errstr );
00323       TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
00324                           std::invalid_argument, errstr );
00325 
00326       // Copy basis vectors from newstate into V
00327       if (newstate.R != R_) {
00328         // copy over the initial residual (unpreconditioned).
00329   MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
00330       }
00331 
00332       // Compute initial direction vectors
00333       // Initially, they are set to the preconditioned residuals
00334       //
00335       if ( lp_->getLeftPrec() != Teuchos::null ) {
00336         lp_->applyLeftPrec( *R_, *Z_ );
00337         if ( lp_->getRightPrec() != Teuchos::null ) {
00338           Teuchos::RCP<MV> tmp1 = MVT::Clone( *Z_, numRHS_ );
00339           lp_->applyRightPrec( *Z_, *tmp1 );
00340           Z_ = tmp1;
00341         }
00342       }
00343       else if ( lp_->getRightPrec() != Teuchos::null ) {
00344   lp_->applyRightPrec( *R_, *Z_ ); 
00345       } 
00346       else {
00347   Z_ = R_;
00348       }
00349       MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00350     }
00351     else {
00352 
00353       TEUCHOS_TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
00354                          "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
00355     }
00356 
00357     // The solver is initialized
00358     initialized_ = true;
00359   }
00360 
00361 
00363   // Iterate until the status test informs us we should stop.
00364   template <class ScalarType, class MV, class OP>
00365   void PseudoBlockCGIter<ScalarType,MV,OP>::iterate()
00366   {
00367     //
00368     // Allocate/initialize data structures
00369     //
00370     if (initialized_ == false) {
00371       initialize();
00372     }
00373 
00374     // Allocate memory for scalars.
00375     int i=0;
00376     std::vector<int> index(1);
00377     std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
00378     Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ );
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 std::vector.
00385     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00386 
00387     // Compute first <r,z> a.k.a. rHz
00388     MVT::MvDot( *R_, *Z_, rHz );
00389 
00390     if ( assertPositiveDefiniteness_ )
00391         for (i=0; i<numRHS_; ++i)
00392             TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
00393                                 CGIterateFailure,
00394                                 "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
00395 
00397     // Iterate until the status test tells us to stop.
00398     //
00399     while (stest_->checkStatus(this) != Passed) {
00400       
00401       // Increment the iteration
00402       iter_++;
00403 
00404       // Multiply the current direction std::vector by A and store in AP_
00405       lp_->applyOp( *P_, *AP_ );
00406       
00407       // Compute alpha := <R_,Z_> / <P_,AP_>
00408       MVT::MvDot( *P_, *AP_, pAp );
00409 
00410       for (i=0; i<numRHS_; ++i) {
00411         if ( assertPositiveDefiniteness_ )
00412             // Check that pAp[i] is a positive number!
00413             TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(pAp[i]) <= zero,
00414                                 CGIterateFailure,
00415                                 "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00416 
00417         alpha(i,i) = rHz[i] / pAp[i];
00418       }
00419 
00420       //
00421       // Update the solution std::vector x := x + alpha * P_
00422       //
00423       MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
00424       lp_->updateSolution();
00425       //
00426       // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
00427       //
00428       for (i=0; i<numRHS_; ++i) {
00429         rHz_old[i] = rHz[i];
00430       }
00431       //
00432       // Compute the new residual R_ := R_ - alpha * AP_
00433       //
00434       MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
00435       //
00436       // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ], 
00437       // and the new direction std::vector p.
00438       //
00439       if ( lp_->getLeftPrec() != Teuchos::null ) {
00440         lp_->applyLeftPrec( *R_, *Z_ );
00441         if ( lp_->getRightPrec() != Teuchos::null ) {
00442           Teuchos::RCP<MV> tmp = MVT::Clone( *Z_, numRHS_ );
00443           lp_->applyRightPrec( *Z_, *tmp );
00444           Z_ = tmp;
00445         }
00446       }
00447       else if ( lp_->getRightPrec() != Teuchos::null ) {
00448         lp_->applyRightPrec( *R_, *Z_ );
00449       } 
00450       else {
00451         Z_ = R_;
00452       }
00453       //
00454       MVT::MvDot( *R_, *Z_, rHz );
00455       if ( assertPositiveDefiniteness_ )
00456           for (i=0; i<numRHS_; ++i)
00457               TEUCHOS_TEST_FOR_EXCEPTION( SCT::real(rHz[i]) < zero,
00458                                   CGIterateFailure,
00459                                   "Belos::PseudoBlockCGIter::iterate(): negative value for r^H*M*r encountered!" );
00460       //
00461       // Update the search directions.
00462       for (i=0; i<numRHS_; ++i) {
00463         beta(i,i) = rHz[i] / rHz_old[i];
00464   index[0] = i;
00465   Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
00466   Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
00467         MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
00468       }
00469       //      
00470     } // end while (sTest_->checkStatus(this) != Passed)
00471   }
00472 
00473 } // end Belos namespace
00474 
00475 #endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines