BelosGCRODRIter.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//                 Belos: Block Linear Solvers Package
//
// Copyright 2004-2016 NTESS and the Belos contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef BELOS_GCRODR_ITER_HPP
#define BELOS_GCRODR_ITER_HPP
 
#include "BelosConfigDefs.hpp"
#include "BelosTypes.hpp"
#include "BelosGCRODRIteration.hpp"
 
#include "BelosLinearProblem.hpp"
#include "BelosMatOrthoManager.hpp"
#include "BelosOutputManager.hpp"
#include "BelosStatusTest.hpp"
#include "BelosOperatorTraits.hpp"
#include "BelosMultiVecTraits.hpp"
#include "BelosDenseMatTraits.hpp"
 
#include "Teuchos_BLAS.hpp"
#include "Teuchos_ScalarTraits.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_TimeMonitor.hpp"
 
namespace Belos {
  
 
  
  template <class ScalarType, class MV, class DM>
  struct GCRODRIterState {
    int curDim;
    
    Teuchos::RCP<MV> V;
 
    Teuchos::RCP<MV> Z;
   
    Teuchos::RCP<MV> U, C;
 
    Teuchos::RCP<DM> H2;
 
    Teuchos::RCP<DM> H;
 
    Teuchos::RCP<DM> B;
 
    GCRODRIterState() : curDim(0), V(Teuchos::null), Z(Teuchos::null), 
      U(Teuchos::null), C(Teuchos::null),
      H2(Teuchos::null), H(Teuchos::null), 
      B(Teuchos::null)
    {}
  };
  
  
 
  
  class GCRODRIterInitFailure : public BelosError {
    public:
      GCRODRIterInitFailure(const std::string& what_arg) : BelosError(what_arg) {}
  };
  
  class GCRODRIterOrthoFailure : public BelosError {
    public:
      GCRODRIterOrthoFailure(const std::string& what_arg) : BelosError(what_arg) {}
  };
  
  
 
  template<class ScalarType, class MV, class OP, class DM>
  class GCRODRIter : virtual public GCRODRIteration<ScalarType,MV,OP,DM> {
    
  public:
    
    //
    // Convenience typedefs
    //
    typedef MultiVecTraits<ScalarType,MV,DM> MVT;
    typedef OperatorTraits<ScalarType,MV,OP> OPT;
    typedef DenseMatTraits<ScalarType,DM>    DMT;
    typedef Teuchos::ScalarTraits<ScalarType> SCT;
    typedef typename SCT::magnitudeType MagnitudeType;
    
 
    
    GCRODRIter( const Teuchos::RCP<LinearProblem<ScalarType,MV,OP,DM> > &problem,
    const Teuchos::RCP<OutputManager<ScalarType> > &printer,
    const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > &tester,
    const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP,DM> > &ortho,
    Teuchos::ParameterList &params );
    
    virtual ~GCRODRIter() {};
    
    
 
    
    void iterate();
    
    void initialize(GCRODRIterState<ScalarType,MV,DM>& newstate);
    
    void initialize() {
      GCRODRIterState<ScalarType,MV,DM> empty;
      initialize(empty);
    }
    
    GCRODRIterState<ScalarType,MV,DM> getState() const {
      GCRODRIterState<ScalarType,MV,DM> state;
      state.curDim = curDim_;
      state.V = V_;
      state.Z = Teuchos::null;
      state.U = U_;
      state.C = C_;
      state.H2 = H2_;
      state.H = H_;
      state.B = B_;
 
      return state;
    }
    
    
    
 
    
    int getNumIters() const { return iter_; }
    
    void resetNumIters( int iter = 0 ) { iter_ = iter; }
    
    Teuchos::RCP<const MV> getNativeResiduals( std::vector<MagnitudeType> *norms ) const;
    
 
    Teuchos::RCP<MV> getCurrentUpdate() const;
    
 
    void updateLSQR( int dim = -1 );
    
