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 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_PSEUDO_BLOCK_CG_ITER_HPP
00030 #define BELOS_PSEUDO_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_SerialDenseMatrix.hpp"
00049 #include "Teuchos_SerialDenseVector.hpp"
00050 #include "Teuchos_ScalarTraits.hpp"
00051 #include "Teuchos_ParameterList.hpp"
00052 #include "Teuchos_TimeMonitor.hpp"
00053 
00065 namespace Belos {
00066   
00067   template<class ScalarType, class MV, class OP>
00068   class PseudoBlockCGIter : 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     PseudoBlockCGIter( 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         Teuchos::ParameterList &params );
00092     
00094     virtual ~PseudoBlockCGIter() {};
00096     
00097     
00099 
00100     
00114     void iterate();
00115     
00136     void initializeCG(CGIterationState<ScalarType,MV> newstate);
00137     
00141     void initialize()
00142     {
00143       CGIterationState<ScalarType,MV> empty;
00144       initializeCG(empty);
00145     }
00146     
00154     CGIterationState<ScalarType,MV> getState() const {
00155       CGIterationState<ScalarType,MV> state;
00156       state.R = R_;
00157       state.P = P_;
00158       state.AP = AP_;
00159       state.Z = Z_;
00160       return state;
00161     }
00162     
00164     
00165     
00167 
00168     
00170     int getNumIters() const { return iter_; }
00171     
00173     void resetNumIters( int iter = 0 ) { iter_ = iter; }
00174     
00177     Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }
00178     
00180 
00182     Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }
00183     
00185     
00187 
00188     
00190     const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
00191     
00193     int getBlockSize() const { return 1; }
00194     
00196     void setBlockSize(int blockSize) { 
00197       TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
00198        "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
00199     }
00200     
00202     bool isInitialized() { return initialized_; }
00203     
00205     
00206   private:
00207     
00208     //
00209     // Classes inputed through constructor that define the linear problem to be solved.
00210     //
00211     const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
00212     const Teuchos::RCP<OutputManager<ScalarType> >          om_;
00213     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
00214     
00215     //
00216     // Algorithmic parameters
00217     //  
00218     // numRHS_ is the current number of linear systems being solved.
00219     int numRHS_;
00220     
00221     // 
00222     // Current solver state
00223     //
00224     // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
00225     // is capable of running; _initialize is controlled  by the initialize() member method
00226     // For the implications of the state of initialized_, please see documentation for initialize()
00227     bool initialized_;
00228     
00229     // Current number of iterations performed.
00230     int iter_;
00231 
00232     // 
00233     // State Storage
00234     //
00235     // Residual
00236     Teuchos::RCP<MV> R_;
00237     //
00238     // Preconditioned residual
00239     Teuchos::RCP<MV> Z_;
00240     //
00241     // Direction vector
00242     Teuchos::RCP<MV> P_;
00243     //
00244     // Operator applied to direction vector
00245     Teuchos::RCP<MV> AP_;
00246 
00247   };
00248   
00250   // Constructor.
00251   template<class ScalarType, class MV, class OP>
00252   PseudoBlockCGIter<ScalarType,MV,OP>::PseudoBlockCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
00253                      const Teuchos::RCP<OutputManager<ScalarType> > &printer,
00254                      const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
00255                      Teuchos::ParameterList &params ):
00256     lp_(problem),
00257     om_(printer),
00258     stest_(tester),
00259     numRHS_(0),
00260     initialized_(false),
00261     iter_(0)
00262   {
00263   }
00264   
00265 
00267   // Initialize this iteration object
00268   template <class ScalarType, class MV, class OP>
00269   void PseudoBlockCGIter<ScalarType,MV,OP>::initializeCG(CGIterationState<ScalarType,MV> newstate)
00270   {
00271     // Check if there is any multivector to clone from.
00272     Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
00273     Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
00274     TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
00275            "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");
00276 
00277     // Get the multivector that is not null.
00278     Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
00279 
00280     // Get the number of right-hand sides we're solving for now.
00281     int numRHS = MVT::GetNumberVecs(*tmp);
00282     numRHS_ = numRHS;
00283   
00284     // Initialize the state storage
00285     // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
00286     if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
00287       R_ = MVT::Clone( *tmp, numRHS_ );
00288       Z_ = MVT::Clone( *tmp, numRHS_ );
00289       P_ = MVT::Clone( *tmp, numRHS_ );
00290       AP_ = MVT::Clone( *tmp, numRHS_ );
00291     }
00292   
00293     // NOTE:  In CGIter R_, the initial residual, is required!!!  
00294     //
00295     std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");
00296 
00297     // Create convenience variables for zero and one.
00298     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00299     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00300 
00301     if (!Teuchos::is_null(newstate.R)) {
00302 
00303       TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
00304                           std::invalid_argument, errstr );
00305       TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
00306                           std::invalid_argument, errstr );
00307 
00308       // Copy basis vectors from newstate into V
00309       if (newstate.R != R_) {
00310         // copy over the initial residual (unpreconditioned).
00311   MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
00312       }
00313 
00314       // Compute initial direction vectors
00315       // Initially, they are set to the preconditioned residuals
00316       //
00317       if ( lp_->getLeftPrec() != Teuchos::null ) {
00318   lp_->applyLeftPrec( *R_, *Z_ ); 
00319   MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
00320       } else {
00321   Z_ = R_;
00322   MVT::MvAddMv( one, *R_, zero, *R_, *P_ );
00323       }
00324     }
00325     else {
00326 
00327       TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
00328                          "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
00329     }
00330 
00331     // The solver is initialized
00332     initialized_ = true;
00333   }
00334 
00335 
00337   // Iterate until the status test informs us we should stop.
00338   template <class ScalarType, class MV, class OP>
00339   void PseudoBlockCGIter<ScalarType,MV,OP>::iterate()
00340   {
00341     //
00342     // Allocate/initialize data structures
00343     //
00344     if (initialized_ == false) {
00345       initialize();
00346     }
00347 
00348     // Allocate memory for scalars.
00349     int i=0;
00350     std::vector<int> index(1);
00351     std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
00352     Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ );
00353 
00354     // Create convenience variables for zero and one.
00355     const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
00356     const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
00357     
00358     // Get the current solution std::vector.
00359     Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();
00360 
00361     // Compute first <r,z> a.k.a. rHz
00362     MVT::MvDot( *R_, *Z_, rHz );
00363     
00365     // Iterate until the status test tells us to stop.
00366     //
00367     while (stest_->checkStatus(this) != Passed) {
00368       
00369       // Increment the iteration
00370       iter_++;
00371 
00372       // Multiply the current direction std::vector by A and store in AP_
00373       lp_->applyOp( *P_, *AP_ );
00374       
00375       // Compute alpha := <R_,Z_> / <P_,AP_>
00376       MVT::MvDot( *P_, *AP_, pAp );
00377       for (i=0; i<numRHS_; ++i) {
00378         alpha(i,i) = rHz[i] / pAp[i];
00379         // Check that alpha is a positive number!
00380         TEST_FOR_EXCEPTION( SCT::real(alpha(i,i)) <= zero, CGIterateFailure,
00381         "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
00382       }
00383       
00384       //
00385       // Update the solution std::vector x := x + alpha * P_
00386       //
00387       MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
00388       lp_->updateSolution();
00389       //
00390       // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
00391       //
00392       for (i=0; i<numRHS_; ++i) {
00393         rHz_old[i] = rHz[i];
00394       }
00395       //
00396       // Compute the new residual R_ := R_ - alpha * AP_
00397       //
00398       MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
00399       //
00400       // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ], 
00401       // and the new direction std::vector p.
00402       //
00403       if ( lp_->getLeftPrec() != Teuchos::null ) {
00404   lp_->applyLeftPrec( *R_, *Z_ );
00405       } else {
00406   Z_ = R_;
00407       }
00408       //
00409       MVT::MvDot( *R_, *Z_, rHz );
00410       //
00411       // Update the search directions.
00412       for (i=0; i<numRHS_; ++i) {
00413         beta(i,i) = rHz[i] / rHz_old[i];
00414   index[0] = i;
00415   Teuchos::RCP<MV> Z_i = MVT::CloneView( *Z_, index );
00416   Teuchos::RCP<MV> P_i = MVT::CloneView( *P_, index );
00417         MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
00418       }
00419       //      
00420     } // end while (sTest_->checkStatus(this) != Passed)
00421   }
00422 
00423 } // end Belos namespace
00424 
00425 #endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */

Generated on Tue Jul 13 09:27:02 2010 for Belos by  doxygen 1.4.7