BelosPseudoBlockCGIter.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
//
//                 Belos: Block Linear Solvers Package
//                 Copyright (2004) Sandia Corporation
//
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
//
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov)
//
// ***********************************************************************
// @HEADER

#ifndef BELOS_PSEUDO_BLOCK_CG_ITER_HPP
#define BELOS_PSEUDO_BLOCK_CG_ITER_HPP

#include "BelosConfigDefs.hpp"
#include "BelosTypes.hpp"
#include "BelosCGIteration.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_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, class OP>
  class PseudoBlockCGIter : virtual public CGIteration<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;
    

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

    
    void iterate();
    
    void initializeCG(CGIterationState<ScalarType,MV> newstate);
    
    void initialize()
    {
      CGIterationState<ScalarType,MV> empty;
      initializeCG(empty);
    }
    
    CGIterationState<ScalarType,MV> getState() const {
      CGIterationState<ScalarType,MV> state;
      state.R = R_;
      state.P = P_;
      state.AP = AP_;
      state.Z = Z_;
      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; }
    
    

    
    const LinearProblem<ScalarType,MV,OP>& getProblem() const { return *lp_; }
    
    int getBlockSize() const { return 1; }
    
    void setBlockSize(int blockSize) { 
      TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,
       "Belos::PseudoBlockCGIter::setBlockSize(): Cannot use a block size that is not one.");
    }
    
    bool isInitialized() { return initialized_; }
    
    
  private:
    
    //
    // 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_;
    
    //
    // Algorithmic parameters
    //  
    // numRHS_ is the current number of linear systems being solved.
    int numRHS_;
    
    // 
    // 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_;
    
    // Current number of iterations performed.
    int iter_;

    // 
    // State Storage
    //
    // Residual
    Teuchos::RCP<MV> R_;
    //
    // Preconditioned residual
    Teuchos::RCP<MV> Z_;
    //
    // Direction vector
    Teuchos::RCP<MV> P_;
    //
    // Operator applied to direction vector
    Teuchos::RCP<MV> AP_;

  };
  
  // Constructor.
  template<class ScalarType, class MV, class OP>
  PseudoBlockCGIter<ScalarType,MV,OP>::PseudoBlockCGIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP> > &problem, 
                     const Teuchos::RCP<OutputManager<ScalarType> > &printer,
                     const Teuchos::RCP<StatusTest<ScalarType,MV,OP> > &tester,
                     Teuchos::ParameterList &params ):
    lp_(problem),
    om_(printer),
    stest_(tester),
    numRHS_(0),
    initialized_(false),
    iter_(0)
  {
  }
  

  // Initialize this iteration object
  template <class ScalarType, class MV, class OP>
  void PseudoBlockCGIter<ScalarType,MV,OP>::initializeCG(CGIterationState<ScalarType,MV> newstate)
  {
    // Check if there is any multivector to clone from.
    Teuchos::RCP<const MV> lhsMV = lp_->getCurrLHSVec();
    Teuchos::RCP<const MV> rhsMV = lp_->getCurrRHSVec();
    TEST_FOR_EXCEPTION((lhsMV==Teuchos::null && rhsMV==Teuchos::null),std::invalid_argument,
           "Belos::PseudoBlockCGIter::initialize(): Cannot initialize state storage!");

    // Get the multivector that is not null.
    Teuchos::RCP<const MV> tmp = ( (rhsMV!=Teuchos::null)? rhsMV: lhsMV );

    // Get the number of right-hand sides we're solving for now.
    int numRHS = MVT::GetNumberVecs(*tmp);
    numRHS_ = numRHS;
  
    // Initialize the state storage
    // If the subspace has not be initialized before or has changed sizes, generate it using the LHS or RHS from lp_.
    if (Teuchos::is_null(R_) || MVT::GetNumberVecs(*R_)!=numRHS_) {
      R_ = MVT::Clone( *tmp, numRHS_ );
      Z_ = MVT::Clone( *tmp, numRHS_ );
      P_ = MVT::Clone( *tmp, numRHS_ );
      AP_ = MVT::Clone( *tmp, numRHS_ );
    }
  
    // NOTE:  In CGIter R_, the initial residual, is required!!!  
    //
    std::string errstr("Belos::BlockPseudoCGIter::initialize(): Specified multivectors must have a consistent length and width.");

    // Create convenience variables for zero and one.
    const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
    const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();

    if (!Teuchos::is_null(newstate.R)) {

      TEST_FOR_EXCEPTION( MVT::GetVecLength(*newstate.R) != MVT::GetVecLength(*R_),
                          std::invalid_argument, errstr );
      TEST_FOR_EXCEPTION( MVT::GetNumberVecs(*newstate.R) != numRHS_,
                          std::invalid_argument, errstr );

      // Copy basis vectors from newstate into V
      if (newstate.R != R_) {
        // copy over the initial residual (unpreconditioned).
  MVT::MvAddMv( one, *newstate.R, zero, *newstate.R, *R_ );
      }

      // Compute initial direction vectors
      // Initially, they are set to the preconditioned residuals
      //
      if ( lp_->getLeftPrec() != Teuchos::null ) {
  lp_->applyLeftPrec( *R_, *Z_ ); 
  MVT::MvAddMv( one, *Z_, zero, *Z_, *P_ );
      } else {
  Z_ = R_;
  MVT::MvAddMv( one, *R_, zero, *R_, *P_ );
      }
    }
    else {

      TEST_FOR_EXCEPTION(Teuchos::is_null(newstate.R),std::invalid_argument,
                         "Belos::CGIter::initialize(): CGStateIterState does not have initial residual.");
    }

    // The solver is initialized
    initialized_ = true;
  }


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

    // Allocate memory for scalars.
    int i=0;
    std::vector<int> index(1);
    std::vector<ScalarType> rHz( numRHS_ ), rHz_old( numRHS_ ), pAp( numRHS_ );
    Teuchos::SerialDenseMatrix<int, ScalarType> alpha( numRHS_,numRHS_ ), beta( numRHS_,numRHS_ );

    // Create convenience variables for zero and one.
    const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
    const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
    
    // Get the current solution std::vector.
    Teuchos::RCP<MV> cur_soln_vec = lp_->getCurrLHSVec();

    // Compute first <r,z> a.k.a. rHz
    MVT::MvDot( *R_, *Z_, rHz );
    
    // Iterate until the status test tells us to stop.
    //
    while (stest_->checkStatus(this) != Passed) {
      
      // Increment the iteration
      iter_++;

      // Multiply the current direction std::vector by A and store in AP_
      lp_->applyOp( *P_, *AP_ );
      
      // Compute alpha := <R_,Z_> / <P_,AP_>
      MVT::MvDot( *P_, *AP_, pAp );
      for (i=0; i<numRHS_; ++i) {
        alpha(i,i) = rHz[i] / pAp[i];
        // Check that alpha is a positive number!
        TEST_FOR_EXCEPTION( SCT::real(alpha(i,i)) <= zero, CGIterateFailure,
        "Belos::PseudoBlockCGIter::iterate(): non-positive value for p^H*A*p encountered!" );
      }
      
      //
      // Update the solution std::vector x := x + alpha * P_
      //
      MVT::MvTimesMatAddMv( one, *P_, alpha, one, *cur_soln_vec );
      lp_->updateSolution();
      //
      // Save the denominator of beta before residual is updated [ old <R_, Z_> ]
      //
      for (i=0; i<numRHS_; ++i) {
        rHz_old[i] = rHz[i];
      }
      //
      // Compute the new residual R_ := R_ - alpha * AP_
      //
      MVT::MvTimesMatAddMv( -one, *AP_, alpha, one, *R_ );
      //
      // Compute beta := [ new <R_, Z_> ] / [ old <R_, Z_> ], 
      // and the new direction std::vector p.
      //
      if ( lp_->getLeftPrec() != Teuchos::null ) {
  lp_->applyLeftPrec( *R_, *Z_ );
      } else {
  Z_ = R_;
      }
      //
      MVT::MvDot( *R_, *Z_, rHz );
      //
      // Update the search directions.
      for (i=0; i<numRHS_; ++i) {
        beta(i,i) = rHz[i] / rHz_old[i];
  index[0] = i;
  Teuchos::RCP<MV> Z_i = MVT::CloneView( *Z_, index );
  Teuchos::RCP<MV> P_i = MVT::CloneView( *P_, index );
        MVT::MvAddMv( one, *Z_i, beta(i,i), *P_i, *P_i );
      }
      //      
    } // end while (sTest_->checkStatus(this) != Passed)
  }

} // end Belos namespace

#endif /* BELOS_PSEUDO_BLOCK_CG_ITER_HPP */

---