BelosBlockFGmresIter.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_BLOCK_FGMRES_ITER_HPP
#define BELOS_BLOCK_FGMRES_ITER_HPP

 
#include "BelosConfigDefs.hpp"
#include "BelosTypes.hpp"
#include "BelosGmresIteration.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"
 
#include <vector>

 
namespace Belos {
 
template<class ScalarType, class MV, class OP, class DM = DefaultDenseMatrix<int,ScalarType>>
class BlockFGmresIter : virtual public GmresIteration<ScalarType,MV,OP,DM> {
 
  public:
 
  //
  // Convenience typedefs
  //
  typedef MultiVecTraits<ScalarType,MV,DM> MVT;
  typedef DenseMatTraits<ScalarType,DM> DMT;
  typedef OperatorTraits<ScalarType,MV,OP> OPT;
  typedef Teuchos::ScalarTraits<ScalarType> SCT;
  typedef typename SCT::magnitudeType MagnitudeType;



  BlockFGmresIter( 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 ~BlockFGmresIter() {};
 



  void iterate();

  void initializeGmres(GmresIterationState<ScalarType,MV,DM>& newstate);

  void initialize()
  {
    GmresIterationState<ScalarType,MV,DM> empty;
    initializeGmres(empty);
  }

  GmresIterationState<ScalarType,MV,DM> getState() const {
    GmresIterationState<ScalarType,MV,DM> state;
    state.curDim = curDim_;
    state.V = V_;
    state.Z = Z_;
    state.H = H_;
    state.R = R_;
    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;


  Teuchos::RCP<MV> getCurrentUpdate() const;


  void updateLSQR( int dim = -1 );

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

  int getMaxSubspaceDim() const { return blockSize_*numBlocks_; }

 



  const LinearProblem<ScalarType,MV,OP,DM>& getProblem() const { return *lp_; }

  int getBlockSize() const { return blockSize_; }

  void setBlockSize(int blockSize) { setSize( blockSize, numBlocks_ ); }

  int getNumBlocks() const { return numBlocks_; }

  void setNumBlocks(int numBlocks) { setSize( blockSize_, numBlocks ); }

  void setSize(int blockSize, int numBlocks);

  bool isInitialized() { return initialized_; }

 
  private:
 
  //
  // Internal methods
  //
  void setStateSize();
 
  //
  // 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
  //
  // blockSize_ is the solver block size.
  // It controls the number of vectors added to the basis on each iteration.
  int blockSize_;
  // numBlocks_ is the size of the allocated space for the Krylov basis, in blocks.
  int numBlocks_;
 
  // Storage for QR factorization of the least squares system.
  std::vector<ScalarType> beta, 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_;
 
  // stateStorageInitialized_ specified that the state storage has be initialized to the current
  // blockSize_ and numBlocks_.  This initialization may be postponed if the linear problem was
  // generated without the right-hand side or solution vectors.
  bool stateStorageInitialized_;
 
  // keepHessenberg_ specifies that the upper Hessenberg matrix should be stored separately
  // from the QR-factored least squares system (R_).  When false, H_ and R_ point to the
  // same object and only the QR-rotated form is available via getState().
  bool keepHessenberg_;
 
  // Current subspace dimension, and number of iterations performed.
  int curDim_, iter_;
 
  //
  // State Storage
  //
  Teuchos::RCP<MV> V_;
  Teuchos::RCP<MV> Z_;
  //
  // Projected matrices
  // H_ : Projected matrix from the Krylov factorization AV = VH + FE^T
  //
  Teuchos::RCP<DM> H_;
  //
  // 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>
  BlockFGmresIter<ScalarType,MV,OP,DM>::
  BlockFGmresIter (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),
    blockSize_(0),
    numBlocks_(0),
    initialized_(false),
    stateStorageInitialized_(false),
    keepHessenberg_(false),
    curDim_(0),
    iter_(0)
  {
    // Find out whether we are saving the Hessenberg matrix separately from R.
    if (om_->isVerbosity(Debug))
      keepHessenberg_ = true;
    else
      keepHessenberg_ = params.get("Keep Hessenberg", false);
 
    // Get the maximum number of blocks allowed for this Krylov subspace
    TEUCHOS_TEST_FOR_EXCEPTION(
      ! params.isParameter ("Num Blocks"), std::invalid_argument,
      "Belos::BlockFGmresIter::constructor: mandatory parameter 'Num Blocks' is not specified.");
    const int nb = params.get<int> ("Num Blocks");
 
    // Set the block size and allocate data.
    const int bs = params.get ("Block Size", 1);
    setSize (bs, nb);
  }

  // Set the block size and make necessary adjustments.
  template <class ScalarType, class MV, class OP, class DM>
  void BlockFGmresIter<ScalarType,MV,OP,DM>::setSize (int blockSize, int numBlocks)
  {
    // This routine only allocates space; it doesn't not perform any computation
    // any change in size will invalidate the state of the solver.
 
    TEUCHOS_TEST_FOR_EXCEPTION(numBlocks <= 0 || blockSize <= 0, std::invalid_argument, "Belos::BlockFGmresIter::setSize was passed a non-positive argument.");
    if (blockSize == blockSize_ && numBlocks == numBlocks_) {
      // do nothing
      return;
    }
 
    if (blockSize!=blockSize_ || numBlocks!=numBlocks_)
      stateStorageInitialized_ = false;
 
    blockSize_ = blockSize;
    numBlocks_ = numBlocks;
 
    initialized_ = false;
    curDim_ = 0;
 
    // Use the current blockSize_ and numBlocks_ to initialize the state storage.
    setStateSize();
 
  }

  // Setup the state storage.
  template <class ScalarType, class MV, class OP, class DM>
  void BlockFGmresIter<ScalarType,MV,OP,DM>::setStateSize ()
  {
    using Teuchos::RCP;
    using Teuchos::rcp;
 
    if (! stateStorageInitialized_) {
      // Check if there is any multivector to clone from.
      RCP<const MV> lhsMV = lp_->getLHS();
      RCP<const MV> rhsMV = lp_->getRHS();
      if (lhsMV == Teuchos::null && rhsMV == Teuchos::null) {
        stateStorageInitialized_ = false;
        return;
      }
      else {
        // blockSize*numBlocks dependent
        //
        int newsd = blockSize_*(numBlocks_+1);
 
        if (blockSize_==1) {
          cs.resize (newsd);
          sn.resize (newsd);
        }
        else {
          beta.resize (newsd);
        }
 
        // Initialize the state storage
        TEUCHOS_TEST_FOR_EXCEPTION(
          blockSize_ * static_cast<ptrdiff_t> (numBlocks_) > MVT::GetGlobalLength (*rhsMV),
          std::invalid_argument, "Belos::BlockFGmresIter::setStateSize(): "
          "Cannot generate a Krylov basis with dimension larger the operator!");
 
        // If the subspace has not be initialized before, generate it using the LHS or RHS from lp_.
        if (V_ == Teuchos::null) {
          // Get the multivector that is not null.
          RCP<const MV> tmp = (rhsMV != Teuchos::null) ? rhsMV : lhsMV;
          TEUCHOS_TEST_FOR_EXCEPTION(
            tmp == Teuchos::null, std::invalid_argument,
            "Belos::BlockFGmresIter::setStateSize(): "
            "linear problem does not specify multivectors to clone from.");
          V_ = MVT::Clone (*tmp, newsd);
        }
        else {
          // Generate V_ by cloning itself ONLY if more space is needed.
          if (MVT::GetNumberVecs (*V_) < newsd) {
            RCP<const MV> tmp = V_;
            V_ = MVT::Clone (*tmp, newsd);
          }
        }
 
        if (Z_ == Teuchos::null) {
          // Get the multivector that is not null.
          RCP<const MV> tmp = (rhsMV != Teuchos::null) ? rhsMV : lhsMV;
          TEUCHOS_TEST_FOR_EXCEPTION(
            tmp == Teuchos::null, std::invalid_argument,
            "Belos::BlockFGmresIter::setStateSize(): "
            "linear problem does not specify multivectors to clone from.");
          Z_ = MVT::Clone (*tmp, newsd);
        }
        else {
          // Generate Z_ by cloning itself ONLY if more space is needed.
          if (MVT::GetNumberVecs (*Z_) < newsd) {
            RCP<const MV> tmp = Z_;
            Z_ = MVT::Clone (*tmp, newsd);
          }
        }
 
        // R_ holds the QR-factored least squares system. Always allocated.
        if (R_ == Teuchos::null) {
          R_ = DMT::Create(newsd, newsd-blockSize_);
        }
        else {
          DMT::Reshape(*R_, newsd, newsd - blockSize_);
        }
 
        // H_ holds the raw (pre-QR) upper Hessenberg matrix.
        // When keepHessenberg_ is false, H_ and R_ point to the same object
        // (matching BlockGmresIter behavior).
        if (keepHessenberg_) {
          if (H_ == Teuchos::null) {
            H_ = DMT::Create(newsd, newsd-blockSize_);
          }
          else {
            DMT::Reshape(*H_, newsd, newsd - blockSize_);
          }
        }
        else {
          H_ = R_;
        }
 
        // Generate z_ only if it doesn't exist, otherwise resize it.
        if (z_ == Teuchos::null) {
          z_ = DMT::Create(newsd, blockSize_);
        }
        else {
          DMT::Reshape(*z_, newsd, blockSize_);
        }
 
        // State storage has now been initialized.
        stateStorageInitialized_ = true;
      }
    }
  }
 
 
  template <class ScalarType, class MV, class OP, class DM>
  Teuchos::RCP<MV>
  BlockFGmresIter<ScalarType,MV,OP,DM>::getCurrentUpdate() const
  {
    Teuchos::RCP<MV> currentUpdate = Teuchos::null;
    if (curDim_ == 0) {
      // If this is the first iteration of the Arnoldi factorization,
      // then there is no update, so return Teuchos::null.
      return currentUpdate;
    }
    else {
      const ScalarType zero = Teuchos::ScalarTraits<ScalarType>::zero ();
      const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one ();
      Teuchos::BLAS<int,ScalarType> blas;
 
      currentUpdate = MVT::Clone (*Z_, blockSize_);
 
      // Make a view and then copy the RHS of the least squares problem.  DON'T OVERWRITE IT!
      DMT::SyncDeviceToHost( *z_ );
      DMT::SyncDeviceToHost( *H_ );
 
      Teuchos::RCP<DM> y = DMT::SubviewCopy(*z_, curDim_, blockSize_);
 
      // Solve the least squares problem.
      blas.TRSM (Teuchos::LEFT_SIDE, Teuchos::UPPER_TRI, Teuchos::NO_TRANS,
                 Teuchos::NON_UNIT_DIAG, curDim_, blockSize_, one,
                 DMT::GetConstRawHostPtr(*R_), DMT::GetStride(*R_), 
                 DMT::GetRawHostPtr(*y), DMT::GetStride(*y));
 
      // Make sure the result goes back to the device
      DMT::SyncHostToDevice( *y );
 
      // Compute the current update.
      std::vector<int> index (curDim_);
      for (int i = 0; i < curDim_; ++i) {
        index[i] = i;
      }
      Teuchos::RCP<const MV> Zjp1 = MVT::CloneView (*Z_, index);
      MVT::MvTimesMatAddMv (one, *Zjp1, *y, zero, *currentUpdate);
    }
    return currentUpdate;
  }
 
 
  template <class ScalarType, class MV, class OP, class DM>
  Teuchos::RCP<const MV>
  BlockFGmresIter<ScalarType,MV,OP,DM>::
  getNativeResiduals (std::vector<MagnitudeType> *norms) const
  {
    // NOTE: Make sure the incoming std::vector is the correct size!
    if (norms != NULL && (int)norms->size() < blockSize_) {
      norms->resize (blockSize_);
    }
 
    if (norms != NULL) {
      Teuchos::BLAS<int, ScalarType> blas;
      DMT::SyncDeviceToHost( *z_ );
      for (int j = 0; j < blockSize_; ++j) {
        (*norms)[j] = blas.NRM2 (blockSize_, &DMT::Value(*z_, curDim_, j), 1);
      }
    }
 
    // FGmres does not return a residual (multi)vector.
    return Teuchos::null;
  }
 
 
  template <class ScalarType, class MV, class OP, class DM>
  void BlockFGmresIter<ScalarType,MV,OP,DM>::
  initializeGmres (GmresIterationState<ScalarType,MV,DM>& newstate)
  {
    using Teuchos::RCP;
    using Teuchos::rcp;
    using std::endl;
 
    // Initialize the state storage if it isn't already.
    if (! stateStorageInitialized_) {
      setStateSize ();
    }
 
    TEUCHOS_TEST_FOR_EXCEPTION(
      ! stateStorageInitialized_, std::invalid_argument,
      "Belos::BlockFGmresIter::initialize(): Cannot initialize state storage!");
 
    // NOTE: In BlockFGmresIter, V and Z are required!!!  Inconsistent
    // multivectors widths and lengths will not be tolerated, and will
    // be treated with exceptions.
    const char errstr[] = "Belos::BlockFGmresIter::initialize(): The given "
      "multivectors must have a consistent length and width.";
 
    if (! newstate.V.is_null () && ! newstate.z.is_null ()) {
 
      // initialize V_,z_, and curDim_
 
      TEUCHOS_TEST_FOR_EXCEPTION(
        MVT::GetGlobalLength(*newstate.V) != MVT::GetGlobalLength(*V_),
        std::invalid_argument, errstr );
      TEUCHOS_TEST_FOR_EXCEPTION(
        MVT::GetNumberVecs(*newstate.V) < blockSize_,
        std::invalid_argument, errstr );
      TEUCHOS_TEST_FOR_EXCEPTION(
        newstate.curDim > blockSize_*(numBlocks_+1),
        std::invalid_argument, errstr );
 
      curDim_ = newstate.curDim;
      const int lclDim = MVT::GetNumberVecs(*newstate.V);
 
      // check size of Z
      TEUCHOS_TEST_FOR_EXCEPTION(
        DMT::GetNumRows(*newstate.z) < curDim_ || DMT::GetNumCols(*newstate.z) < blockSize_,
        std::invalid_argument, errstr);
 
      // copy basis vectors from newstate into V
      if (newstate.V != V_) {
        // only copy over the first block and print a warning.
        if (curDim_ == 0 && lclDim > blockSize_) {
          std::ostream& warn = om_->stream (Warnings);
          warn << "Belos::BlockFGmresIter::initialize(): the solver was "
               << "initialized with a kernel of " << lclDim << endl
               << "The block size however is only " << blockSize_ << endl
               << "The last " << lclDim - blockSize_
               << " vectors will be discarded." << endl;
        }
        std::vector<int> nevind (curDim_ + blockSize_);
        for (int i = 0; i < curDim_ + blockSize_; ++i) {
          nevind[i] = i;
        }
        RCP<const MV> newV = MVT::CloneView (*newstate.V, nevind);
        RCP<MV> lclV = MVT::CloneViewNonConst (*V_, nevind);
        MVT::Assign(*newV, *lclV);
 
        // done with local pointers
        lclV = Teuchos::null;
      }
 
      // put data into z_, make sure old information is not still hanging around.
      if (newstate.z != z_) {
        DMT::PutScalar(*z_);
        RCP<const DM> newZ = DMT::SubviewConst(*newstate.z, curDim_ + blockSize_, blockSize_);
        RCP<DM> lclz = DMT::Subview(*z_, curDim_ + blockSize_, blockSize_);
        DMT::Assign(*lclz, *newZ);
        lclz = Teuchos::null; // done with local pointers
      }
    }
    else {
      TEUCHOS_TEST_FOR_EXCEPTION(
        newstate.V == Teuchos::null,std::invalid_argument,
        "Belos::BlockFGmresIter::initialize(): BlockFGmresStateIterState does not have initial kernel V_0.");
 
      TEUCHOS_TEST_FOR_EXCEPTION(
        newstate.z == Teuchos::null,std::invalid_argument,
        "Belos::BlockFGmresIter::initialize(): BlockFGmresStateIterState does not have initial norms z_0.");
    }
 
    // the solver is initialized
    initialized_ = true;
  }
 
 
  template <class ScalarType, class MV, class OP, class DM>
  void BlockFGmresIter<ScalarType,MV,OP,DM>::iterate()
  {
    using Teuchos::RCP;
    using Teuchos::rcp;
 
    // Allocate/initialize data structures
    if (initialized_ == false) {
      initialize();
    }
 
    // Compute the current search dimension.
    const int searchDim = blockSize_ * numBlocks_;
 
    // Iterate until the status test tells us to stop.
    // Raise an exception if a computed block is not full rank.
    while (stest_->checkStatus (this) != Passed && curDim_+blockSize_ <= searchDim) {
      ++iter_;
 
      // F can be found at the curDim_ block, but the next block is at curDim_ + blockSize_.
      const int lclDim = curDim_ + blockSize_;
 
      // Get the current part of the basis.
      std::vector<int> curind (blockSize_);
      for (int i = 0; i < blockSize_; ++i) {
        curind[i] = lclDim + i;
      }
      RCP<MV> Vnext = MVT::CloneViewNonConst (*V_, curind);
 
      // Get a view of the previous vectors.
      // This is used for orthogonalization and for computing V^H K H.
      for (int i = 0; i < blockSize_; ++i) {
        curind[i] = curDim_ + i;
      }
      RCP<const MV> Vprev = MVT::CloneView (*V_, curind);
      RCP<MV> Znext = MVT::CloneViewNonConst (*Z_, curind);
 
      // Compute the next (multi)vector in the Krylov basis:  Znext = M*Vprev
      lp_->applyRightPrec (*Vprev, *Znext);
      Vprev = Teuchos::null;
 
      // Compute the next (multi)vector in the Krylov basis:  Vnext = A*Znext
      lp_->applyOp (*Znext, *Vnext);
      Znext = Teuchos::null;
 
      // Remove all previous Krylov basis vectors from Vnext
      // Get a view of all the previous vectors
      std::vector<int> prevind (lclDim);
      for (int i = 0; i < lclDim; ++i) {
        prevind[i] = i;
      }
      Vprev = MVT::CloneView (*V_, prevind);
      Teuchos::Array<RCP<const MV> > AVprev (1, Vprev);
 
      // Get a view of the part of the Hessenberg matrix needed to hold the ortho coeffs.
      // Ortho always writes into H_ (the raw Hessenberg).
      RCP<DM> subH = DMT::Subview(*H_, lclDim, blockSize_, 0, curDim_);
      Teuchos::Array<RCP<DM> > AsubH;
      AsubH.append (subH);
 
      // Get a view of the part of the Hessenberg matrix needed to hold the norm coeffs.
      RCP<DM> subH2 = DMT::Subview(*H_, blockSize_, blockSize_, lclDim, curDim_);
      const int rank = ortho_->projectAndNormalize (*Vnext, AsubH, subH2, AVprev);
      TEUCHOS_TEST_FOR_EXCEPTION(
        rank != blockSize_, GmresIterationOrthoFailure,
        "Belos::BlockFGmresIter::iterate(): After orthogonalization, the new "
        "basis block does not have full rank.  It contains " << blockSize_
        << " vector" << (blockSize_ != 1 ? "s" : "")
        << ", but its rank is " << rank << ".");
 
      // If keeping the Hessenberg separately, copy the new columns into R_
      // before updateLSQR() overwrites them with the QR factorization.
      if (keepHessenberg_) {
        // Copy over the orthogonalization coefficients.
        RCP<DM> subR = DMT::Subview(*R_,lclDim,blockSize_,0,curDim_ );
        DMT::Assign(*subR,*subH);
 
        // Copy over the lower diagonal block of the Hessenberg matrix.
        RCP<DM> subR2 = DMT::Subview(*R_,blockSize_,blockSize_,lclDim,curDim_ );
        DMT::Assign(*subR2,*subH2);
      }
 
      //
      // V has been extended, and H has been extended.
      //
      // Update the QR factorization of the upper Hessenberg matrix (applied to R_).
      //
      updateLSQR ();
      //
      // Update basis dim and release all pointers.
      //
      Vnext = Teuchos::null;
      curDim_ += blockSize_;
    } // end while (statusTest == false)
  }
 
 
  template<class ScalarType, class MV, class OP, class DM>
  void BlockFGmresIter<ScalarType,MV,OP,DM>::updateLSQR (int dim)
  {
    typedef Teuchos::ScalarTraits<ScalarType> STS;
    typedef Teuchos::ScalarTraits<MagnitudeType> STM;
 
    const ScalarType zero = STS::zero ();
    const ScalarType two = (STS::one () + STS::one());
    ScalarType sigma, mu, vscale, maxelem;
    Teuchos::BLAS<int, ScalarType> blas;
 
    // Get correct dimension based on input 'dim'.  Remember that
    // orthogonalization 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;
    }
 
    // Apply previous transformations, and compute new transformation
    // to reduce upper Hessenberg system to upper triangular form.
    // The type of transformation we use depends the block size.  We
    // use Givens rotations for a block size of 1, and Householder
    // reflectors otherwise.
    DMT::SyncDeviceToHost( *H_ );
    DMT::SyncDeviceToHost( *z_ );
 
    if (blockSize_ == 1) {
 
      // QR factorization of upper Hessenberg matrix using Givens rotations
      for (int 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]);
    }
    else {
      // QR factorization of least-squares system using Householder reflectors.
      for (int j = 0; j < blockSize_; ++j) {
        // Apply previous Householder reflectors to new block of Hessenberg matrix
        for (int i = 0; i < curDim + j; ++i) {
          sigma = blas.DOT (blockSize_, &DMT::Value(*R_,i+1,i), 1, &DMT::Value(*R_,i+1,curDim+j), 1);
          sigma += DMT::ValueConst(*R_,i,curDim+j);
          sigma *= beta[i];
          blas.AXPY (blockSize_, ScalarType(-sigma), &DMT::Value(*R_,i+1,i), 1, &DMT::Value(*R_,i+1,curDim+j), 1);
          DMT::Value(*R_,i,curDim+j) -= sigma;
        }
 
        // Compute new Householder reflector
        const int maxidx = blas.IAMAX (blockSize_+1, &DMT::Value(*R_,curDim+j,curDim+j), 1);
        maxelem = DMT::ValueConst(*R_,curDim + j + maxidx - 1, curDim + j);
        for (int i = 0; i < blockSize_ + 1; ++i) {
          DMT::Value(*R_,curDim+j+i,curDim+j) /= maxelem;
        }
        sigma = blas.DOT (blockSize_, &DMT::Value(*R_,curDim + j + 1, curDim + j), 1,
                          &DMT::Value(*R_,curDim + j + 1, curDim + j), 1);
        if (sigma == zero) {
          beta[curDim + j] = zero;
        } else {
          mu = STS::squareroot (DMT::Value(*R_,curDim+j,curDim+j)*DMT::Value(*R_,curDim+j,curDim+j)+sigma);
          if (STS::real (DMT::Value(*R_,curDim + j, curDim + j)) < STM::zero ()) {
            vscale = DMT::ValueConst(*R_,curDim+j,curDim+j) - mu;
          } else {
            vscale = -sigma / (DMT::Value(*R_,curDim+j, curDim+j) + mu);
          }
          beta[curDim+j] = two * vscale * vscale / (sigma + vscale*vscale);
          DMT::Value(*R_,curDim+j, curDim+j) = maxelem*mu;
          for (int i = 0; i < blockSize_; ++i) {
            DMT::Value(*R_,curDim+j+1+i,curDim+j) /= vscale;
          }
        }
 
        // Apply new Householder reflector to the right-hand side.
        for (int i = 0; i < blockSize_; ++i) {
          sigma = blas.DOT (blockSize_, &DMT::Value(*R_,curDim+j+1,curDim+j),
                            1, &DMT::Value(*z_,curDim+j+1,i), 1);
          sigma += DMT::ValueConst(*z_,curDim+j,i);
          sigma *= beta[curDim+j];
          blas.AXPY (blockSize_, ScalarType(-sigma), &DMT::Value(*R_,curDim+j+1,curDim+j),
                     1, &DMT::Value(*z_,curDim+j+1,i), 1);
          DMT::Value(*z_,curDim+j,i) -= sigma;
        }
      }
    } // end if (blockSize_ == 1)
 
    DMT::SyncHostToDevice( *H_ );
    DMT::SyncHostToDevice( *z_ );
 
    // If the least-squares problem is updated wrt "dim" then update curDim_.
    if (dim >= curDim_ && dim < getMaxSubspaceDim ()) {
      curDim_ = dim + blockSize_;
    }
  } // end updateLSQR()
 
} // namespace Belos
 
#endif /* BELOS_BLOCK_FGMRES_ITER_HPP */