BelosPCPGIter.hpp

Go to the documentation of this file.

//@HEADER
// ************************************************************************
//
//                 Belos: Block Linear Solvers Package
//                  Copyright 2004 Sandia Corporation
//
// Under the terms of Contract DE-AC04-94AL85000 with Sandia Corporation,
// the U.S. Government retains certain rights in this software.
//
// Redistribution and use in source and binary forms, with or without
// modification, are permitted provided that the following conditions are
// met:
//
// 1. Redistributions of source code must retain the above copyright
// notice, this list of conditions and the following disclaimer.
//
// 2. Redistributions in binary form must reproduce the above copyright
// notice, this list of conditions and the following disclaimer in the
// documentation and/or other materials provided with the distribution.
//
// 3. Neither the name of the Corporation nor the names of the
// contributors may be used to endorse or promote products derived from
// this software without specific prior written permission.
//
// THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
// EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
// IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
// PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
// CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
// EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
// PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
// PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
// LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
// NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
// SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
//
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ************************************************************************
//@HEADER

#ifndef BELOS_PCPG_ITER_HPP
#define BELOS_PCPG_ITER_HPP

#include "BelosConfigDefs.hpp"
#include "BelosTypes.hpp"

#include "BelosLinearProblem.hpp"
#include "BelosMatOrthoManager.hpp"
#include "BelosOutputManager.hpp"
#include "BelosStatusTest.hpp"
#include "BelosOperatorTraits.hpp"
#include "BelosMultiVecTraits.hpp"

#include "Teuchos_BLAS.hpp"
#include "Teuchos_LAPACK.hpp"
#include "Teuchos_SerialDenseMatrix.hpp"
#include "Teuchos_SerialDenseVector.hpp"
#include "Teuchos_ScalarTraits.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_TimeMonitor.hpp"

namespace Belos {
  

  
  template <class ScalarType, class MV>
  struct PCPGIterState {
    int curDim;
    int prevUdim;

    Teuchos::RCP<MV> R;

    Teuchos::RCP<MV> Z;

    Teuchos::RCP<MV> P;

    Teuchos::RCP<MV> AP;

    Teuchos::RCP<MV> U;

    Teuchos::RCP<MV> C;

    Teuchos::RCP<const Teuchos::SerialDenseMatrix<int,ScalarType> > D;

    PCPGIterState() : curDim(0), 
                      prevUdim(0), 
                      R(Teuchos::null), Z(Teuchos::null), 
                      P(Teuchos::null), AP(Teuchos::null),
                      U(Teuchos::null), C(Teuchos::null),
          D(Teuchos::null)
    {}
  };
  
  

  
  class PCPGIterInitFailure : public BelosError {public:
    PCPGIterInitFailure(const std::string& what_arg) : BelosError(what_arg)
    {}};

