Belos_StatusTest_GenResNorm_MP_Vector.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_STATUS_TEST_GEN_RESNORM_MP_VECTOR_H
#define BELOS_STATUS_TEST_GEN_RESNORM_MP_VECTOR_H
 
#include "Stokhos_Sacado_Kokkos_MP_Vector.hpp"
#include "BelosStatusTestGenResNorm.hpp"
 
namespace Belos {
 
template <class Storage, class MV, class OP>
class StatusTestGenResNorm<Sacado::MP::Vector<Storage>, MV, OP> :
    public StatusTestResNorm<Sacado::MP::Vector<Storage>,MV,OP> {
 
 public:
 
  // Convenience typedefs
  typedef Sacado::MP::Vector<Storage> ScalarType;
  typedef Teuchos::ScalarTraits<ScalarType> SCT;
  typedef typename SCT::magnitudeType MagnitudeType;
  typedef MultiVecTraits<ScalarType,MV>  MVT;
 
 
 
  enum ResType {Implicit, 
                Explicit  
  };
 
 
 
 
  StatusTestGenResNorm( MagnitudeType Tolerance, int quorum = -1, bool showMaxResNormOnly = false );
 
  virtual ~StatusTestGenResNorm();
 
 
 
 
  int defineResForm( ResType TypeOfResidual, NormType TypeOfNorm);
 
 
  int defineScaleForm( ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one());
 
 
  int setTolerance(MagnitudeType tolerance) {tolerance_ = tolerance; return(0);}
 
  int setQuorum(int quorum) {quorum_ = quorum; return(0);}
 
  int setShowMaxResNormOnly(bool showMaxResNormOnly) {showMaxResNormOnly_ = showMaxResNormOnly; return(0);}
 
 
 
 
  StatusType checkStatus(Iteration<ScalarType,MV,OP>* iSolver);
 
  StatusType getStatus() const {return(status_);};
 
 
 
  void reset();
 
 
 
 
  void print(std::ostream& os, int indent = 0) const;
 
  void printStatus(std::ostream& os, StatusType type) const;
 
 
 
  Teuchos::RCP<MV> getSolution() { if (restype_==Implicit) { return Teuchos::null; } else { return curSoln_; } }
 
  int getQuorum() const { return quorum_; }
 
  bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }
 
