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

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

using Teuchos::RCP;
using Teuchos::rcp;
using Teuchos::null;
using Teuchos::rcp_const_cast;
using Teuchos::ParameterList;

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 RCP<const OP> &A, 
      const RCP<MV> &X, 
      const RCP<const MV> &B
      );
    

    LinearProblem(const LinearProblem<ScalarType,MV,OP>& Problem);
    

    virtual ~LinearProblem(void);
    

    

    void setOperator(const RCP<const OP> &A) { A_ = A; isSet_=false; }
    

    void setLHS(const RCP<MV> &X) { X_ = X; isSet_=false; }
    

    void setRHS(const RCP<const MV> &B) { B_ = B; isSet_=false; }
    

    void setLeftPrec(const RCP<const OP> &LP) {  LP_ = LP; }
    

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

    void setCurrLS();
    

    void setLSIndex(std::vector<int>& index); 
    

    void setHermitian(){ isHermitian_ = true; }
   

    void setLabel(const std::string& label) { label_ = label; }


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

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

    

    

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

    

    
    RCP<const OP> getOperator() const { return(A_); }
    
    RCP<MV> getLHS() const { return(X_); }
    
    RCP<const MV> getRHS() const { return(B_); }
    

    RCP<const MV> getInitResVec() const { return(R0_); }
    

    RCP<const MV> getInitPrecResVec() const { return(PR0_); }
    

    RCP<MV> getCurrLHSVec();
    

    RCP<MV> getCurrRHSVec();
    
    RCP<const OP> getLeftPrec() const { return(LP_); };
    
    RCP<const OP> getRightPrec() const { return(RP_); };
    

    const std::vector<int> getLSIndex() const { return(rhsIndex_); }

    /* This can be used by status test classes to determine if the solver manager has advanced 
       and is solving another linear system.
    */
    int getLSNumber() const { return(lsNum_); }

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


    

    

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

    bool isProblemSet() const { return(isSet_); }
    
    bool isHermitian() const { return(isHermitian_); }
    
    bool isLeftPrec() const { return(LP_!=null); }

    bool isRightPrec() const { return(RP_!=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:
    
    RCP<const OP> A_;
    
    RCP<MV> X_;
    
    RCP<MV> curX_;
    
    RCP<const MV> B_;
    
    RCP<MV> curB_;
    
    RCP<MV> R0_;
   
    RCP<MV> PR0_;
 
    RCP<const OP> LP_;  
    
    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 RCP<const OP> &A, 
             const RCP<MV> &X, 
             const 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(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_ = null;
    curX_ = 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_);
      curB_ = MVT::Clone( *B_, blocksize_ );
      MVT::MvRandom(*curB_);
      //
      // Now put in the part of B into curB 
      RCP<const MV> tptr = MVT::CloneView( *B_, vldIndex );
      MVT::SetBlock( *tptr, newIndex, *curB_ );
      //
      // Now put in the part of X into curX
      tptr = MVT::CloneView( *X_, vldIndex );
      MVT::SetBlock( *tptr, newIndex, *curX_ );
      //
      solutionUpdated_ = false;
    }
    else {
      curX_ = MVT::CloneView( *X_, rhsIndex_ );
      curB_ = rcp_const_cast<MV>(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++;
        } 
      }
      RCP<MV> tptr = MVT::CloneView( *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_ = null;
    curB_ = null;
    rhsIndex_.resize(0);
  }
  

  template <class ScalarType, class MV, class OP>
  RCP<MV> LinearProblem<ScalarType,MV,OP>::updateSolution( const RCP<MV>& update, 
                 bool updateLP,
                 ScalarType scale )
  { 
    RCP<MV> newSoln;
    if (update != null) {
      if (updateLP == true) {
  if (RP_!=null) {
    //
    // Apply the right preconditioner before computing the current solution.
    RCP<MV> TrueUpdate = MVT::Clone( *update, MVT::GetNumberVecs( *update ) );
    {
      Teuchos::TimeMonitor PrecTimer(*timerPrec_);
      OPT::Apply( *RP_, *update, *TrueUpdate ); 
    }
    MVT::MvAddMv( 1.0, *curX_, scale, *TrueUpdate, *curX_ ); 
  } 
  else {
    MVT::MvAddMv( 1.0, *curX_, scale, *update, *curX_ ); 
  }
  solutionUpdated_ = true; 
  newSoln = curX_;
      }
      else {
  newSoln = MVT::Clone( *update, MVT::GetNumberVecs( *update ) );
  if (RP_!=null) {
    //
    // Apply the right preconditioner before computing the current solution.
    RCP<MV> trueUpdate = MVT::Clone( *update, MVT::GetNumberVecs( *update ) );
    {
      Teuchos::TimeMonitor PrecTimer(*timerPrec_);
      OPT::Apply( *RP_, *update, *trueUpdate ); 
    }
    MVT::MvAddMv( 1.0, *curX_, scale, *trueUpdate, *newSoln ); 
  } 
  else {
    MVT::MvAddMv( 1.0, *curX_, scale, *update, *newSoln ); 
  }
      }
    }
    else {
      newSoln = curX_;
    }
    return newSoln;
  }
  

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

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

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

    // Check the validity of the linear problem object.
    // If no operator A exists, then throw an std::exception.
    if (A_ == null || X_ == null || B_ == null) {
      isSet_ = false;
      return isSet_;
    }

    // Initialize the state booleans
    solutionUpdated_ = false;
    
    // Compute the initial residuals.
    if (R0_==null || MVT::GetNumberVecs( *R0_ )!=MVT::GetNumberVecs( *X_ )) {
      R0_ = MVT::Clone( *X_, MVT::GetNumberVecs( *X_ ) );
    }
    computeCurrResVec( &*R0_, &*X_, &*B_ );

    if (LP_!=null) {
      if (PR0_==null || MVT::GetNumberVecs( *PR0_ )!=MVT::GetNumberVecs( *X_ )) {
        PR0_ = MVT::Clone( *X_, MVT::GetNumberVecs( *X_ ) );
      }
      {
        Teuchos::TimeMonitor PrecTimer(*timerPrec_);
        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>
  RCP<MV> LinearProblem<ScalarType,MV,OP>::getCurrLHSVec()
  {
    if (isSet_) {
      return curX_;
    }
    else {
      return Teuchos::null;
    }
  }
  
  template <class ScalarType, class MV, class OP>
  RCP<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
  {
    RCP<MV> ytemp = MVT::Clone( y, MVT::GetNumberVecs( y ) );
    bool leftPrec = LP_!=null;
    bool rightPrec = RP_!=null;
    //
    // No preconditioning.
    // 
    if (!leftPrec && !rightPrec){ 
      Teuchos::TimeMonitor OpTimer(*timerOp_);
      OPT::Apply( *A_, x, y );
    }
    //
    // Preconditioning is being done on both sides
    //
    else if( leftPrec && rightPrec ) 
      {
        {
          Teuchos::TimeMonitor PrecTimer(*timerPrec_);
    OPT::Apply( *RP_, x, y );   
        }
        {
          Teuchos::TimeMonitor OpTimer(*timerOp_);
    OPT::Apply( *A_, y, *ytemp );
        }
        {
          Teuchos::TimeMonitor PrecTimer(*timerPrec_);
    OPT::Apply( *LP_, *ytemp, y );
        }
      }
    //
    // Preconditioning is only being done on the left side
    //
    else if( leftPrec ) 
      {
        {
          Teuchos::TimeMonitor OpTimer(*timerOp_);
    OPT::Apply( *A_, x, *ytemp );
        }
        {
          Teuchos::TimeMonitor PrecTimer(*timerPrec_);
    OPT::Apply( *LP_, *ytemp, y );
        }
      }
    //
    // Preconditioning is only being done on the right side
    //
    else 
      {
        {
          Teuchos::TimeMonitor PrecTimer(*timerPrec_);
    OPT::Apply( *RP_, x, *ytemp );
        }
        {
          Teuchos::TimeMonitor OpTimer(*timerOp_);
          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()) {
      Teuchos::TimeMonitor OpTimer(*timerOp_);
      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_!=null) {
      Teuchos::TimeMonitor PrecTimer(*timerPrec_);
      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_!=null) {
      Teuchos::TimeMonitor PrecTimer(*timerPrec_);
      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_!=null)
      {
        RCP<MV> R_temp = MVT::Clone( *X, MVT::GetNumberVecs( *X ) );
              {
                Teuchos::TimeMonitor OpTimer(*timerOp_);
          OPT::Apply( *A_, *X, *R_temp );
              }
        MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
              {
                Teuchos::TimeMonitor PrecTimer(*timerPrec_);
          OPT::Apply( *LP_, *R_temp, *R );
              }
      }
    else 
      {
              {
                Teuchos::TimeMonitor OpTimer(*timerOp_);
          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.
  RCP<const MV> localB, localX;
  if (B)
    localB = rcp( B, false );
  else
    localB = curB_;
  
  if (X)
    localX = rcp( X, false );
  else
    localX = curX_;
  
  if (LP_!=null)
    {
      RCP<MV> R_temp = MVT::Clone( *localX, MVT::GetNumberVecs( *localX ) );
            {
              Teuchos::TimeMonitor OpTimer(*timerOp_);
        OPT::Apply( *A_, *localX, *R_temp );
            }
      MVT::MvAddMv( -1.0, *R_temp, 1.0, *localB, *R_temp );
            {
              Teuchos::TimeMonitor PrecTimer(*timerPrec_);
        OPT::Apply( *LP_, *R_temp, *R );
            }
    }
  else 
    {
            {
              Teuchos::TimeMonitor OpTimer(*timerOp_);
          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
  {
          {
            Teuchos::TimeMonitor OpTimer(*timerOp_);
      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.
  RCP<const MV> localB, localX;
  if (B)
    localB = rcp( B, false );
  else
    localB = curB_;
  
  if (X)
    localX = rcp( X, false );
  else
    localX = curX_;
    
        {
          Teuchos::TimeMonitor OpTimer(*timerOp_);
    OPT::Apply( *A_, *localX, *R );
        }
  MVT::MvAddMv( -1.0, *R, 1.0, *localB, *R );
      }    
    }
  }
  
} // end Belos namespace

#endif /* BELOS_LINEAR_PROBLEM_HPP */

---