  class PCPGIterateFailure : public BelosError {public:
    PCPGIterateFailure(const std::string& what_arg) : BelosError(what_arg)
    {}};


  
  class PCPGIterOrthoFailure : public BelosError {public:
    PCPGIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg)
    {}};
  
  class PCPGIterLAPACKFailure : public BelosError {public:
    PCPGIterLAPACKFailure(const std::string& what_arg) : BelosError(what_arg)
    {}};
  
  
  
  template<class ScalarType, class MV, class OP>
  class PCPGIter : virtual public Iteration<ScalarType,MV,OP> {
    
  public:
    
    //
    // Convenience typedefs
    //
    typedef MultiVecTraits<ScalarType,MV> MVT;
    typedef OperatorTraits<ScalarType,MV,OP> OPT;
    typedef Teuchos::ScalarTraits<ScalarType> SCT;
    typedef typename SCT::magnitudeType MagnitudeType;
    

    
    PCPGIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
    const Teuchos::RCP<OutputManager<ScalarType> > &printer,
    const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
    const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
    Teuchos::ParameterList &params );
    
    virtual ~PCPGIter() {};
    
    

    
    void iterate();
    
    void initialize(PCPGIterState<ScalarType,MV> newstate);
    
    void initialize()
    {
      PCPGIterState<ScalarType,MV> empty;
      initialize(empty);
    }
    
    PCPGIterState<ScalarType,MV> getState() const {
      PCPGIterState<ScalarType,MV> state;
      state.Z = Z_;         // CG state
      state.P = P_;
      state.AP = AP_;
      state.R = R_;
      state.U = U_;         // seed state 
      state.C = C_; 
      state.D = D_;
      state.curDim = curDim_;
      state.prevUdim = prevUdim_;
      return state;
    }
    
    
    

    
    int getNumIters() const { return iter_; }
    
    void resetNumIters( int iter = 0 ) { iter_ = iter; }
    
    Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const { return R_; }


    Teuchos::RCP<MV> getCurrentUpdate() const { return Teuchos::null; }

    int getCurSubspaceDim() const { 
      if (!initialized_) return 0;
      return curDim_;
    };

    int getPrevSubspaceDim() const { 
      if (!initialized_) return 0;
      return prevUdim_;
    };
    
    
    

    
    const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
    
    int getBlockSize() const { return 1; }
    
    int getNumRecycledBlocks() const { return savedBlocks_; }

    
    void setBlockSize(int blockSize) {
      TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
       "Belos::PCPGIter::setBlockSize(): Cannot use a block size that is not one.");
    }

    void setSize( int savedBlocks );

    bool isInitialized() { return initialized_; }

    void resetState(); 
    
    
  private:
    
    //
    // Internal methods
    //
    void setStateSize();
    
    //
    // Classes inputed through constructor that define the linear problem to be solved.
    //
    const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> >    lp_;
    const Teuchos::RCP<OutputManager<ScalarType> >          om_;
    const Teuchos::RCP<StatusTest<ScalarType,MV,OP> >       stest_;
    const Teuchos::RCP<OrthoManager<ScalarType,MV> >        ortho_;

    //
    // Algorithmic parameters
    // savedBlocks_ is the number of blocks allocated for the reused subspace
    int savedBlocks_; 
    //
    // 
    // Current solver state
    //
    // initialized_ specifies that the basis vectors have been initialized and the iterate() routine
    // is capable of running; _initialize is controlled  by the initialize() member method
    // For the implications of the state of initialized_, please see documentation for initialize()
    bool initialized_;
    
    // stateStorageInitialized_ indicates that the state storage has be initialized to the current
    // savedBlocks_.  State storage initialization may be postponed if the linear problem was
    // generated without either the right-hand side or solution vectors.
    bool stateStorageInitialized_;

    // keepDiagonal_ specifies that the iteration must keep the diagonal matrix of pivots 
    bool keepDiagonal_;

    // initDiagonal_ specifies that the iteration will reinitialize the diagonal matrix by zeroing
    // out all entries before an iteration is started.
    bool initDiagonal_;
    
    // Current subspace dimension
    int curDim_;

    // Dimension of seed space used to solve current linear system
    int prevUdim_;
    
    // Number of iterations performed
    int iter_;
    // 
    // State Storage    ... of course this part is different for CG
    //
    // Residual
    Teuchos::RCP<MV> R_;
    //
    // Preconditioned residual
    Teuchos::RCP<MV> Z_;
    //
    // Direction std::vector
    Teuchos::RCP<MV> P_;
    //
    // Operator applied to direction std::vector
    Teuchos::RCP<MV> AP_;
    //
    // Recycled subspace vectors.
    Teuchos::RCP<MV> U_;
    // 
    // C = A * U,  linear system is Ax=b 
    Teuchos::RCP<MV> C_;
    //
    // Projected matrices
    // D_ : Diagonal matrix of pivots D = P'AP 
    Teuchos::RCP<Teuchos::SerialDenseMatrix<int,ScalarType> > D_;
  };
  
  // Constructor.
  template<class ScalarType, class MV, class OP>
  PCPGIter<ScalarType,MV,OP>::PCPGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
             const Teuchos::RCP<OutputManager<ScalarType> > &printer,
             const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
             const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP> > &ortho,
             Teuchos::ParameterList &params ):
    lp_(problem),
    om_(printer),
    stest_(tester),
    ortho_(ortho),
    savedBlocks_(0),
    initialized_(false),
    stateStorageInitialized_(false),
    keepDiagonal_(false), 
    initDiagonal_(false),
    curDim_(0),
    prevUdim_(0),
    iter_(0)
  {
    // Get the maximum number of blocks allowed for this Krylov subspace

    TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Saved Blocks"), std::invalid_argument,
                       "Belos::PCPGIter::constructor: mandatory parameter \"Saved Blocks\" is not specified.");
    int rb = Teuchos::getParameter<int>(params, "Saved Blocks");

    // Find out whether we are saving the Diagonal matrix.
    keepDiagonal_ = params.get("Keep Diagonal", false);

    // Find out whether we are initializing the Diagonal matrix.
    initDiagonal_ = params.get("Initialize Diagonal", false);

    // Set the number of blocks and allocate data
    setSize( rb );
  }
  
  // Set the block size and adjust as necessary
  template <class ScalarType, class MV, class OP>
  void PCPGIter<ScalarType,MV,OP>::setSize( int savedBlocks )
  {
    // allocate space only; perform no computation
    // Any change in size invalidates the state of the solver as implemented here.

    TEUCHOS_TEST_FOR_EXCEPTION(savedBlocks <= 0, std::invalid_argument, "Belos::PCPGIter::setSize() was passed a non-positive argument for \"Num Saved Blocks\".");

    if ( savedBlocks_ != savedBlocks) {
      stateStorageInitialized_ = false;
      savedBlocks_ = savedBlocks;
      initialized_ = false;
      curDim_ = 0;
      prevUdim_ = 0;
      setStateSize(); // Use the current savedBlocks_ to initialize the state storage.    
    }
  }

  // Enable the reuse of a single solver object for completely different linear systems
  template <class ScalarType, class MV, class OP>
  void PCPGIter<ScalarType,MV,OP>::resetState()
  {
      stateStorageInitialized_ = false;
      initialized_ = false;
      curDim_ = 0;
      prevUdim_ = 0;
      setStateSize();
  }

  // Setup the state storage.  Called by either initialize or, if savedBlocks_ changes, setSize.
  template <class ScalarType, class MV, class OP>
  void PCPGIter<ScalarType,MV,OP>::setStateSize ()
  {
    if (!stateStorageInitialized_) {

      // Check if there is any multivector to clone from.
      Teuchos::RCP<const MV> lhsMV = lp_->getLHS();
      Teuchos::RCP<const MV> rhsMV = lp_->getRHS();
      if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
  return;  // postpone exception 
      }
      else {
  
  // blockSize*recycledBlocks dependent
        int newsd = savedBlocks_ ; //int newsd = blockSize_* savedBlocks_ ;
  //
  // Initialize the CG state storage
        // If the subspace is not initialized, generate it using the LHS or RHS from lp_.
  // Generate CG state only if it does not exist, otherwise resize it.
        if (Z_ == Teuchos::null) {
          Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
          Z_ = MVT::Clone( *tmp, 1 );
        }
        if (P_ == Teuchos::null) {
          Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
          P_ = MVT::Clone( *tmp, 1 );
        }
        if (AP_ == Teuchos::null) {
          Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
          AP_ = MVT::Clone( *tmp, 1 );
        }

  if (C_ == Teuchos::null) {        

    // Get the multivector that is not null. 
    Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
    TEUCHOS_TEST_FOR_EXCEPTION(tmp == Teuchos::null,std::invalid_argument,
           "Belos::PCPGIter::setStateSize(): linear problem does not specify multivectors to clone from.");
    TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
           "Belos::PCPGIter::setStateSize(): prevUdim not zero and C is null.");
    C_ = MVT::Clone( *tmp, savedBlocks_ );
  }
  else {
    // Generate C_ by cloning itself ONLY if more space is needed.
    if (MVT::GetNumberVecs(*C_) < savedBlocks_ ) {
      Teuchos::RCP<const MV> tmp = C_;
      C_ = MVT::Clone( *tmp, savedBlocks_ );
    }
  }
  if (U_ == Teuchos::null) {        
    Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );
    TEUCHOS_TEST_FOR_EXCEPTION( 0 != prevUdim_,std::invalid_argument,
           "Belos::PCPGIter::setStateSize(): prevUdim not zero and U is null.");
    U_ = MVT::Clone( *tmp, savedBlocks_ );
  }
  else {
    // Generate U_ by cloning itself ONLY if more space is needed.
    if (MVT::GetNumberVecs(*U_) < savedBlocks_ ) {
      Teuchos::RCP<const MV> tmp = U_;
      U_ = MVT::Clone( *tmp, savedBlocks_ );
    }
  }
        if (keepDiagonal_) {
          if (D_ == Teuchos::null) {
            D_ = Teuchos::rcp( new Teuchos::SerialDenseMatrix<int,ScalarType>() );
          }
          if (initDiagonal_) {
            D_->shape( newsd, newsd );
          }
          else {
            if (D_->numRows() < newsd || D_->numCols() < newsd) {
              D_->shapeUninitialized( newsd, newsd );
            }
          }
        }
  // State storage has now been initialized.
  stateStorageInitialized_ = true;
      } // if there is a vector to clone from
    } // if !stateStorageInitialized_ 
  } // end of setStateSize

  // Initialize the iteration object
  template <class ScalarType, class MV, class OP>
  void PCPGIter<ScalarType,MV,OP>::initialize(PCPGIterState<ScalarType,MV> newstate)
  {

    TEUCHOS_TEST_FOR_EXCEPTION(!stateStorageInitialized_,std::invalid_argument,
           "Belos::PCPGIter::initialize(): Cannot initialize state storage!");
    
    // Requirements: R_ and consistent multivectors widths and lengths
    //
    std::string errstr("Belos::PCPGIter::initialize(): Specified multivectors must have a consistent length and width.");

    if (newstate.R != Teuchos::null){ 

      R_ = newstate.R; // SolverManager::R_ == newstate.R == Iterator::R_
      if (newstate.U == Teuchos::null){ 
        prevUdim_ = 0;
        newstate.U = U_;
        newstate.C = C_;
      }
      else {
        prevUdim_ =  newstate.curDim;
        if (newstate.C == Teuchos::null){  // Stub for new feature
          std::vector<int> index(prevUdim_);
          for (int i=0; i< prevUdim_; ++i)  
            index[i] = i; 
          Teuchos::RCP<const MV> Ukeff = MVT::CloneView( *newstate.U, index );
          newstate.C = MVT::Clone( *newstate.U, prevUdim_ ); 
          Teuchos::RCP<MV> Ckeff = MVT::CloneViewNonConst( *newstate.C, index );    
          lp_->apply( *Ukeff, *Ckeff );
        }
        curDim_ = prevUdim_ ;
      }

      // Initialize the state storage if not already allocated in the constructor
      if (!stateStorageInitialized_) 
        setStateSize();

      //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.V) != MVT::GetVecLength(*V_), std::invalid_argument, errstr );
      //TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.V) < 1, std::invalid_argument, errstr );

      newstate.prevUdim =  prevUdim_; // big change in functionality from GCRODR 
      newstate.curDim =  curDim_; 

      //TEUCHOS_TEST_FOR_EXCEPTION(newstate.z->numRows() < curDim_ || newstate.z->numCols() < 1, std::invalid_argument, errstr);

      std::vector<int> zero_index(1);
      zero_index[0] = 0;
      if ( lp_->getLeftPrec() != Teuchos::null ) { // Compute the initial search direction 
        lp_->applyLeftPrec( *R_, *Z_ );
        MVT::SetBlock( *Z_,  zero_index , *P_ );  // P(:,zero_index) := Z
      } else {                                      
        Z_ = R_;
        MVT::SetBlock( *R_, zero_index, *P_ );
      }

      std::vector<int> nextind(1);
      nextind[0] = curDim_;

      MVT::SetBlock( *P_,  nextind, *newstate.U ); // U(:,curDim_ ) := P_

      ++curDim_;
      newstate.curDim = curDim_; 

      TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.U) != savedBlocks_ ,
                          std::invalid_argument, errstr );
      if (newstate.U != U_) { // Why this is needed?
  U_ = newstate.U;
      }

      TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.C) != savedBlocks_ ,
                          std::invalid_argument, errstr );
      if (newstate.C != C_) {
  C_ = newstate.C;
      }
    }
    else {

      TEUCHOS_TEST_FOR_EXCEPTION(newstate.R == Teuchos::null,std::invalid_argument,
                         "Belos::PCPGIter::initialize(): PCPGStateIterState does not have initial kernel R_0.");
    }

    // the solver is initialized
    initialized_ = true;
  }


  // Iterate until the status test informs us we should stop.
  template <class ScalarType, class MV, class OP>
  void PCPGIter<ScalarType,MV,OP>::iterate()
  {
    //
    // Allocate/initialize data structures
    //
    if (initialized_ == false) {
      initialize();
    }
    bool debug = false;

    // Allocate memory for scalars.
    Teuchos::SerialDenseMatrix<int,ScalarType> alpha( 1, 1 );
    Teuchos::SerialDenseMatrix<int,ScalarType> pAp( 1, 1 );
    Teuchos::SerialDenseMatrix<int,ScalarType> beta( 1, 1 );
    Teuchos::SerialDenseMatrix<int,ScalarType> rHz( 1, 1 ), rHz_old( 1, 1 );

    if( iter_ != 0 )
      std::cout << " Iterate Warning: begin from nonzero iter_ ?" << std::endl;  //DMD

    // GenOrtho Project Stubs
    std::vector<int> prevInd;
    Teuchos::RCP<const MV> Uprev;
    Teuchos::RCP<const MV> Cprev;
    Teuchos::SerialDenseMatrix<int,ScalarType> CZ;

    if( prevUdim_ ){
      prevInd.resize( prevUdim_ );
      for( int i=0; i<prevUdim_ ; i++) prevInd[i] = i;
      CZ.reshape( prevUdim_ , 1 );
      Uprev = MVT::CloneView(*U_, prevInd);
      Cprev = MVT::CloneView(*C_, prevInd);
    }

    // Get the current solution std::vector.
    Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();

    // Check that the current solution std::vector only has one column.
    TEUCHOS_TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*cur_soln_vec) != 1, PCPGIterInitFailure,
                        "Belos::CGIter::iterate(): current linear system has more than one std::vector!" );

    //Check that the input is correctly set up 
    TEUCHOS_TEST_FOR_EXCEPTION( curDim_  != prevUdim_ + 1, PCPGIterInitFailure,
                        "Belos::CGIter::iterate(): mistake in initialization !" );


    const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
    const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();


    std::vector<int> curind(1);
    std::vector<ScalarType> rnorm(MVT::GetNumberVecs(*cur_soln_vec));
    if (prevUdim_ > 0){                 // A-orthonalize P=Z to Uprev
      Teuchos::RCP<MV> P; 
      curind[0] = curDim_ - 1;          // column = dimension - 1 
      P = MVT::CloneViewNonConst(*U_,curind); 
      MVT::MvTransMv( one, *Cprev, *P, CZ );
      MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P );       // P -= U*(C'Z)

      if( debug ){
        MVT::MvTransMv( one, *Cprev, *P, CZ );
        std::cout << " Input CZ post ortho " << std::endl;
        CZ.print( std::cout );
      }
      if( curDim_ == savedBlocks_ ){
        std::vector<int> zero_index(1);
        zero_index[0] = 0;
        MVT::SetBlock( *P, zero_index, *P_ );
      }
      P = Teuchos::null;
    }

    // Compute first <r,z> a.k.a. rHz
    MVT::MvTransMv( one, *R_, *Z_, rHz );

    // iterate until the status test is satisfied
    //
    while (stest_->checkStatus(this) != Passed ) {
      Teuchos::RCP<const MV> P; 
      Teuchos::RCP<MV> AP;
      iter_++;                          // The next iteration begins.
      //std::vector<int> curind(1);
      curind[0] = curDim_ - 1;          // column = dimension - 1 
      if( debug ){
        MVT::MvNorm(*R_, rnorm);
        std::cout << iter_ << "  " << curDim_ <<  "   " << rnorm[0] << std::endl;
      }
      if( prevUdim_ + iter_ < savedBlocks_ ){
        P = MVT::CloneView(*U_,curind); 
        AP = MVT::CloneViewNonConst(*C_,curind); 
        lp_->applyOp( *P, *AP );
        MVT::MvTransMv( one, *P, *AP, pAp );
      }else{
        if( prevUdim_ + iter_ == savedBlocks_ ){
          AP = MVT::CloneViewNonConst(*C_,curind); 
          lp_->applyOp( *P_, *AP );
          MVT::MvTransMv( one, *P_, *AP, pAp );
        }else{
          lp_->applyOp( *P_, *AP_ );
          MVT::MvTransMv( one, *P_, *AP_, pAp );
        }
      }

      if( keepDiagonal_  && prevUdim_ + iter_ <= savedBlocks_ )
        (*D_)(iter_ -1 ,iter_ -1 ) = pAp(0,0);

      // positive pAp required 
      TEUCHOS_TEST_FOR_EXCEPTION( pAp(0,0) <= zero, PCPGIterateFailure,
                          "Belos::CGIter::iterate(): non-positive value for p^H*A*p encountered!" );

      // alpha := <R_,Z_> / <P,AP>
      alpha(0,0) = rHz(0,0) / pAp(0,0);

      // positive alpha required 
      TEUCHOS_TEST_FOR_EXCEPTION( alpha(0,0) <= zero, PCPGIterateFailure,
                          "Belos::CGIter::iterate(): non-positive value for alpha encountered!" );

      // solution update  x += alpha * P
      if( curDim_ < savedBlocks_ ){
         MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P, *cur_soln_vec );
      }else{
         MVT::MvAddMv( one, *cur_soln_vec, alpha(0,0), *P_, *cur_soln_vec );
      }
      //lp_->updateSolution(); ... does nothing.
      //
      // The denominator of beta is saved before residual is updated [ old <R_, Z_> ].
      //
      rHz_old(0,0) = rHz(0,0);
      //
      // residual update R_ := R_ - alpha * AP
      //
      if( prevUdim_ + iter_ <= savedBlocks_ ){
         MVT::MvAddMv( one, *R_, -alpha(0,0), *AP, *R_ );
         AP = Teuchos::null;
      }else{
         MVT::MvAddMv( one, *R_, -alpha(0,0), *AP_, *R_ );
      }
      //
      // update beta := [ new <R_, Z_> ] / [ old <R_, Z_> ] and the search direction p.
      //
      if ( lp_->getLeftPrec() != Teuchos::null ) {
        lp_->applyLeftPrec( *R_, *Z_ );
      } else {
        Z_ = R_;
      }
      //
      MVT::MvTransMv( one, *R_, *Z_, rHz );
      //
      beta(0,0) = rHz(0,0) / rHz_old(0,0);
      //
      if( curDim_ < savedBlocks_ ){
         curDim_++;                                                         // update basis dim
         curind[0] = curDim_ - 1;
         Teuchos::RCP<MV> Pnext = MVT::CloneViewNonConst(*U_,curind);
         MVT::MvAddMv( one, *Z_, beta(0,0), *P, *Pnext );
         if( prevUdim_ ){ // Deflate seed space 
             MVT::MvTransMv( one, *Cprev, *Z_, CZ );
             MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *Pnext ); // Pnext -= U*(C'Z)
             if( debug ){
               std::cout << " Check CZ " << std::endl;
               MVT::MvTransMv( one, *Cprev, *Pnext, CZ );
               CZ.print( std::cout );
             }
         }
         P = Teuchos::null;
         if( curDim_ == savedBlocks_ ){
           std::vector<int> zero_index(1);
           zero_index[0] = 0;
           MVT::SetBlock( *Pnext, zero_index, *P_ );
         }
         Pnext = Teuchos::null;
      }else{
         MVT::MvAddMv( one, *Z_, beta(0,0), *P_, *P_ );
         if( prevUdim_ ){ // Deflate seed space
             MVT::MvTransMv( one, *Cprev, *Z_, CZ );
             MVT::MvTimesMatAddMv( -one, *Uprev, CZ, one, *P_ );       // P_ -= U*(C'Z)

             if( debug ){
               std::cout << " Check CZ " << std::endl;
               MVT::MvTransMv( one, *Cprev, *P_, CZ );
               CZ.print( std::cout );
             }
         }
      }
      // CGB: 5/26/2010
      // this RCP<const MV> P was previously a variable outside the loop. however, it didn't appear to be see any use between
      // loop iterations. therefore, I moved it inside to avoid scoping errors with previously used variables named P.
      // to ensure that this wasn't a bug, I verify below that we have set P == null, i.e., that we are not going to use it again
      // same for AP
      TEUCHOS_TEST_FOR_EXCEPTION( AP != Teuchos::null || P != Teuchos::null, std::logic_error, "Loop recurrence violated. Please contact Belos team.");
    } // end coupled two-term recursion
    if( prevUdim_ + iter_ < savedBlocks_ ) --curDim_; // discard negligible search direction
  }

} // end Belos namespace

#endif /* BELOS_PCPG_ITER_HPP */