    int getCurSubspaceDim() const { 
      if (!initialized_) return 0;
      return curDim_;
    };
    
    int getMaxSubspaceDim() const { return numBlocks_; }
    
    
    
 
    
    const LinearProblem<ScalarType,MV,OP,DM>& getProblem() const { return *lp_; }
    
    int getNumBlocks() const { return numBlocks_; }
    
    void setNumBlocks(int numBlocks) { setSize( recycledBlocks_, numBlocks ); };
    
    int getBlockSize() const { return 1; }
    
    void setBlockSize(int blockSize) {
      TEUCHOS_TEST_FOR_EXCEPTION(blockSize!=1,std::invalid_argument,"Belos::GCRODRIter::setBlockSize(): Cannot use a block size that is not one.");
    }
 
    void setSize( int recycledBlocks, int numBlocks ) {
      // only call resize if size changed
      if ( recycledBlocks_ != recycledBlocks )
        recycledBlocks_ = recycledBlocks;
      if ( numBlocks_ != numBlocks ) {
        numBlocks_ = numBlocks;
        cs_.resize( numBlocks_+1 );
        sn_.resize( numBlocks_+1 );
        z_ = DMT::Create( numBlocks_+1, 1, false );
        R_ = DMT::Create( numBlocks_+1, numBlocks, false );
      }
    }
 
    bool isInitialized() { return initialized_; }
    
    
  private:
    
    //
    // Internal methods
    //
    
    //
    // Classes inputed through constructor that define the linear problem to be solved.
    //
    const Teuchos::RCP<LinearProblem<ScalarType,MV,OP,DM> >    lp_;
    const Teuchos::RCP<OutputManager<ScalarType> >          om_;
    const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> >       stest_;
    const Teuchos::RCP<OrthoManager<ScalarType,MV,DM> >        ortho_;
    
    //
    // Algorithmic parameters
    //  
    // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
    int numBlocks_; 
 
    // recycledBlocks_ is the size of the allocated space for the recycled subspace, in blocks.
    int recycledBlocks_; 
   
    // Storage for QR factorization of the least squares system.
    std::vector<ScalarType> sn_;
    std::vector<MagnitudeType> cs_;
 
    // 
    // 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 subspace dimension, and number of iterations performed.
    int curDim_, iter_;
   
    // Pointer to the (0,0) position of the upper Hessenberg matrix.
    int ptrH00_; 
    // 
    // State Storage
    //
    // Krylov vectors.
    Teuchos::RCP<MV> V_;
    //
    // Recycled subspace vectors.
    Teuchos::RCP<MV> U_, C_;
    //
    // Global projected matrix, including H_ and B_
    Teuchos::RCP<DM> H2_;
    //
    // Projected matrices
    // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
    Teuchos::RCP<DM> H_;
    //
    // B_ : Projected matrix from the recycled subspace B = C^H*A*V
    Teuchos::RCP<DM> B_;
    //
    // QR decomposition of Projected matrices for solving the least squares system HY = B.
    // R_: Upper triangular reduction of H
    // z_: Q applied to right-hand side of the least squares system
    Teuchos::RCP<DM> R_;
    Teuchos::RCP<DM> z_;  
  };
  
