BelosStatusTestResNorm.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_STATUS_TEST_RESNORM_H
#define BELOS_STATUS_TEST_RESNORM_H

#include "BelosStatusTest.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosMultiVecTraits.hpp"

namespace Belos {

template <class ScalarType, class MV, class OP>
class StatusTestResNorm: public StatusTest<ScalarType,MV,OP> {

 public:

  // Convenience typedefs
  typedef Teuchos::ScalarTraits<ScalarType> SCT;
  typedef typename SCT::magnitudeType MagnitudeType;
  typedef MultiVecTraits<ScalarType,MV>  MVT;



  enum ResType {Implicit, 
    Explicit  
  };
  enum ScaleType {NormOfRHS,     
      NormOfInitRes, 
      NormOfPrecInitRes, 
      None,          
      UserProvided   
  };




  StatusTestResNorm( MagnitudeType Tolerance, bool showMaxResNormOnly = false );

  virtual ~StatusTestResNorm();




  int DefineResForm( ResType TypeOfResidual, NormType TypeOfNorm);
  

  int DefineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());


  int ResetTolerance(MagnitudeType Tolerance) {tolerance_ = Tolerance; return(0);};





  StatusType CheckStatus(IterativeSolver<ScalarType,MV,OP>* iSolver);

  StatusType GetStatus() const {return(status_);};


 
  void Reset();





  bool ResidualVectorRequired() const { return(resvecrequired_); };




  ostream& Print(ostream& os, int indent = 0) const;



  MagnitudeType GetTolerance() const {return(tolerance_);};
  
  const std::vector<MagnitudeType>* GetTestValue() const {return(&testvector_);};

  const std::vector<MagnitudeType>* GetResNormValue() const {return(&resvector_);};

  const std::vector<MagnitudeType>* GetScaledNormValue() const {return(&scalevector_);};




  StatusType FirstcallCheckStatusSetup(IterativeSolver<ScalarType,MV,OP>* iSolver);

 protected:

 private:


  
  MagnitudeType tolerance_;

  bool showMaxResNormOnly_;
 
  ResType restype_;
  
  NormType resnormtype_;
  
  ScaleType scaletype_;
  
  NormType scalenormtype_;

  MagnitudeType scalevalue_;

  std::vector<MagnitudeType> scalevector_;
  
  std::vector<MagnitudeType> resvector_;

  std::vector<MagnitudeType> testvector_;
  
  StatusType status_;
  
  int cur_rhs_num_;

  int cur_blksz_;

  int cur_num_rhs_;

  int numrhs_;

  bool resvecrequired_;

  bool firstcallCheckStatus_;

  bool firstcallDefineResForm_;

  bool firstcallDefineScaleForm_;


};

template <class ScalarType, class MV, class OP>
StatusTestResNorm<ScalarType,MV,OP>::StatusTestResNorm( MagnitudeType Tolerance, bool showMaxResNormOnly )
  : tolerance_(Tolerance),
    showMaxResNormOnly_(showMaxResNormOnly),
    restype_(Implicit),
    resnormtype_(TwoNorm),  
    scaletype_(NormOfInitRes),
    scalenormtype_(TwoNorm),
    scalevalue_(1.0),
    status_(Unchecked),
    cur_rhs_num_(0),
    cur_blksz_(0),
    numrhs_(0),
    resvecrequired_(false),
    firstcallCheckStatus_(true),
    firstcallDefineResForm_(true),
    firstcallDefineScaleForm_(true)
{
  // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
  // the implicit residual vector.
}

template <class ScalarType, class MV, class OP>
StatusTestResNorm<ScalarType,MV,OP>::~StatusTestResNorm() 
{}

template <class ScalarType, class MV, class OP>
void StatusTestResNorm<ScalarType,MV,OP>::Reset() 
{
  status_ = Unchecked;
  cur_rhs_num_ = 0;
  cur_blksz_ = 0;
  numrhs_ = 0;
  firstcallCheckStatus_ = true;
}

template <class ScalarType, class MV, class OP>
int StatusTestResNorm<ScalarType,MV,OP>::DefineResForm( ResType TypeOfResidual, NormType TypeOfNorm )
{    
  assert( firstcallDefineResForm_ );
  firstcallDefineResForm_ = false;
    
  restype_ = TypeOfResidual;
  resnormtype_ = TypeOfNorm;
    
  // These conditions force the residual vector to be computed
  if (restype_==Explicit)
    resvecrequired_ = true;
    
  return(0);
}

template <class ScalarType, class MV, class OP> 
int StatusTestResNorm<ScalarType,MV,OP>::DefineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
                                                         MagnitudeType ScaleValue )
{
    
  assert( firstcallDefineScaleForm_ );
  firstcallDefineScaleForm_ = false;
    
  scaletype_ = TypeOfScaling;
  scalenormtype_ = TypeOfNorm;
  scalevalue_ = ScaleValue;
    
  return(0);
}

