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"

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

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



    LinearProblem(void);
    

    LinearProblem(const RefCountPtr<const OP> &A, 
      const RefCountPtr<MV> &X, 
      const RefCountPtr<const MV> &B
      );
    

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

    virtual ~LinearProblem(void);
    

    

    void SetOperator(const RefCountPtr<const OP> &A) { A_ = A; };
    

    void SetLHS(const RefCountPtr<MV> &X);
    

    void SetRHS(const RefCountPtr<const MV> &B) { B_ = B; };
    

    void SetLeftPrec(const RefCountPtr<const OP> &LP) {  LP_ = LP; Left_Prec_ = true; };
    

    void SetRightPrec(const RefCountPtr<const OP> &RP) { RP_ = RP; Right_Prec_ = true; };
    
    void SetParameterList(const RefCountPtr<ParameterList> &PL) { PL_ = PL; };

    void SetBlockSize(int blocksize) { default_blocksize_ = blocksize; blocksize_ = blocksize; };
    

    void SetCurrLSVec();
    

    void AssertSymmetric(){ operatorSymmetric_ = true; };
    

    void SolutionUpdated( const MV* SolnUpdate = 0,
                          ScalarType scale = Teuchos::ScalarTraits<ScalarType>::one() );
    
    

    

    void Reset( const RefCountPtr<MV> &newX = null, const RefCountPtr<const MV> &newB = null );
    

    
    RefCountPtr<const OP> GetOperator() const { return(A_); };
    
    RefCountPtr<MV> GetLHS() const { return(X_); };
    
    RefCountPtr<const MV> GetRHS() const { return(B_); };
    

    const MV& GetInitResVec();
    

    const MV& GetCurrResVec( const MV* CurrSoln = 0 );
    

    RefCountPtr<MV> GetCurrLHSVec();
    

    RefCountPtr<MV> GetCurrRHSVec();
    
    RefCountPtr<const OP> GetLeftPrec() const { return(LP_); };
    
    RefCountPtr<const OP> GetRightPrec() const { return(RP_); };
    
    RefCountPtr<ParameterList> GetParameterList() const { return(PL_); };

    int GetBlockSize() const { return( default_blocksize_ ); };
    

    int GetCurrBlockSize() const { return( blocksize_ ); };


    int GetNumToSolve() const { return( num_to_solve_ ); };
    

    int GetRHSIndex() const { return( rhs_index_ ); };
    

    bool IsSolutionUpdated() const { return(solutionUpdated_); };
    
    bool IsOperatorSymmetric() const { return(operatorSymmetric_); };
    
    

    

    ReturnType Apply( const MV& x, MV& y );
    

    ReturnType ApplyOp( const MV& x, MV& y );
    

    ReturnType ApplyLeftPrec( const MV& x, MV& y );
    

    ReturnType ApplyRightPrec( const MV& x, MV& y );
    

    ReturnType ComputeResVec( MV* R, const MV* X = 0, const MV* B = 0 );
    
    
  private:
    
    void SetUpBlocks();
    
    RefCountPtr<const OP> A_;
    
    RefCountPtr<MV> X_;
    
    RefCountPtr<MV> CurX_;
    
    RefCountPtr<const MV> B_;
    
    RefCountPtr<MV> CurB_;
    
    RefCountPtr<MV> R_;
    
    RefCountPtr<MV> R0_;
    
    RefCountPtr<const OP> LP_;  
    
    RefCountPtr<const OP> RP_;
    
    RefCountPtr<ParameterList> PL_;

    int default_blocksize_;
    
    int blocksize_;

    int num_to_solve_;
    
    int rhs_index_;
    
    bool Left_Prec_;
    bool Right_Prec_;
    bool Left_Scale_;
    bool Right_Scale_;
    bool operatorSymmetric_;
    bool solutionUpdated_;    
    bool solutionFinal_;
    bool initresidsComputed_;
    
    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) : 
    default_blocksize_(1),
    blocksize_(1),
    num_to_solve_(0),
    rhs_index_(0),  
    Left_Prec_(false),
    Right_Prec_(false),
    Left_Scale_(false),
    Right_Scale_(false),
    operatorSymmetric_(false),
    solutionUpdated_(false),
    solutionFinal_(true),
    initresidsComputed_(false)
  {
  }
  
  template <class ScalarType, class MV, class OP>
  LinearProblem<ScalarType,MV,OP>::LinearProblem(const RefCountPtr<const OP> &A, 
             const RefCountPtr<MV> &X, 
             const RefCountPtr<const MV> &B
             ) :
    A_(A),
    X_(X),
    B_(B),
    default_blocksize_(1),
    blocksize_(1),
    num_to_solve_(1),
    rhs_index_(0),
    Left_Prec_(false),
    Right_Prec_(false),
    Left_Scale_(false),
    Right_Scale_(false),
    operatorSymmetric_(false),
    solutionUpdated_(false),
    solutionFinal_(true),
    initresidsComputed_(false)
  {
    R0_ = MVT::Clone( *X_, MVT::GetNumberVecs( *X_ ) );
  }
  
  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_),
    R_(Problem.R_),
    R0_(Problem.R0_),
    LP_(Problem.LP_),
    RP_(Problem.RP_),
    PL_(Problem.PL_),
    default_blocksize_(Problem.default_blocksize_),
    blocksize_(Problem.blocksize_),
    num_to_solve_(Problem.num_to_solve_),
    rhs_index_(Problem.rhs_index_),
    Left_Prec_(Problem.Left_Prec_),
    Right_Prec_(Problem.Right_Prec_),
    Left_Scale_(Problem.Left_Scale_),
    Right_Scale_(Problem.Right_Scale_),
    operatorSymmetric_(Problem.operatorSymmetric_),
    solutionUpdated_(Problem.solutionUpdated_),
    solutionFinal_(Problem.solutionFinal_),
    initresidsComputed_(Problem.initresidsComputed_)
  {
  }
  
  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>::SetUpBlocks()
  {
    // Compute the new block linear system.
    // ( first clean up old linear system )
    if (CurB_.get()) CurB_ = null;
    if (CurX_.get()) CurX_ = null;
    if (R_.get()) R_ = null;
    //
    // Determine how many linear systems are left to solve for and populate LHS and RHS vector.
    // If the number of linear systems left are less than the current blocksize, then
    // we create a multivector and copy the left over LHS and RHS vectors into them.
    // The rest of the multivector is populated with random vectors (RHS) or zero vectors (LHS).
    //
    num_to_solve_ = MVT::GetNumberVecs(*X_) - rhs_index_;
    //
    // Return the NULL pointer if we don't have any more systems to solve for.
    if ( num_to_solve_ <= 0 ) { return; }  
    //
    int i;
    std::vector<int> index( num_to_solve_ );
    for ( i=0; i<num_to_solve_; i++ ) { index[i] = rhs_index_ + i; }
    //
/*    if ( num_to_solve_ < default_blocksize_ )
      blocksize_ = num_to_solve_;
    else
      blocksize_ = default_blocksize_;
*/
    //
    if ( num_to_solve_ < blocksize_ ) 
      {
  std::vector<int> index2(num_to_solve_);
  for (i=0; i<num_to_solve_; i++) {
    index2[i] = i;
  }
  //
  // First create multivectors of blocksize and 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_);
  R_ = MVT::Clone( *X_, blocksize_);
  //
  RefCountPtr<const MV> tptr = MVT::CloneView( *B_, index );
  MVT::SetBlock( *tptr, index2, *CurB_ );
  //
  RefCountPtr<MV> tptr2 = MVT::CloneView( *X_, index );
  MVT::SetBlock( *tptr2, index2, *CurX_ );
      } else { 
  //
  // If the number of linear systems left are more than or equal to the current blocksize, then
  // we create a view into the LHS and RHS.
  //
  num_to_solve_ = blocksize_;
  index.resize( num_to_solve_ );
  for ( i=0; i<num_to_solve_; i++ ) { index[i] = rhs_index_ + i; }
  CurX_ = MVT::CloneView( *X_, index );
  CurB_ = rcp_const_cast<MV>(MVT::CloneView( *B_, index ));
  R_ = MVT::Clone( *X_, num_to_solve_ );
  //
      }
    //
    // Compute the current residual.
    // 
    if (R_.get()) {
      OPT::Apply( *A_, *CurX_, *R_ );
      MVT::MvAddMv( 1.0, *CurB_, -1.0, *R_, *R_ );
      solutionUpdated_ = false;
    }
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::SetLHS(const RefCountPtr<MV> &X)
  {
    X_ = X; 
    R0_ = MVT::Clone( *X_, MVT::GetNumberVecs( *X_ ) ); 
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::SetCurrLSVec() 
  { 
    int i;
    //
    // We only need to copy the solutions back if the linear systems of
    // interest are less than the block size.
    //
    if (num_to_solve_ < blocksize_) {
      //
      std::vector<int> index( num_to_solve_ );
      //
      RefCountPtr<MV> tptr;
      //
      // Get a view of the current solutions and correction vector.
      //
      for (i=0; i<num_to_solve_; i++) { 
  index[i] = i; 
      }
      tptr = MVT::CloneView( *CurX_, index );
      //
      // Copy the correction vector to the solution vector.
      //
      for (i=0; i<num_to_solve_; i++) { 
  index[i] = rhs_index_ + i; 
      }
      MVT::SetBlock( *tptr, index, *X_ );
    }
    //
    // Get the linear problem ready to determine the next linear system.
    //
    solutionFinal_ = true; 
    rhs_index_ += num_to_solve_; 
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::SolutionUpdated( const MV* SolnUpdate, ScalarType scale )
  { 
    if (SolnUpdate) {
      if (Right_Prec_) {
  //
  // Apply the right preconditioner before computing the current solution.
  RefCountPtr<MV> TrueUpdate = MVT::Clone( *SolnUpdate, MVT::GetNumberVecs( *SolnUpdate ) );
  OPT::Apply( *RP_, *SolnUpdate, *TrueUpdate ); 
  MVT::MvAddMv( 1.0, *CurX_, scale, *TrueUpdate, *CurX_ ); 
      } else {
  MVT::MvAddMv( 1.0, *CurX_, scale, *SolnUpdate, *CurX_ ); 
      }
    }
    solutionUpdated_ = true; 
  }
  
  template <class ScalarType, class MV, class OP>
  void LinearProblem<ScalarType,MV,OP>::Reset( const RefCountPtr<MV> &newX, const RefCountPtr<const MV> &newB )
  {
    solutionUpdated_ = false;
    solutionFinal_ = true;
    initresidsComputed_ = false;
    rhs_index_ = 0;
    
    X_ = newX;
    B_ = newB;
    GetInitResVec();
  }
  
  template <class ScalarType, class MV, class OP>
  const MV& LinearProblem<ScalarType,MV,OP>::GetInitResVec() 
  {
    // Compute the initial residual if it hasn't been computed
    // and all the components of the linear system are there.
    // The left preconditioner will be applied if it exists, resulting
    // in a preconditioned residual.
    if (!initresidsComputed_ && A_.get() && X_.get() && B_.get()) 
      {
  if (R0_.get()) R0_ = null;
  R0_ = MVT::Clone( *X_, MVT::GetNumberVecs( *X_ ) );
  OPT::Apply( *A_, *X_, *R0_ );
  MVT::MvAddMv( 1.0, *B_, -1.0, *R0_, *R0_ );
  initresidsComputed_ = true;
      }
    return (*R0_);
  }
  
  template <class ScalarType, class MV, class OP>
  const MV& LinearProblem<ScalarType,MV,OP>::GetCurrResVec( const MV* CurrSoln ) 
  {
    // Compute the residual of the current linear system.
    // This should be used if the solution has been updated.
    // Alternatively, if the current solution has been computed by GMRES
    // this can be passed in and the current residual will be updated using
    // it.
    //
    if (solutionUpdated_) 
      {
  OPT::Apply( *A_, *GetCurrLHSVec(), *R_ );
  MVT::MvAddMv( 1.0, *GetCurrRHSVec(), -1.0, *R_, *R_ ); 
  solutionUpdated_ = false;
      }
    else if (CurrSoln) 
      {
  OPT::Apply( *A_, *CurrSoln, *R_ );
  MVT::MvAddMv( 1.0, *GetCurrRHSVec(), -1.0, *R_, *R_ ); 
      }
    return (*R_);
  }
  
  template <class ScalarType, class MV, class OP>
  RefCountPtr<MV> LinearProblem<ScalarType,MV,OP>::GetCurrLHSVec()
  {
    if (solutionFinal_) {
      solutionFinal_ = false; // make sure we don't populate the current linear system again.
      SetUpBlocks();
    }
    return CurX_; 
  }
  
  template <class ScalarType, class MV, class OP>
  RefCountPtr<MV> LinearProblem<ScalarType,MV,OP>::GetCurrRHSVec()
  {
    if (solutionFinal_) {
      solutionFinal_ = false; // make sure we don't populate the current linear system again.
      SetUpBlocks();
    }
    return CurB_;
  }
  
  template <class ScalarType, class MV, class OP>
  ReturnType LinearProblem<ScalarType,MV,OP>::Apply( const MV& x, MV& y )
  {
    RefCountPtr<MV> ytemp = MVT::Clone( y, MVT::GetNumberVecs( y ) );
    //
    // No preconditioning.
    // 
    if (!Left_Prec_ && !Right_Prec_){ OPT::Apply( *A_, x, y );}
    //
    // Preconditioning is being done on both sides
    //
    else if( Left_Prec_ && Right_Prec_ ) 
      {
  OPT::Apply( *RP_, x, y );   
  OPT::Apply( *A_, y, *ytemp );
  OPT::Apply( *LP_, *ytemp, y );
      }
    //
    // Preconditioning is only being done on the left side
    //
    else if( Left_Prec_ ) 
      {
  OPT::Apply( *A_, x, *ytemp );
  OPT::Apply( *LP_, *ytemp, y );
      }
    //
    // Preconditioning is only being done on the right side
    //
    else 
      {
  OPT::Apply( *RP_, x, *ytemp );
  OPT::Apply( *A_, *ytemp, y );
      }  
    return Ok;
  }
  
  template <class ScalarType, class MV, class OP>
  ReturnType LinearProblem<ScalarType,MV,OP>::ApplyOp( const MV& x, MV& y )
  {
    if (A_.get())
      return ( OPT::Apply( *A_,x, y) );   
    else
      return Undefined;
  }
  
  template <class ScalarType, class MV, class OP>
  ReturnType LinearProblem<ScalarType,MV,OP>::ApplyLeftPrec( const MV& x, MV& y )
  {
    if (Left_Prec_)
      return ( OPT::Apply( *LP_,x, y) );
    else 
      return Undefined;
  }
  
  template <class ScalarType, class MV, class OP>
  ReturnType LinearProblem<ScalarType,MV,OP>::ApplyRightPrec( const MV& x, MV& y )
  {
    if (Right_Prec_)
      return ( OPT::Apply( *RP_,x, y) );
    else
      return Undefined;
  }
  
  template <class ScalarType, class MV, class OP>
  ReturnType LinearProblem<ScalarType,MV,OP>::ComputeResVec( MV* R, const MV* X, const MV* B )
  {
    if (X && B) // The entries are specified, so compute the residual of Op(A)X = B
      {
  if (Left_Prec_)
    {
      RefCountPtr<MV> R_temp = MVT::Clone( *X, MVT::GetNumberVecs( *X ) );
      OPT::Apply( *A_, *X, *R_temp );
      MVT::MvAddMv( -1.0, *R_temp, 1.0, *B, *R_temp );
      OPT::Apply( *LP_, *R_temp, *R );
    }
  else 
    {
      OPT::Apply( *A_, *X, *R );
      MVT::MvAddMv( -1.0, *R, 1.0, *B, *R );
    }
      }
    else { 
      // One of the entries is not specified, so just use the linear system information we have.
      // Later we may want to check to see which multivec is not specified, and use what is specified.
      if (Left_Prec_)
  {
    RefCountPtr<MV> R_temp = MVT::Clone( *X_, MVT::GetNumberVecs( *X_ ) );
    OPT::Apply( *A_, *X_, *R_temp );
    MVT::MvAddMv( -1.0, *R_temp, 1.0, *B_, *R_temp );
    OPT::Apply( *LP_, *R_temp, *R );
  }
      else 
  {
    OPT::Apply( *A_, *X_, *R );
    MVT::MvAddMv( -1.0, *R, 1.0, *B_, *R );
  }
    }    
    return Ok;
  }
  
} // end Belos namespace

#endif /* BELOS_LINEAR_PROBLEM_HPP */

---