  std::vector<int> convIndices() { return ind_; }
 
  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_);};
 
  bool getLOADetected() const { return false; }
 
  const std::vector<int> getEnsembleIterations() const { return ensemble_iterations; }
 
 
 
 
  StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);
 
 
  std::string description() const
  {
    std::ostringstream oss;
    oss << "Belos::StatusTestGenResNorm<>: " << resFormStr();
    oss << ", tol = " << tolerance_;
    return oss.str();
  }
 
 protected:
 
 private:
 
 
 
  std::string resFormStr() const
  {
    std::ostringstream oss;
    oss << "(";
    oss << ((resnormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
    oss << ((restype_==Explicit) ? " Exp" : " Imp");
    oss << " Res Vec) ";
 
    // If there is no residual scaling, return current string.
    if (scaletype_!=None)
    {
      // Insert division sign.
      oss << "/ ";
 
      // Determine output string for scaling, if there is any.
      if (scaletype_==UserProvided)
        oss << " (User Scale)";
      else {
        oss << "(";
        oss << ((scalenormtype_==OneNorm) ? "1-Norm" : (resnormtype_==TwoNorm) ? "2-Norm" : "Inf-Norm");
        if (scaletype_==NormOfInitRes || scaletype_==NormOfFullInitRes || scaletype_==NormOfFullScaledInitRes)
          oss << " Res0";
        else if (scaletype_==NormOfPrecInitRes || scaletype_==NormOfFullPrecInitRes || scaletype_==NormOfFullScaledPrecInitRes)
          oss << " Prec Res0";
        else
          oss << " RHS ";
        oss << ")";
      }
    }
 
    return oss.str();
  }
 
 
 
 
  MagnitudeType tolerance_;
 
  int quorum_;
 
  bool showMaxResNormOnly_;
 
  ResType restype_;
 
  NormType resnormtype_;
 
  ScaleType scaletype_;
 
  NormType scalenormtype_;
 
  MagnitudeType scalevalue_;
 
  std::vector<MagnitudeType> scalevector_;
 
  std::vector<MagnitudeType> resvector_;
 
  std::vector<MagnitudeType> testvector_;
 
  std::vector<int> ind_;
 
  Teuchos::RCP<MV> curSoln_;
 
  StatusType status_;
 
  int curBlksz_;
 
  int curNumRHS_;
 
  std::vector<int> curLSIdx_;
 
  int curLSNum_;
 
  int numrhs_;
 
  bool firstcallCheckStatus_;
 
  bool firstcallDefineResForm_;
 
  bool firstcallDefineScaleForm_;
 
  std::vector<int> ensemble_converged;
 
  std::vector<int> ensemble_iterations;
 
 
};
 
template <class StorageType, class MV, class OP>
StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::
StatusTestGenResNorm (MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly)
  : tolerance_(Tolerance),
    quorum_(quorum),
    showMaxResNormOnly_(showMaxResNormOnly),
    restype_(Implicit),
    resnormtype_(TwoNorm),
    scaletype_(NormOfInitRes),
    scalenormtype_(TwoNorm),
    scalevalue_(Teuchos::ScalarTraits<MagnitudeType>::one ()),
    status_(Undefined),
    curBlksz_(0),
    curNumRHS_(0),
    curLSNum_(0),
    numrhs_(0),
    firstcallCheckStatus_(true),
    firstcallDefineResForm_(true),
    firstcallDefineScaleForm_(true),
    ensemble_converged(StorageType::static_size, 0),
    ensemble_iterations(StorageType::static_size, 0)
{
  // This constructor will compute the residual ||r_i||/||r0_i|| <= tolerance using the 2-norm of
  // the implicit residual std::vector.
}
 
template <class StorageType, class MV, class OP>
StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::~StatusTestGenResNorm()
{}
 
template <class StorageType, class MV, class OP>
void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::reset()
{
  status_ = Undefined;
  curBlksz_ = 0;
  curLSNum_ = 0;
  curLSIdx_.resize(0);
  numrhs_ = 0;
  ind_.resize(0);
  firstcallCheckStatus_ = true;
  curSoln_ = Teuchos::null;
  ensemble_converged = std::vector<int>(StorageType::static_size, 0);
  ensemble_iterations = std::vector<int>(StorageType::static_size, 0);
}
 
template <class StorageType, class MV, class OP>
int StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::defineResForm( ResType TypeOfResidual, NormType TypeOfNorm )
{
  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
        "StatusTestGenResNorm::defineResForm(): The residual form has already been defined.");
  firstcallDefineResForm_ = false;
 
  restype_ = TypeOfResidual;
  resnormtype_ = TypeOfNorm;
 
  return(0);
}
 
template <class StorageType, class MV, class OP>
int StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
                                                         MagnitudeType ScaleValue )
{
  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
        "StatusTestGenResNorm::defineScaleForm(): The scaling type has already been defined.");
  firstcallDefineScaleForm_ = false;
 
  scaletype_ = TypeOfScaling;
  scalenormtype_ = TypeOfNorm;
  scalevalue_ = ScaleValue;
 
  return(0);
}
 
