BelosLinearProblem.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_LINEAR_PROBLEM_HPP
#define BELOS_LINEAR_PROBLEM_HPP

#include "BelosMultiVecTraits.hpp"
#include "BelosOperatorTraits.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_TimeMonitor.hpp"

namespace Belos {


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

  template <class ScalarType, class MV, class OP>
  class LinearProblem {
  public:
    


    LinearProblem (void);
    
    LinearProblem (const Teuchos::RCP<const OP> &A, 
       const Teuchos::RCP<MV> &X, 
       const Teuchos::RCP<const MV> &B);
    
    LinearProblem (const LinearProblem<ScalarType,MV,OP>& Problem);
    
    virtual ~LinearProblem (void) {}

    

    
    void setOperator (const Teuchos::RCP<const OP> &A) { 
      A_ = A; 
      isSet_=false; 
    }
    
    void setLHS (const Teuchos::RCP<MV> &X) { 
      X_ = X; 
      isSet_=false; 
    }
    
    void setRHS (const Teuchos::RCP<const MV> &B) { 
      B_ = B; 
      isSet_=false; 
    }

    void setInitResVec(const Teuchos::RCP<const MV> &R0) {
      R0_user_ = R0;
      isSet_=false;
    }

    void setInitPrecResVec(const Teuchos::RCP<const MV> &PR0) {
      PR0_user_ = PR0;
      isSet_=false;
    }
    
    void setLeftPrec(const Teuchos::RCP<const OP> &LP) {  LP_ = LP; }

    void setRightPrec(const Teuchos::RCP<const OP> &RP) { RP_ = RP; }

    virtual void setCurrLS ();

    virtual void setLSIndex (const std::vector<int>& index); 
    
    void setHermitian( bool isSym = true ) { isHermitian_ = isSym; }
   
    void setLabel (const std::string& label);

    virtual 
    Teuchos::RCP<MV> 
    updateSolution (const Teuchos::RCP<MV>& update = Teuchos::null,
        bool updateLP = false,
        ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one());    

    virtual
    Teuchos::RCP<MV> 
    updateSolution( const Teuchos::RCP<MV>& update = Teuchos::null,
                    ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() ) const
    { return const_cast<LinearProblem<ScalarType,MV,OP> *>(this)->updateSolution( update, false, scale ); }

    

    
    virtual
    bool 
    setProblem (const Teuchos::RCP<MV> &newX = Teuchos::null, 
    const Teuchos::RCP<const MV> &newB = Teuchos::null);

    

    
    Teuchos::RCP<const OP> getOperator() const { return(A_); }
    
    Teuchos::RCP<MV> getLHS() const { return(X_); }
    
    Teuchos::RCP<const MV> getRHS() const { return(B_); }
    
    Teuchos::RCP<const MV> getInitResVec() const;
    
    Teuchos::RCP<const MV> getInitPrecResVec() const;
    
    Teuchos::RCP<MV> getCurrLHSVec();
    
    Teuchos::RCP<const MV> getCurrRHSVec();
    
    Teuchos::RCP<const OP> getLeftPrec() const { return(LP_); };
    
    Teuchos::RCP<const OP> getRightPrec() const { return(RP_); };
    
    const std::vector<int>& getLSIndex() const { return(rhsIndex_); }

    int getLSNumber() const { return(lsNum_); }

    Teuchos::Array<Teuchos::RCP<Teuchos::Time> > getTimers() const {
      return Teuchos::tuple(timerOp_,timerPrec_);
    }


    

    
    bool isSolutionUpdated() const { return(solutionUpdated_); }

    bool isProblemSet() const { return(isSet_); }

    bool isHermitian() const { return(isHermitian_); }
    
    bool isLeftPrec() const { return(LP_!=Teuchos::null); }

    bool isRightPrec() const { return(RP_!=Teuchos::null); }
 
    

    

    virtual void apply( const MV& x, MV& y ) const;
    
    virtual void applyOp( const MV& x, MV& y ) const;
    
    virtual void applyLeftPrec( const MV& x, MV& y ) const;

    virtual void applyRightPrec( const MV& x, MV& y ) const;
    

    virtual void computeCurrResVec( MV* R , const MV* X = 0, const MV* B = 0 ) const;


    virtual void computeCurrPrecResVec( MV* R, const MV* X = 0, const MV* B = 0 ) const;
    
  
  protected:

    Teuchos::RCP<const OP> A_;
    
    Teuchos::RCP<MV> X_;
    
    Teuchos::RCP<MV> curX_;
    
    Teuchos::RCP<const MV> B_;
    