template <class ScalarType, class MV, class OP>
StatusType StatusTestResNorm<ScalarType,MV,OP>::CheckStatus( IterativeSolver<ScalarType,MV,OP>* iSolver )
{
  int i;
  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
  RefCountPtr<LinearProblem<ScalarType,MV,OP> > lp = iSolver->GetLinearProblem();
  // Compute scaling term (done once for each block that's being solved)
  if (firstcallCheckStatus_) {
    StatusType status = FirstcallCheckStatusSetup(iSolver);
    if(status==Failed) {
      status_ = Failed;
      return(status_);
    }
  }
  //
  // This section computes the norm of the residual vector
  //
  if ( cur_rhs_num_ != lp->GetRHSIndex() || cur_blksz_ != lp->GetCurrBlockSize() || cur_num_rhs_ != lp->GetNumToSolve() ) {
    //
    // We have moved on to the next rhs block
    //
    cur_rhs_num_ = lp->GetRHSIndex();
    cur_blksz_ = lp->GetCurrBlockSize();
    cur_num_rhs_ = lp->GetNumToSolve();
    //
  } else {
    //
    // We are in the same rhs block, return if we are converged
    //
    if (status_==Converged) { return status_; }
  }
  if (restype_==Implicit) {
    //
    // Get the native residual norms from the solver for this block of right-hand sides.
    // If the residual is returned in multivector form, use the resnormtype to compute the residual norms.
    // Otherwise the native residual is assumed to be stored in the resvector_.
    //
    std::vector<MagnitudeType> tmp_resvector( cur_blksz_ );
    RefCountPtr<const MV> residMV = iSolver->GetNativeResiduals( &tmp_resvector );     
    if ( residMV.get() != NULL ) { 
      tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
      MVT::MvNorm( *residMV, &tmp_resvector, resnormtype_ );    
      for (i=0; i<MVT::GetNumberVecs( *residMV ) && i<cur_num_rhs_; i++)
        resvector_[i+cur_rhs_num_] = tmp_resvector[i]; 
    } else {
      for (i=0; i<cur_num_rhs_; i++)
        resvector_[i+cur_rhs_num_] = tmp_resvector[i];
    }
  }
  else if (restype_==Explicit) {
    //
    // Request the true residual for this block of right-hand sides.
    // See if the linear problem manager has been updated before
    // asking for the true residual from the solver.
    //
    //
    if ( lp->IsSolutionUpdated() ) {
      const MV &cur_res = lp->GetCurrResVec();
      std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( cur_res ) );
      MVT::MvNorm( cur_res, &tmp_resvector, resnormtype_ );
      for (i=0; i<MVT::GetNumberVecs( cur_res ) && i<cur_num_rhs_; i++)
        resvector_[i+cur_rhs_num_] = tmp_resvector[i];
    } else {
      RefCountPtr<const MV> cur_soln = iSolver->GetCurrentSoln();
      const MV &cur_res = lp->GetCurrResVec( &*cur_soln );
      std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( cur_res ) );
      MVT::MvNorm( cur_res, &tmp_resvector, resnormtype_ );
      for (i=0; i<MVT::GetNumberVecs( cur_res ) && i<cur_num_rhs_; i++)
        resvector_[i+cur_rhs_num_] = tmp_resvector[i];      
    }
  }
  //
  // Compute the new linear system residuals for testing.
  // (if any of them don't meet the tolerance or are NaN, then we exit with that status)
  //
  if ( scalevector_.size() > 0 ) {
    for (i = cur_rhs_num_; i < (cur_rhs_num_ + cur_num_rhs_); i++) {
     
      // Scale the vector accordingly
      if ( scalevector_[i] != zero ) {
        // Don't intentionally divide by zero.
        testvector_[ i ] = resvector_[ i ] / scalevector_[ i ] / scalevalue_;
      } else {
        testvector_[ i ] = resvector_[ i ] / scalevalue_;
      }
    }
  }
  else {
    for (i = cur_rhs_num_; i < (cur_rhs_num_ + cur_num_rhs_); i++)
      testvector_[ i ] = resvector_[ i ] / scalevalue_;
  } 

  // Check status of new linear system residuals
  status_ = Converged; // This may be set to unconverged or NaN.
  for (i = cur_rhs_num_; i < (cur_rhs_num_ + cur_num_rhs_); i++) {
    // Check if any of the residuals are larger than the tolerance.
    if (testvector_[ i ] > tolerance_) {
      status_ = Unconverged;
      return(status_);
    } else if (testvector_[ i ] <= tolerance_) { 
      // do nothing.
    } else {
      status_ = NaN;            
      return(status_); // Return immediately if we detect a NaN.
    }
  } 
  
  // Return the current status
  return status_;
}