template <class StorageType, class MV, class OP>
StatusType StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::checkStatus( Iteration<ScalarType,MV,OP>* iSolver )
{
  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
  const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
  // 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 std::vector
  //
  if ( curLSNum_ != lp.getLSNumber() ) {
    //
    // We have moved on to the next rhs block
    //
    curLSNum_ = lp.getLSNumber();
    curLSIdx_ = lp.getLSIndex();
    curBlksz_ = (int)curLSIdx_.size();
    int validLS = 0;
    for (int i=0; i<curBlksz_; ++i) {
      if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
        validLS++;
    }
    curNumRHS_ = validLS;
    curSoln_ = Teuchos::null;
    //
  } else {
    //
    // We are in the same rhs block, return if we are converged
    //
    if (status_==Passed) { 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( curBlksz_ );
    Teuchos::RCP<const MV> residMV = iSolver->getNativeResiduals( &tmp_resvector );
    if ( residMV != Teuchos::null ) {
      tmp_resvector.resize( MVT::GetNumberVecs( *residMV ) );
      MVT::MvNorm( *residMV, tmp_resvector, resnormtype_ );
      typename std::vector<int>::iterator p = curLSIdx_.begin();
      for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
        // Check if this index is valid
        if (*p != -1)
          resvector_[*p] = tmp_resvector[i];
      }
    } else {
      typename std::vector<int>::iterator p = curLSIdx_.begin();
      for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
        // Check if this index is valid
        if (*p != -1)
          resvector_[*p] = tmp_resvector[i];
      }
    }
  }
  else if (restype_==Explicit) {
    //
    // Request the true residual for this block of right-hand sides.
    //
    Teuchos::RCP<MV> cur_update = iSolver->getCurrentUpdate();
    curSoln_ = lp.updateSolution( cur_update );
    Teuchos::RCP<MV> cur_res = MVT::Clone( *curSoln_, MVT::GetNumberVecs( *curSoln_ ) );
    lp.computeCurrResVec( &*cur_res, &*curSoln_ );
    std::vector<MagnitudeType> tmp_resvector( MVT::GetNumberVecs( *cur_res ) );
    MVT::MvNorm( *cur_res, tmp_resvector, resnormtype_ );
    typename std::vector<int>::iterator p = curLSIdx_.begin();
    for (int i=0; p<curLSIdx_.end(); ++p, ++i) {
      // Check if this index is valid
      if (*p != -1)
        resvector_[*p] = 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 ) {
    typename std::vector<int>::iterator p = curLSIdx_.begin();
    for (; p<curLSIdx_.end(); ++p) {
      // Check if this index is valid
      if (*p != -1) {
        // Scale the std::vector accordingly
        if ( scalevector_[ *p ] != zero ) {
          // Don't intentionally divide by zero.
          testvector_[ *p ] = resvector_[ *p ] / scalevector_[ *p ] / scalevalue_;
        } else {
          testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
        }
      }
    }
  }
  else {
    typename std::vector<int>::iterator p = curLSIdx_.begin();
    for (; p<curLSIdx_.end(); ++p) {
      // Check if this index is valid
      if (*p != -1)
        testvector_[ *p ] = resvector_[ *p ] / scalevalue_;
    }
  }
 
  // Check status of new linear system residuals and see if we have the quorum.
  int have = 0;
  ind_.resize( curLSIdx_.size() );
  typename std::vector<int>::iterator p = curLSIdx_.begin();
  for (; p<curLSIdx_.end(); ++p) {
    // Check if this index is valid
    if (*p != -1) {
      // Check if any of the residuals are larger than the tolerance.
      const int ensemble_size = StorageType::static_size;
      bool all_converged = true;
      for (int i=0; i<ensemble_size; ++i) {
        if (!ensemble_converged[i]) {
          if (testvector_[ *p ].coeff(i) > tolerance_) {
            ++ensemble_iterations[i];
            all_converged = false;
          }
          else if (testvector_[ *p ].coeff(i) <= tolerance_) {
            ensemble_converged[i] = 1;
          }
          else {
            // Throw an std::exception if a NaN is found.
            status_ = Failed;
            TEUCHOS_TEST_FOR_EXCEPTION(true,StatusTestError,"StatusTestGenResNorm::checkStatus(): NaN has been detected.");
          }
        }
      }
      if (all_converged) {
        ind_[have] = *p;
        have++;
      }
    }
  }
  ind_.resize(have);
  int need = (quorum_ == -1) ? curNumRHS_: quorum_;
  status_ = (have >= need) ? Passed : Failed;
 
  // Return the current status
  return status_;
}
 