    Teuchos::RCP<const MV> curB_;
    
    Teuchos::RCP<MV> R0_;
   
    Teuchos::RCP<MV> PR0_;

    Teuchos::RCP<const MV> R0_user_;

    Teuchos::RCP<const MV> PR0_user_;
 
    Teuchos::RCP<const OP> LP_;  
    
    Teuchos::RCP<const OP> RP_;
    
    mutable Teuchos::RCP<MV> Y_temp_;

    mutable Teuchos::RCP<Teuchos::Time> timerOp_, timerPrec_;

    int blocksize_;

    int num2Solve_;
    
    std::vector<int> rhsIndex_;    

    int lsNum_;



    bool isSet_;

    bool isHermitian_;

    bool solutionUpdated_;    

   
    std::string label_;

  private:
 
    typedef MultiVecTraits<ScalarType,MV>  MVT;
    typedef OperatorTraits<ScalarType,MV,OP>  OPT;
  };
  
  //--------------------------------------------
  //  Constructor Implementations
  //--------------------------------------------
  
  template <class ScalarType, class MV, class OP>
  LinearProblem<ScalarType,MV,OP>::LinearProblem(void) : 
    blocksize_(0),
    num2Solve_(0),
    rhsIndex_(0),
    lsNum_(0),
    isSet_(false),
    isHermitian_(false),
    solutionUpdated_(false),
    label_("Belos")
  {
  }
  
  template <class ScalarType, class MV, class OP>
  LinearProblem<ScalarType,MV,OP>::LinearProblem(const Teuchos::RCP<const OP> &A, 
             const Teuchos::RCP<MV> &X, 
             const Teuchos::RCP<const MV> &B
             ) :
    A_(A),
    X_(X),
    B_(B),
    blocksize_(0),
    num2Solve_(0),
    rhsIndex_(0),
    lsNum_(0),
    isSet_(false),
    isHermitian_(false),
    solutionUpdated_(false),
    label_("Belos")
  {
  }
  
  template <class ScalarType, class MV, class OP>
  LinearProblem<ScalarType,MV,OP>::LinearProblem(const LinearProblem<ScalarType,MV,OP>& Problem) :
    A_(Problem.A_),
    X_(Problem.X_),
    curX_(Problem.curX_),
    B_(Problem.B_),
    curB_(Problem.curB_),
    R0_(Problem.R0_),
    PR0_(Problem.PR0_),
    R0_user_(Problem.R0_user_),
    PR0_user_(Problem.PR0_user_),
    LP_(Problem.LP_),
    RP_(Problem.RP_),
    timerOp_(Problem.timerOp_),
    timerPrec_(Problem.timerPrec_),
    blocksize_(Problem.blocksize_),
    num2Solve_(Problem.num2Solve_),
    rhsIndex_(Problem.rhsIndex_),
    lsNum_(Problem.lsNum_),
    isSet_(Problem.isSet_),
    isHermitian_(Problem.isHermitian_),
    solutionUpdated_(Problem.solutionUpdated_),
    label_(Problem.label_)
  {
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::setLSIndex(const std::vector<int>& index)
  {
    // Set new linear systems using the indices in index.
    rhsIndex_ = index;
    
    // Compute the new block linear system.
    // ( first clean up old linear system )
    curB_ = Teuchos::null;
    curX_ = Teuchos::null;
   
    // Create indices for the new linear system.
    int validIdx = 0, ivalidIdx = 0;
    blocksize_ = rhsIndex_.size();
    std::vector<int> vldIndex( blocksize_ );
    std::vector<int> newIndex( blocksize_ );
    std::vector<int> iIndex( blocksize_ );
    for (int i=0; i<blocksize_; ++i) {
      if (rhsIndex_[i] > -1) {
        vldIndex[validIdx] = rhsIndex_[i];
        newIndex[validIdx] = i;
        validIdx++;
      }
      else {
        iIndex[ivalidIdx] = i;
        ivalidIdx++;
      }
    }
    vldIndex.resize(validIdx);
    newIndex.resize(validIdx);   
    iIndex.resize(ivalidIdx);
    num2Solve_ = validIdx;

    // Create the new linear system using index
    if (num2Solve_ != blocksize_) {
      newIndex.resize(num2Solve_);
      vldIndex.resize(num2Solve_);
      //
      // First create multivectors of blocksize.
      // Fill the RHS with random vectors LHS with zero vectors.
      curX_ = MVT::Clone( *X_, blocksize_ );
      MVT::MvInit(*curX_);
      Teuchos::RCP<MV> tmpCurB = MVT::Clone( *B_, blocksize_ );
      MVT::MvRandom(*tmpCurB);
      //
      // Now put in the part of B into curB 
      Teuchos::RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
      MVT::SetBlock( *tptr, newIndex, *tmpCurB );
      curB_ = tmpCurB;
      //
      // Now put in the part of X into curX
      tptr = MVT::CloneView( *X_, vldIndex );
      MVT::SetBlock( *tptr, newIndex, *curX_ );
      //
      solutionUpdated_ = false;
    }
    else {
      curX_ = MVT::CloneViewNonConst( *X_, rhsIndex_ );
      curB_ = MVT::CloneView( *B_, rhsIndex_ );
    }
    //
    // Increment the number of linear systems that have been loaded into this object.
    //
    lsNum_++;
  }


  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::setCurrLS() 
  { 
    //
    // We only need to copy the solutions back if the linear systems of
    // interest are less than the block size.
    //
    if (num2Solve_ < blocksize_) {
      //
      // Get a view of the current solutions and correction std::vector.
      //
      int validIdx = 0;
      std::vector<int> newIndex( num2Solve_ );
      std::vector<int> vldIndex( num2Solve_ );
      for (int i=0; i<blocksize_; ++i) {
        if ( rhsIndex_[i] > -1 ) { 
          vldIndex[validIdx] = rhsIndex_[i];
    newIndex[validIdx] = i;
          validIdx++;
        } 
      }
      Teuchos::RCP<MV> tptr = MVT::CloneViewNonConst( *curX_, newIndex );
      MVT::SetBlock( *tptr, vldIndex, *X_ );
    }
    //
    // Clear the current vectors of this linear system so that any future calls
    // to get the vectors for this system return null pointers.
    //
    curX_ = Teuchos::null;
    curB_ = Teuchos::null;
    rhsIndex_.resize(0);
  }
  

  template <class ScalarType, class MV, class OP>
  Teuchos::RCP<MV> 
  LinearProblem<ScalarType,MV,OP>::
  updateSolution (const Teuchos::RCP<MV>& update, 
      bool updateLP,
      ScalarType scale)
  { 
    using Teuchos::RCP;
    using Teuchos::null;

    RCP<MV> newSoln;
    if (update.is_null())
      { // The caller didn't supply an update vector, so we assume
  // that the current solution curX_ is unchanged, and return a
  // pointer to it.
  newSoln = curX_;
      }
    else // the update vector is NOT null.
      { 
  if (updateLP) // The caller wants us to update the linear problem.
    { 
      if (RP_.is_null()) 
        { // There is no right preconditioner.
    // curX_ := curX_ + scale * update.
    MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ ); 
        }
      else 
        { // There is indeed a right preconditioner, so apply it
    // before computing the new solution.
    RCP<MV> rightPrecUpdate = 
      MVT::Clone (*update, MVT::GetNumberVecs (*update));
                applyRightPrec( *update, *rightPrecUpdate ); 
    // curX_ := curX_ + scale * rightPrecUpdate.
    MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *curX_ ); 
        } 
      solutionUpdated_ = true; 
      newSoln = curX_;
    }
  else 
    { // Rather than updating our stored current solution curX_,
      // we make a copy and compute the new solution in the
      // copy, without modifying curX_.
      newSoln = MVT::Clone (*update, MVT::GetNumberVecs (*update));
      if (RP_.is_null())
        { // There is no right preconditioner.
    // newSoln := curX_ + scale * update.
    MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln ); 
        }
      else 
        { // There is indeed a right preconditioner, so apply it
    // before computing the new solution.
    RCP<MV> rightPrecUpdate = 
      MVT::Clone (*update, MVT::GetNumberVecs (*update));
                applyRightPrec( *update, *rightPrecUpdate ); 
    // newSoln := curX_ + scale * rightPrecUpdate.
    MVT::MvAddMv( 1.0, *curX_, scale, *rightPrecUpdate, *newSoln ); 
        } 
    }
      }
    return newSoln;
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::setLabel(const std::string& label)
  {
    if (label != label_) {
      label_ = label;
      // Create new timers if they have already been created.
      if (timerOp_ != Teuchos::null) {
        std::string opLabel = label_ + ": Operation Op*x";
#ifdef BELOS_TEUCHOS_TIME_MONITOR
        timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
#endif
      }
      if (timerPrec_ != Teuchos::null) {
        std::string precLabel = label_ + ": Operation Prec*x";
#ifdef BELOS_TEUCHOS_TIME_MONITOR
        timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
#endif
      }
    }
  }

  template <class ScalarType, class MV, class OP>
  bool 
  LinearProblem<ScalarType,MV,OP>::
  setProblem (const Teuchos::RCP<MV> &newX, 
        const Teuchos::RCP<const MV> &newB)
  {
    // Create timers if the haven't been created yet.
    if (timerOp_ == Teuchos::null) {
      std::string opLabel = label_ + ": Operation Op*x";
#ifdef BELOS_TEUCHOS_TIME_MONITOR
      timerOp_ = Teuchos::TimeMonitor::getNewCounter( opLabel );
#endif
    }
    if (timerPrec_ == Teuchos::null) {
      std::string precLabel = label_ + ": Operation Prec*x";
#ifdef BELOS_TEUCHOS_TIME_MONITOR
      timerPrec_ = Teuchos::TimeMonitor::getNewCounter( precLabel );
#endif
    }

    Y_temp_ = Teuchos::null;

    // Set the linear system using the arguments newX and newB
    if (newX != Teuchos::null)
      X_ = newX;
    if (newB != Teuchos::null)
      B_ = newB;

    // Invalidate the current linear system indices and multivectors.
    rhsIndex_.resize(0);
    curX_ = Teuchos::null;
    curB_ = Teuchos::null;

    // If we didn't set a matrix A, a left-hand side X, or a
    // right-hand side B, then we didn't set the problem.
    if (A_ == Teuchos::null || X_ == Teuchos::null || B_ == Teuchos::null) {
      isSet_ = false;
      return isSet_;
    }

    // Reset whether the solution has been updated.  (We're just
    // setting the problem now, so of course we haven't updated the
    // solution yet.)
    solutionUpdated_ = false;
    
    // Compute the initial residuals.
    if(Teuchos::is_null(R0_user_)) {
      if (R0_==Teuchos::null || MVT::GetNumberVecs( *R0_ )!=MVT::GetNumberVecs( *B_ )) {
        R0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
      }
      computeCurrResVec( &*R0_, &*X_, &*B_ );

      if (LP_!=Teuchos::null) {
        if (PR0_==Teuchos::null || (PR0_==R0_) || (MVT::GetNumberVecs(*PR0_)!=MVT::GetNumberVecs(*B_))) {
          PR0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
        }
        applyLeftPrec( *R0_, *PR0_ );
      } 
      else {
        PR0_ = R0_;
      }
    }
    else { // User specified the residuals
      // If the user did not specify the right sized residual, create one and set R0_user_ to null.
      if (MVT::GetNumberVecs( *R0_user_ )!=MVT::GetNumberVecs( *B_ )) {
        Teuchos::RCP<MV> helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
        computeCurrResVec( &*helper, &*X_, &*B_ );
        R0_user_ = Teuchos::null;
        R0_ = helper;
      }

      if (LP_!=Teuchos::null) {
        // If the user provided preconditioned residual is the wrong size or pointing at
        // the wrong object, create one and set the PR0_user_ to null.
        if (PR0_user_==Teuchos::null || (PR0_user_==R0_) || (PR0_user_==R0_user_) 
          || (MVT::GetNumberVecs(*PR0_user_)!=MVT::GetNumberVecs(*B_))) {
          Teuchos::RCP<MV> helper = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
   
          // Get the initial residual from getInitResVec because R0_user_ may be null from above.
          applyLeftPrec( *getInitResVec(), *helper );
          PR0_user_ = Teuchos::null;
          PR0_ = helper;
        }
      }
      else {
        // The preconditioned initial residual vector is the residual vector.
        // NOTE:  R0_user_ could be null if incompatible.
        if (R0_user_!=Teuchos::null)  
        {
          PR0_user_ = R0_user_;
        }
        else
        {
          PR0_user_ = Teuchos::null;
          PR0_ = R0_;
        }
      }
    }

    // The problem has been set and is ready for use.
    isSet_ = true;
    
    // Return isSet.
    return isSet_;
  }

  template <class ScalarType, class MV, class OP>
  Teuchos::RCP<const MV> LinearProblem<ScalarType,MV,OP>::getInitResVec() const 
  {
    if(Teuchos::nonnull(R0_user_)) {
      return R0_user_;
    }
    return(R0_); 
  }

  template <class ScalarType, class MV, class OP>
  Teuchos::RCP<const MV> LinearProblem<ScalarType,MV,OP>::getInitPrecResVec() const 
  { 
    if(Teuchos::nonnull(PR0_user_)) {
      return PR0_user_;
    }
    return(PR0_); 
  }
  
  template <class ScalarType, class MV, class OP>
  Teuchos::RCP<MV> LinearProblem<ScalarType,MV,OP>::getCurrLHSVec()
  {
    if (isSet_) {
      return curX_;
    }
    else {
      return Teuchos::null;
    }
  }
  
  template <class ScalarType, class MV, class OP>
  Teuchos::RCP<const MV> LinearProblem<ScalarType,MV,OP>::getCurrRHSVec()
  {
    if (isSet_) {
      return curB_;
    }
    else {
      return Teuchos::null;
    }
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::apply( const MV& x, MV& y ) const
  {
    using Teuchos::null;
    using Teuchos::RCP;

    const bool leftPrec = LP_ != null;
    const bool rightPrec = RP_ != null;

    if(leftPrec || rightPrec) {
      const auto num_vecs_y = MVT::GetNumberVecs(y);
      if(!Y_temp_ || MVT::GetNumberVecs(*Y_temp_) != num_vecs_y) {
        Y_temp_ = MVT::Clone (y, num_vecs_y);
      }
    }

    //
    // No preconditioning.
    // 
    if (!leftPrec && !rightPrec) { 
      applyOp( x, y );
    }
    //
    // Preconditioning is being done on both sides
    //
    else if( leftPrec && rightPrec ) 
    {
      applyRightPrec( x, y );
      applyOp( y, *Y_temp_ );
      applyLeftPrec( *Y_temp_, y );
    }
    //
    // Preconditioning is only being done on the left side
    //
    else if( leftPrec ) 
    {
      applyOp( x, *Y_temp_ );
      applyLeftPrec( *Y_temp_, y );
    }
    //
    // Preconditioning is only being done on the right side
    //
    else 
    {
      applyRightPrec( x, *Y_temp_ );
      applyOp( *Y_temp_, y );
    }  
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::applyOp( const MV& x, MV& y ) const {
    if (A_.get()) {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
      Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
      OPT::Apply( *A_,x, y);   
    }
    else {
      MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
        Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
    }
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::applyLeftPrec( const MV& x, MV& y ) const {
    if (LP_!=Teuchos::null) {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
      Teuchos::TimeMonitor PrecTimer(*timerPrec_);
#endif
      OPT::Apply( *LP_,x, y);
    }
    else {
      MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
        Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
    }
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::applyRightPrec( const MV& x, MV& y ) const {
    if (RP_!=Teuchos::null) {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
      Teuchos::TimeMonitor PrecTimer(*timerPrec_);
#endif
      OPT::Apply( *RP_,x, y);
    }
    else {
      MVT::MvAddMv( Teuchos::ScalarTraits<ScalarType>::one(), x, 
        Teuchos::ScalarTraits<ScalarType>::zero(), x, y );
    }
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::computeCurrPrecResVec( MV* R, const MV* X, const MV* B ) const {

    if (R) {
      if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
      {
  if (LP_!=Teuchos::null)
  {
    Teuchos::RCP<MV> R_temp = MVT::Clone( *B, MVT::GetNumberVecs( *B ) );
          applyOp( *X, *R_temp );
    MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
          applyLeftPrec( *R_temp, *R );
  }
  else 
  {
          applyOp( *X, *R );
    MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
  }
      }
      else { 
  // The solution and right-hand side may not be specified, check and use which ones exist.
  Teuchos::RCP<const MV> localB, localX;
  if (B)
    localB = Teuchos::rcp( B, false );
  else
    localB = curB_;
  
  if (X)
    localX = Teuchos::rcp( X, false );
  else
    localX = curX_;
  
  if (LP_!=Teuchos::null)
  {
    Teuchos::RCP<MV> R_temp = MVT::Clone( *localB, MVT::GetNumberVecs( *localB ) );
          applyOp( *localX, *R_temp );
    MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
          applyLeftPrec( *R_temp, *R );
  }
  else 
  {
          applyOp( *localX, *R );
    MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
  }
      }    
    } 
  }
  
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::computeCurrResVec( MV* R, const MV* X, const MV* B ) const {

    if (R) {
      if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
      {
        applyOp( *X, *R );
  MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
      }
      else { 
  // The solution and right-hand side may not be specified, check and use which ones exist.
  Teuchos::RCP<const MV> localB, localX;
  if (B)
    localB = Teuchos::rcp( B, false );
  else
    localB = curB_;
  
  if (X)
    localX = Teuchos::rcp( X, false );
  else
    localX = curX_;

        applyOp( *localX, *R );   
  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
      }    
    }
  }
  
} // end Belos namespace

#endif /* BELOS_LINEAR_PROBLEM_HPP */