template <class ScalarType, class MV, class OP>
ostream& StatusTestResNorm<ScalarType,MV,OP>::Print(ostream& os, int indent) const
{
  for (int j = 0; j < indent; j ++)
    os << ' ';
  this->PrintStatus(os, status_);
  os << "(";
  os << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
  os << ((restype_==Explicit) ? " Exp" : " Imp");
  os << " Res Vec) ";
  if (scaletype_!=None)
    os << "/ ";
  if (scaletype_==UserProvided)
    os << " (User Scale)";
  else {
    os << "(";
    os << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
    if (scaletype_==NormOfInitRes)
      os << " Res0";
    else if (scaletype_==NormOfPrecInitRes)
      os << " Prec Res0";
    else
      os << " RHS ";
    os << ")";
  }
  if (status_==Unchecked)
    os << " Unchecked ( tol = " << tolerance_ << " ) "<<endl;
  else {
    os << endl;
    if(showMaxResNormOnly_) {
      const MagnitudeType maxRelRes = *std::max_element(
        testvector_.begin()+cur_rhs_num_,testvector_.begin()+cur_rhs_num_+cur_num_rhs_
        );
      for (int j = 0; j < indent + 13; j ++)
        os << ' ';
      os << "max{residual["<<cur_rhs_num_<<"..."<<cur_rhs_num_+cur_num_rhs_-1<<"]} = " << maxRelRes
         << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << endl;
    }
    else {
      for ( int i=0; i<numrhs_; i++ ) {
        for (int j = 0; j < indent + 13; j ++)
          os << ' ';
        os << "residual [ " << i << " ] = " << testvector_[ i ];
        os << ((testvector_[i]<tolerance_) ? " < " : (testvector_[i]==tolerance_) ? " == " : (testvector_[i]>tolerance_) ? " > " : " "  ) << tolerance_ << endl;
      }
    }
  }
  os << endl;
  return os;
}

template <class ScalarType, class MV, class OP>
StatusType StatusTestResNorm<ScalarType,MV,OP>::FirstcallCheckStatusSetup( IterativeSolver<ScalarType,MV,OP>* iSolver )
{
  int i;
  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
  MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
  RefCountPtr<LinearProblem<ScalarType,MV,OP> > lp = iSolver->GetLinearProblem();
  // Compute scaling term (done once for each block that's being solved)
  if (firstcallCheckStatus_) {
    //
    // Get some current solver information.
    //
    firstcallCheckStatus_ = false;
    cur_rhs_num_ = lp->GetRHSIndex();
    cur_blksz_ = lp->GetCurrBlockSize();
    cur_num_rhs_ = lp->GetNumToSolve();
    //
    if (scaletype_== NormOfRHS) {
      const MV& rhs = *(lp->GetRHS());
      numrhs_ = MVT::GetNumberVecs( rhs );
      scalevector_.resize( numrhs_ );
      resvector_.resize( numrhs_ ); 
      testvector_.resize( numrhs_ );
      MVT::MvNorm( rhs, &scalevector_, scalenormtype_ );
    }
    else if (scaletype_==NormOfInitRes) {
      const MV &init_res = lp->GetInitResVec();
      numrhs_ = MVT::GetNumberVecs( init_res );
      scalevector_.resize( numrhs_ );
      resvector_.resize( numrhs_ ); 
      testvector_.resize( numrhs_ );
      MVT::MvNorm( init_res, &scalevector_, scalenormtype_ );
    }
    else if (scaletype_==NormOfPrecInitRes) {
      const MV& init_res = lp->GetInitResVec();
      numrhs_ = MVT::GetNumberVecs( init_res );
      scalevector_.resize( numrhs_ );
      resvector_.resize( numrhs_ ); 
      testvector_.resize( numrhs_ );
      RefCountPtr<MV> prec_init_res = MVT::Clone( init_res, numrhs_ );
      if (lp->ApplyLeftPrec( init_res, *prec_init_res ) != Undefined)
        MVT::MvNorm( *prec_init_res, &scalevector_, scalenormtype_ );
      else 
        MVT::MvNorm( init_res, &scalevector_, scalenormtype_ );
    }

    // Initialize the testvector.
    for (i=0; i<numrhs_; i++) { testvector_[i] = one; }

    // Return an error if the scaling is zero.
    if (scalevalue_ == zero) {
      return Failed;
    }
  }
  return Unchecked;
}

} // end namespace Belos

#endif /* BELOS_STATUS_TEST_RESNORM_H */

---