  // Constructor.
  template<class ScalarType, class MV, class OP, class DM>
  GCRODRIter<ScalarType,MV,OP,DM>::GCRODRIter(const Teuchos::RCP<LinearProblem<ScalarType,MV,OP,DM> > &problem,
             const Teuchos::RCP<OutputManager<ScalarType> > &printer,
             const Teuchos::RCP<StatusTest<ScalarType,MV,OP,DM> > &tester,
             const Teuchos::RCP<MatOrthoManager<ScalarType,MV,OP,DM> > &ortho,
             Teuchos::ParameterList &params ):
    lp_(problem), om_(printer), stest_(tester), ortho_(ortho) {
    numBlocks_      = 0;
    recycledBlocks_ = 0;
    initialized_    = false;
    curDim_         = 0;
    iter_           = 0;
    ptrH00_         = 0;
    V_              = Teuchos::null;
    U_              = Teuchos::null;
    C_              = Teuchos::null;
    H2_             = Teuchos::null;
    H_              = Teuchos::null;
    B_              = Teuchos::null;
 
    // Get the maximum number of blocks allowed for this Krylov subspace
    TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Num Blocks"), std::invalid_argument, "Belos::GCRODRIter::constructor: mandatory parameter \"Num Blocks\" is not specified.");
    int nb = Teuchos::getParameter<int>(params, "Num Blocks");
 
    TEUCHOS_TEST_FOR_EXCEPTION(!params.isParameter("Recycled Blocks"), std::invalid_argument,"Belos::GCRODRIter::constructor: mandatory parameter \"Recycled Blocks\" is not specified.");
    int rb = Teuchos::getParameter<int>(params, "Recycled Blocks");
 
    TEUCHOS_TEST_FOR_EXCEPTION(nb <= 0, std::invalid_argument, "Belos::GCRODRIter() was passed a non-positive argument for \"Num Blocks\".");
    TEUCHOS_TEST_FOR_EXCEPTION(rb >= nb, std::invalid_argument, "Belos::GCRODRIter() the number of recycled blocks is larger than the allowable subspace.");
 
