Belos Package Browser (Single Doxygen Collection) Development
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     // 
00246     // State Storage
00247     //
00248     // Residual
00249     Teuchos::RCP<MV> R_;
00250     //
00251     // Preconditioned residual
00252     Teuchos::RCP<MV> Z_;
00253     //
00254     // Direction vector
00255     Teuchos::RCP<MV> P_;
00256     //
00257     // Operator applied to direction vector
00258     Teuchos::RCP<MV> AP_;
00259 
00260   };
00261   
00263   // Constructor.
00264   template<class ScalarType, class MV, class OP>
00265   PseudoBlockCGIter<ScalarType,MV,OP>::PseudoBlockCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00266                      const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00267                      const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00268                      Teuchos::ParameterList &params ):
00269     lp_(problem),
00270     om_(printer),
00271     stest_(tester),
00272     numRHS_(0),
00273     initialized_(false),
00274     iter_(0)
00275   {
00276   }
00277   
00278 
00280   // Initialize this iteration object
00281   template <class ScalarType, class MV, class OP>
00282   void PseudoBlockCGIter<ScalarType,MV,OP>::initializeCG(CGIterationState<ScalarType,MV> newstate)
00283   {
00284     // Check if there is any multivector to clone from.
00285     Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
00286     Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
00287     TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
00288            "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
00289 
00290     // Get the multivector that is not null.
00291     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00292 
00293     // Get the number of right-hand sides we're solving for now.
00294     int numRHS = MVT::GetNumberVecs(*tmp);
00295     numRHS_ = numRHS;
00296   
00297     // Initialize the state storage
00298     // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
00299     if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
00300       R_ = MVT::Clone( *tmp, numRHS_ );
00301       Z_ = MVT::Clone( *tmp, numRHS_ );
00302       P_ = MVT::Clone( *tmp, numRHS_ );
00303       AP_ = MVT::Clone( *tmp, numRHS_ );
00304     }
00305   
00306     // NOTE:  In CGIter R_, the initial residual, is required!!!  
00307     //
00308     std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
00309 
00310     // Create convenience variables for zero and one.
00311     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00312     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00313 
00314     if (!Teuchos::is_null(newstate.R)) {
00315 
00316       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
00317                           std::invalid_argument, errstr );
00318       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
00319                           std::invalid_argument, errstr );
00320 
00321       // Copy basis vectors from newstate into V
00322       if (newstate.R != R_) {
00323         // copy over the initial residual (unpreconditioned).
00324   MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
00325       }
00326 
00327       // Compute initial direction vectors
00328       // Initially, they are set to the preconditioned residuals
00329       //
00330       if ( lp_->getLeftPrec() != Teuchos::null ) {
00331   lp_->applyLeftPrec( *R_, *Z_ ); 
00332   MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00333       } else {
00334   Z_ = R_;
00335   MVT::MvAddMv( one, *R_, zero, *R_, *P_ );
00336       }
00337     }
00338     else {
00339 
00340       TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
00341                          "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
00342     }
00343 
00344     // The solver is initialized
00345     initialized_ = true;
00346   }
00347 
00348 
00350   // Iterate until the status test informs us we should stop.
00351   template <class ScalarType, class MV, class OP>
00352   void PseudoBlockCGIter<ScalarType,MV,OP>::iterate()
00353   {
00354     //
00355     // Allocate/initialize data structures
00356     //
00357     if (initialized_ == false) {
00358       initialize();
00359     }
00360 
00361     // Allocate memory for scalars.
00362     int i=0;
00363     std::vector<int> index(1);
00364     std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
00365     Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ );
00366 
00367     // Create convenience variables for zero and one.
00368     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00369     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00370     
00371     // Get the current solution std::vector.
00372     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00373 
00374     // Compute first <r,z> a.k.a. rHz
00375     MVT::MvDot( *R_, *Z_, rHz );
00376     
00378     // Iterate until the status test tells us to stop.
00379     //
00380     while (stest_->checkStatus(this) != Passed) {
00381       
00382       // Increment the iteration
00383       iter_++;
00384 
00385       // Multiply the current direction std::vector by A and store in AP_
00386       lp_->applyOp( *P_, *AP_ );
00387       
00388       // Compute alpha := <R_,Z_> / <P_,AP_>
00389       MVT::MvDot( *P_, *AP_, pAp );
00390       for (i=0; i<numRHS_; ++i) {
00391         alpha(i,i) = rHz[i] / pAp[i];
00392         // Check that alpha is a positive number!
00393         TEST_FOR_EXCEPTION( SCT::real(alpha(i,i)) <= zero, CGIterateFailure,
00394         "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00395       }
00396       
00397       //
00398       // Update the solution std::vector x := x + alpha * P_
00399       //
00400       MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
00401       lp_->updateSolution();
00402       //
00403       // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
00404       //
00405       for (i=0; i<numRHS_; ++i) {
00406         rHz_old[i] = rHz[i];
00407       }
00408       //
00409       // Compute the new residual R_ := R_ - alpha * AP_
00410       //
00411       MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
00412       //
00413       // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ], 
00414       // and the new direction std::vector p.
00415       //
00416       if ( lp_->getLeftPrec() != Teuchos::null ) {
00417   lp_->applyLeftPrec( *R_, *Z_ );
00418       } else {
00419   Z_ = R_;
00420       }
00421       //
00422       MVT::MvDot( *R_, *Z_, rHz );
00423       //
00424       // Update the search directions.
00425       for (i=0; i<numRHS_; ++i) {
00426         beta(i,i) = rHz[i] / rHz_old[i];
00427   index[0] = i;
00428   Teuchos::RCP<const MV> Z_i = MVT::CloneView( *Z_, index );
00429   Teuchos::RCP<MV> P_i = MVT::CloneViewNonConst( *P_, index );
00430         MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
00431       }
00432       //      
00433     } // end while (sTest_->checkStatus(this) != Passed)
00434   }
00435 
00436 } // end Belos namespace
00437 
00438 #endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines