BelosLinearProblem.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_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 setLeftPrec(const ::Teuchos::RCP<const OP> &LP) {  LP_ = LP; }

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

    void setCurrLS ();

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

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

    ::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 ); }

    

    
    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 { return(R0_); }
    
    ::Teuchos::RCP<const MV> getInitPrecResVec() const { return(PR0_); }
    
    ::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); }
 
    

    

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

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

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


    void computeCurrPrecResVec( MV* R, const MV* X = 0, const MV* B = 0 ) const;
    
    
  private:
    
    ::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 OP> LP_;  
    
    ::Teuchos::RCP<const OP> RP_;
    
    mutable ::Teuchos::RCP<::Teuchos::Time> timerOp_, timerPrec_;

    int blocksize_;

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

    int lsNum_;



    bool Left_Scale_;

    bool Right_Scale_;

    bool isSet_;

    bool isHermitian_;

    bool solutionUpdated_;    

   
    std::string label_;
 
    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),
    Left_Scale_(false),
    Right_Scale_(false),
    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),
    Left_Scale_(false),
    Right_Scale_(false),
    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_),
    LP_(Problem.LP_),
    RP_(Problem.RP_),
    timerOp_(Problem.timerOp_),
    timerPrec_(Problem.timerPrec_),
    blocksize_(Problem.blocksize_),
    num2Solve_(Problem.num2Solve_),
    rhsIndex_(Problem.rhsIndex_),
    lsNum_(Problem.lsNum_),
    Left_Scale_(Problem.Left_Scale_),
    Right_Scale_(Problem.Right_Scale_),
    isSet_(Problem.isSet_),
    isHermitian_(Problem.isHermitian_),
    solutionUpdated_(Problem.solutionUpdated_),
    label_(Problem.label_)
  {
  }
  
  template <class ScalarType, class MV, class OP>
  LinearProblem<ScalarType,MV,OP>::~LinearProblem(void)
  {}
  
  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));
    {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
      ::Teuchos::TimeMonitor PrecTimer (*timerPrec_);
#endif
      OPT::Apply( *RP_, *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));
    {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
      ::Teuchos::TimeMonitor PrecTimer(*timerPrec_);
#endif
      OPT::Apply( *RP_, *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
    }

    // 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 (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 || MVT::GetNumberVecs( *PR0_ )!=MVT::GetNumberVecs( *B_ )) {
        PR0_ = MVT::Clone( *B_, MVT::GetNumberVecs( *B_ ) );
      }
      {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
        ::Teuchos::TimeMonitor PrecTimer(*timerPrec_);
#endif
        OPT::Apply( *LP_, *R0_, *PR0_ );
      }
    } 
    else {
      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<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;

    // We only need a temporary vector for intermediate results if
    // there is a left or right preconditioner.  We really should just
    // keep this temporary vector around instead of allocating it each
    // time.
    RCP<MV> ytemp = (leftPrec || rightPrec) ? MVT::Clone (y, MVT::GetNumberVecs (y)) : null;

    //
    // No preconditioning.
    // 
    if (!leftPrec && !rightPrec){ 
#ifdef BELOS_TEUCHOS_TIME_MONITOR
      ::Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
      OPT::Apply( *A_, x, y );
    }
    //
    // Preconditioning is being done on both sides
    //
    else if( leftPrec && rightPrec ) 
      {
        {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
          ::Teuchos::TimeMonitor PrecTimer(*timerPrec_);
#endif
    OPT::Apply( *RP_, x, y );   
        }
        {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
          ::Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
    OPT::Apply( *A_, y, *ytemp );
        }
        {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
          ::Teuchos::TimeMonitor PrecTimer(*timerPrec_);
#endif
    OPT::Apply( *LP_, *ytemp, y );
        }
      }
    //
    // Preconditioning is only being done on the left side
    //
    else if( leftPrec ) 
      {
        {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
          ::Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
    OPT::Apply( *A_, x, *ytemp );
        }
        {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
          ::Teuchos::TimeMonitor PrecTimer(*timerPrec_);
#endif
    OPT::Apply( *LP_, *ytemp, y );
        }
      }
    //
    // Preconditioning is only being done on the right side
    //
    else 
      {
        {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
          ::Teuchos::TimeMonitor PrecTimer(*timerPrec_);
#endif
    OPT::Apply( *RP_, x, *ytemp );
        }
        {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
          ::Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
          OPT::Apply( *A_, *ytemp, 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
      return ( 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
      return ( 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 ) );
              {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
                ::Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
          OPT::Apply( *A_, *X, *R_temp );
              }
        MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
              {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
                ::Teuchos::TimeMonitor PrecTimer(*timerPrec_);
#endif
          OPT::Apply( *LP_, *R_temp, *R );
              }
      }
    else 
      {
              {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
                ::Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
          OPT::Apply( *A_, *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 ) );
            {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
              ::Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
        OPT::Apply( *A_, *localX, *R_temp );
            }
      MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
            {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
              ::Teuchos::TimeMonitor PrecTimer(*timerPrec_);
#endif
        OPT::Apply( *LP_, *R_temp, *R );
            }
    }
  else 
    {
            {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
              ::Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
          OPT::Apply( *A_, *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
  {
          {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
            ::Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
      OPT::Apply( *A_, *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_;
    
        {
#ifdef BELOS_TEUCHOS_TIME_MONITOR
          ::Teuchos::TimeMonitor OpTimer(*timerOp_);
#endif
    OPT::Apply( *A_, *localX, *R );
        }
  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
      }    
    }
  }
  
} // end Belos namespace

#endif /* BELOS_LINEAR_PROBLEM_HPP */