    numBlocks_ = nb;
    recycledBlocks_ = rb;
    cs_.resize( numBlocks_+1 );
    sn_.resize( numBlocks_+1 );
    z_ = DMT::Create( numBlocks_+1, 1, false );
    R_ = DMT::Create( numBlocks_+1, numBlocks_, false );
 
  }
  
  // Get the current update from this subspace.
  template<class ScalarType, class MV, class OP, class DM>
  Teuchos::RCP<MV> GCRODRIter<ScalarType,MV,OP,DM>::getCurrentUpdate() const {
    //
    // If this is the first iteration of the Arnoldi factorization, 
    // there is no update, so return Teuchos::null. 
    //
    Teuchos::RCP<MV> currentUpdate = Teuchos::null;
    if (curDim_==0) { 
      return currentUpdate; 
    } else {
      const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
      const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
      Teuchos::BLAS<int,ScalarType> blas;
      currentUpdate = MVT::Clone( *V_, 1 );
      //
      //  Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
      //
      Teuchos::RCP<DM> y = DMT::SubviewCopy(*z_, curDim_, 1);
      DMT::SyncDeviceToHost( *y );
      DMT::SyncDeviceToHost( *R_ );
      //
      //  Solve the least squares problem.
      //
      blas.TRSM( Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
     Teuchos::NON_UNIT_DIAG, curDim_, 1, one,  
     DMT::GetConstRawHostPtr(*R_), DMT::GetStride(*R_), 
     DMT::GetRawHostPtr(*y), DMT::GetStride(*y) );
      //
      //  Compute the current update from the Krylov basis; V(:,1:curDim_)*y.
      //
      std::vector<int> index(curDim_);
      for ( int i=0; i<curDim_; i++ ) index[i] = i;
      Teuchos::RCP<const MV> Vjp1 = MVT::CloneView( *V_, index );
      MVT::MvTimesMatAddMv( one, *Vjp1, *y, zero, *currentUpdate );
      //
      //  Add in portion of update from recycled subspace U; U(:,1:recycledBlocks_)*B*y.
      //
      if (U_ != Teuchos::null) {
  Teuchos::RCP<DM> z = DMT::Create( recycledBlocks_, 1 );
  DMT::SyncDeviceToHost( *H2_ );
        blas.GEMM( Teuchos::NO_TRANS, Teuchos::NO_TRANS, recycledBlocks_, 1, curDim_, one, 
       DMT::GetConstRawHostPtr(*B_), DMT::GetStride(*B_),
       DMT::GetConstRawHostPtr(*y), DMT::GetStride(*y),
       zero, DMT::GetRawHostPtr(*z), DMT::GetStride(*z));
        DMT::SyncHostToDevice( *z );
        MVT::MvTimesMatAddMv( -one, *U_, *z, one, *currentUpdate );
      }
    }
    std::vector<MagnitudeType> normUpdate( 1 );
    MVT::MvNorm( *currentUpdate, normUpdate );
 
    return currentUpdate;
  }
 
 
  // Get the native residuals stored in this iteration.  
  // Note: This method does not return a MultiVector of the residual vectors, only return Teuchos::null
  template<class ScalarType, class MV, class OP, class DM>
  Teuchos::RCP<const MV> GCRODRIter<ScalarType,MV,OP,DM>::getNativeResiduals( std::vector<MagnitudeType> *norms ) const {
    //
    // NOTE: Make sure the incoming std::vector is the correct size!
    //
    if ( norms && (int)norms->size()==0 )                         
      norms->resize( 1 );                                          
    
    if (norms) {
      DMT::SyncDeviceToHost( *z_ );
      const ScalarType curNativeResid = DMT::ValueConst(*z_,curDim_,0);
      (*norms)[0] = SCT::magnitude (curNativeResid);
    }
    return Teuchos::null;
  }
  
  
 
  // Initialize this iteration object
  template<class ScalarType, class MV, class OP, class DM>
  void GCRODRIter<ScalarType,MV,OP,DM>::initialize(GCRODRIterState<ScalarType,MV,DM>& newstate) {
   
    if (newstate.V != Teuchos::null && newstate.H2 != Teuchos::null) {
      curDim_ = newstate.curDim;
      V_      = newstate.V;
      U_      = newstate.U;
      C_      = newstate.C;
      H2_     = newstate.H2;
      // There is no recycled space, so this iteration is creating the first one.
      if (newstate.U == Teuchos::null) {
        ptrH00_ = recycledBlocks_+1;
        H_ = DMT::Subview( *H2_, numBlocks_+1, numBlocks_, ptrH00_, ptrH00_ );
        B_ = Teuchos::null;
      }
      else {        
        ptrH00_ = recycledBlocks_;
        H_ = DMT::Subview( *H2_, numBlocks_+1, numBlocks_, ptrH00_, ptrH00_ );
        B_ = DMT::Subview( *H2_, recycledBlocks_, numBlocks_, 0, ptrH00_ );
      }
    }
    else {
      TEUCHOS_TEST_FOR_EXCEPTION(newstate.V == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): GCRODRIterState does not have V initialized.");
      TEUCHOS_TEST_FOR_EXCEPTION(newstate.H2 == Teuchos::null,std::invalid_argument,"Belos::GCRODRIter::initialize(): GCRODRIterState does not have H2 initialized.");
    }
 
    // the solver is initialized
    initialized_ = true;
 
  }
 
 
  // Iterate until the status test informs us we should stop.
  template<class ScalarType, class MV, class OP, class DM>
  void GCRODRIter<ScalarType,MV,OP,DM>::iterate() {
 
    TEUCHOS_TEST_FOR_EXCEPTION( initialized_ == false, GCRODRIterInitFailure,"Belos::GCRODRIter::iterate(): GCRODRIter class not initialized." );
 
    // Force call to setsize to ensure internal storage is correct dimension
    setSize( recycledBlocks_, numBlocks_ );
    
    Teuchos::RCP<MV> Vnext;
    Teuchos::RCP<const MV> Vprev;
    std::vector<int> curind(1);
 
    // z_ must be zeroed out in order to compute Givens rotations correctly
    DMT::PutScalar(*z_);
 
    // Orthonormalize the new V_0
    curind[0] = 0;
    Vnext = MVT::CloneViewNonConst(*V_,curind);
    // Orthonormalize first column
    // First, get a view of the first element of z_ to hold the orthonormalization coefficients
    Teuchos::RCP<DM> z0 = DMT::Subview( *z_, 1, 1 );
    int rank = ortho_->normalize( *Vnext, z0 );
    TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
 
    std::vector<int> prevind(numBlocks_+1);
 
    // iterate until the status test tells us to stop.
    //
    // also break if our basis is full
    //
    if (U_ == Teuchos::null) { // iterate without recycle space (because none is available yet)
      while (stest_->checkStatus(this) != Passed && curDim_+1 <= numBlocks_) {
        iter_++;
        int lclDim = curDim_ + 1;
 
        // Get next index in basis
        curind[0] = lclDim;
        Vnext = MVT::CloneViewNonConst(*V_,curind);
 
        // Get previous index in basis
        curind[0] = curDim_;
        Vprev = MVT::CloneView(*V_,curind);
 
        // Compute the next vector in the Krylov basis:  Vnext = Op*Vprev
        lp_->apply(*Vprev,*Vnext);
 
        // Remove all previous Krylov basis vectors from Vnext and put coefficients in H and R.
 
        // Get a view of all the previous vectors
        prevind.resize(lclDim);
        for (int i=0; i<lclDim; i++) { prevind[i] = i; }
        Vprev = MVT::CloneView(*V_,prevind);
        Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
 
        // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
  Teuchos::RCP<DM> subH = DMT::Subview(*H2_, lclDim, 1, ptrH00_, ptrH00_+curDim_);
        Teuchos::Array<Teuchos::RCP<DM> > AsubH( 1, subH );
 
        // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
        Teuchos::RCP<DM> subR = DMT::Subview(*H2_, 1, 1, ptrH00_+lclDim, ptrH00_+curDim_);
 
        // Project out the previous Krylov vectors and normalize the next vector.
  rank = ortho_->projectAndNormalize(*Vnext, AsubH, subR, AVprev);
 
  // Copy over the coefficients to R just in case we run into an error.
  Teuchos::RCP<DM> subR2 = DMT::Subview(*R_, lclDim+1, 1, 0, curDim_);
        Teuchos::RCP<const DM> subH2 = DMT::SubviewConst(*H2_, lclDim+1, 1, ptrH00_, ptrH00_+curDim_);
  DMT::Assign(*subR2, *subH2);
  subR2 = Teuchos::null;
 
        TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
 
        // Update the QR factorization of the upper Hessenberg matrix
        updateLSQR();
 
        // Update basis dim
        curDim_++;
      } // end while (statusTest == false)
    }
    else { // iterate with recycle space
      while (stest_->checkStatus(this) != Passed && curDim_+1 <= numBlocks_) {
        iter_++;
        int lclDim = curDim_ + 1; 
 
        // Get the current part of the basis.
        curind[0] = lclDim;
        Vnext = MVT::CloneViewNonConst(*V_,curind);
 
        // Get a view of the previous vectors.
        // This is used for orthogonalization and for computing V^H K H
        curind[0] = curDim_;
        Vprev = MVT::CloneView(*V_,curind);
 
        // Compute the next std::vector in the Krylov basis:  Vnext = Op*Vprev
        lp_->apply(*Vprev,*Vnext);
        Vprev = Teuchos::null;
 
  DMT::SyncHostToDevice(*H2_);
 
        // First, remove the recycled subspace (C) from Vnext and put coefficients in B.
        Teuchos::Array<Teuchos::RCP<const MV> > C(1, C_);
        Teuchos::RCP<DM> subB = DMT::Subview(*H2_, recycledBlocks_, 1, 0, ptrH00_+curDim_);
        Teuchos::RCP<DM> tmpB = DMT::Create(recycledBlocks_, 1);
  Teuchos::Array<Teuchos::RCP<DM> > AsubB( 1, subB );
 
        // Project out the recycled subspace.
        ortho_->project( *Vnext, AsubB, C );
 
        // Now, remove all previous Krylov basis vectors from Vnext and put coefficients in H and R.          
        // Get a view of all the previous vectors
        prevind.resize(lclDim);
        for (int i=0; i<lclDim; i++) { prevind[i] = i; }
        Vprev = MVT::CloneView(*V_,prevind);
        Teuchos::Array<Teuchos::RCP<const MV> > AVprev(1, Vprev);
  
        // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
        Teuchos::RCP<DM> subH = DMT::Subview(*H2_, lclDim, 1, ptrH00_, ptrH00_+curDim_);
        Teuchos::Array<Teuchos::RCP<DM> > AsubH(1, subH);
      
        // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
        Teuchos::RCP<DM> subR = DMT::Subview(*H2_, 1, 1, ptrH00_+lclDim, ptrH00_+curDim_);
 
        // Project out the previous Krylov vectors and normalize the next vector.
        rank = ortho_->projectAndNormalize(*Vnext, AsubH, subR, AVprev);
 
        // Copy over the coefficients to R just in case we run into an error.
  Teuchos::RCP<DM> subR2 = DMT::Subview(*R_, lclDim+1, 1, 0, curDim_);
        Teuchos::RCP<const DM> subH2 = DMT::SubviewConst(*H2_, lclDim+1, 1, ptrH00_, ptrH00_+curDim_);
  DMT::Assign(*subR2, *subH2);
  subR2 = Teuchos::null;
      
        TEUCHOS_TEST_FOR_EXCEPTION(rank != 1,GCRODRIterOrthoFailure, "Belos::GCRODRIter::iterate(): couldn't generate basis of full rank.");
 
        // Update the QR factorization of the upper Hessenberg matrix
        updateLSQR();
 
        // Update basis dim
        curDim_++;
      
      } // end while (statusTest == false)
    } // end if (U_ == Teuchos::null)
   
  }
 
  
  template<class ScalarType, class MV, class OP, class DM>
  void GCRODRIter<ScalarType,MV,OP,DM>::updateLSQR( int dim ) {
 
    int i;
    const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero();
   
    // Get correct dimension based on input "dim"
    // Remember that ortho failures result in an exit before updateLSQR() is called.
    // Therefore, it is possible that dim == curDim_.
    int curDim = curDim_;
    if ( (dim >= curDim_) && (dim < getMaxSubspaceDim()) )
      curDim = dim;
    
    Teuchos::BLAS<int, ScalarType> blas;
    //
    // Apply previous transformations and compute new transformation to reduce upper-Hessenberg
    // system to upper-triangular form.
    //
    // QR factorization of Least-Squares system with Givens rotations
    //
    DMT::SyncDeviceToHost(*R_);
    DMT::SyncDeviceToHost(*z_);
    //
    for (i=0; i<curDim; i++) {
      //
      // Apply previous Givens rotations to new column of Hessenberg matrix
      //
      blas.ROT( 1, &(DMT::Value(*R_,i,curDim)), 1, &(DMT::Value(*R_,i+1, curDim)), 1, &cs_[i], &sn_[i] );
 
    }
    //
    // Calculate new Givens rotation
    //
    blas.ROTG( &(DMT::Value(*R_,curDim,curDim)), &(DMT::Value(*R_,curDim+1,curDim)), &cs_[curDim], &sn_[curDim] );
    DMT::Value(*R_,curDim+1,curDim) = zero;
    //
    // Update RHS w/ new transformation
    //
    blas.ROT( 1, &(DMT::Value(*z_,curDim,0)), 1, &(DMT::Value(*z_,curDim+1,0)), 1, &cs_[curDim], &sn_[curDim] );
    //
    DMT::SyncHostToDevice(*R_);
    DMT::SyncHostToDevice(*z_);
    //
  } // end updateLSQR()
 
} // end Belos namespace
 
#endif /* BELOS_GCRODR_ITER_HPP */