template <class StorageType, class MV, class OP>
void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::print(std::ostream& os, int indent) const
{
  for (int j = 0; j < indent; j ++)
    os << ' ';
  printStatus(os, status_);
  os << resFormStr();
  if (status_==Undefined)
    os << ", tol = " << tolerance_ << std::endl;
  else {
    os << std::endl;
    if(showMaxResNormOnly_ && curBlksz_ > 1) {
      const MagnitudeType maxRelRes = *std::max_element(
        testvector_.begin()+curLSIdx_[0],testvector_.begin()+curLSIdx_[curBlksz_-1]
        );
      for (int j = 0; j < indent + 13; j ++)
        os << ' ';
      os << "max{residual["<<curLSIdx_[0]<<"..."<<curLSIdx_[curBlksz_-1]<<"]} = " << maxRelRes
         << ( maxRelRes <= tolerance_ ? " <= " : " > " ) << tolerance_ << std::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_ << std::endl;
      }
    }
  }
  os << std::endl;
}
 
template <class StorageType, class MV, class OP>
void StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::printStatus(std::ostream& os, StatusType type) const
{
  os << std::left << std::setw(13) << std::setfill('.');
  switch (type) {
  case  Passed:
    os << "Converged";
    break;
  case  Failed:
    os << "Unconverged";
    break;
  case  Undefined:
  default:
    os << "**";
    break;
  }
  os << std::left << std::setfill(' ');
    return;
}
 
template <class StorageType, class MV, class OP>
StatusType StatusTestGenResNorm<Sacado::MP::Vector<StorageType>,MV,OP>::firstCallCheckStatusSetup( Iteration<ScalarType,MV,OP>* iSolver )
{
  int i;
  MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
  MagnitudeType one = Teuchos::ScalarTraits<MagnitudeType>::one();
  const LinearProblem<ScalarType,MV,OP>& lp = iSolver->getProblem();
  // Compute scaling term (done once for each block that's being solved)
  if (firstcallCheckStatus_) {
    //
    // Get some current solver information.
    //
    firstcallCheckStatus_ = false;
 
    if (scaletype_== NormOfRHS) {
      Teuchos::RCP<const MV> rhs = lp.getRHS();
      numrhs_ = MVT::GetNumberVecs( *rhs );
      scalevector_.resize( numrhs_ );
      MVT::MvNorm( *rhs, scalevector_, scalenormtype_ );
    }
    else if (scaletype_==NormOfInitRes || scaletype_==NormOfFullInitRes || scaletype_==NormOfFullScaledInitRes) {
      Teuchos::RCP<const MV> init_res = lp.getInitResVec();
      numrhs_ = MVT::GetNumberVecs( *init_res );
      scalevector_.resize( numrhs_ );
      MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
    }
    else if (scaletype_==NormOfPrecInitRes || scaletype_==NormOfFullPrecInitRes || scaletype_==NormOfFullScaledPrecInitRes) {
      Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
      numrhs_ = MVT::GetNumberVecs( *init_res );
      scalevector_.resize( numrhs_ );
      MVT::MvNorm( *init_res, scalevector_, scalenormtype_ );
    }
    else {
      numrhs_ = MVT::GetNumberVecs( *(lp.getRHS()) );
    }
 
    resvector_.resize( numrhs_ );
    testvector_.resize( numrhs_ );
 
    curLSNum_ = lp.getLSNumber();
    curLSIdx_ = lp.getLSIndex();
    curBlksz_ = (int)curLSIdx_.size();
    int validLS = 0;
    for (i=0; i<curBlksz_; ++i) {
      if (curLSIdx_[i] > -1 && curLSIdx_[i] < numrhs_)
        validLS++;
    }
    curNumRHS_ = validLS;
    //
    // 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 Undefined;
}
 
} // end namespace Belos
 
#endif /* BELOS_STATUS_TEST_RESNORM_H */