BelosStatusTestImpResNorm.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_IMPRESNORM_H
#define BELOS_STATUS_TEST_IMPRESNORM_H

#include "BelosStatusTestResNorm.hpp"
#include "BelosLinearProblem.hpp"
#include "BelosMultiVecTraits.hpp"
#include "Teuchos_as.hpp"

namespace Belos {

template <class ScalarType, class MV, class OP>
class StatusTestImpResNorm: public StatusTestResNorm<ScalarType,MV,OP> {
public:
  typedef typename Teuchos::ScalarTraits<ScalarType>::magnitudeType MagnitudeType;

private:

  typedef Teuchos::ScalarTraits<ScalarType> STS;
  typedef Teuchos::ScalarTraits<MagnitudeType> STM;
  typedef MultiVecTraits<ScalarType,MV> MVT;

public:


  StatusTestImpResNorm (MagnitudeType Tolerance,
                        int quorum = -1,
                        bool showMaxResNormOnly = false);

  virtual ~StatusTestImpResNorm();





  int defineResForm( 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() { return curSoln_; }

  int getQuorum() const { return quorum_; }

  bool getShowMaxResNormOnly() { return showMaxResNormOnly_; }

  std::vector<int> convIndices() { return ind_; }

  MagnitudeType getTolerance () const {
    return tolerance_;
  }

  MagnitudeType getCurrTolerance () const {
    return currTolerance_;
  }

  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 lossDetected_; }



  StatusType firstCallCheckStatusSetup(Iteration<ScalarType,MV,OP>* iSolver);


  std::string description() const
  {
    std::ostringstream oss;
    oss << "Belos::StatusTestImpResNorm<>: " << 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 << " 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)
          oss << " Res0";
        else if (scaletype_==NormOfPrecInitRes)
          oss << " Prec Res0";
        else
          oss << " RHS ";
        oss << ")";
      }
    }

    return oss.str();
  }



  MagnitudeType tolerance_, currTolerance_;

  int quorum_;

  bool showMaxResNormOnly_;

  NormType resnormtype_;

  ScaleType scaletype_;

  NormType scalenormtype_;

  MagnitudeType scalevalue_;

  std::vector<MagnitudeType> scalevector_;

  std::vector<MagnitudeType> resvector_;

  std::vector<MagnitudeType> testvector_;

  Teuchos::RCP<MV> curSoln_;

  std::vector<int> ind_;

  StatusType status_;

  int curBlksz_;

  int curNumRHS_;

  std::vector<int> curLSIdx_;

  int curLSNum_;

  int numrhs_;

  bool firstcallCheckStatus_;

  bool firstcallDefineResForm_;

  bool firstcallDefineScaleForm_;

  bool lossDetected_;

};

template <class ScalarType, class MV, class OP>
StatusTestImpResNorm<ScalarType,MV,OP>::
StatusTestImpResNorm (MagnitudeType Tolerance, int quorum, bool showMaxResNormOnly)
  : tolerance_(Tolerance),
    currTolerance_(Tolerance),
    quorum_(quorum),
    showMaxResNormOnly_(showMaxResNormOnly),
    resnormtype_(TwoNorm),
    scaletype_(NormOfInitRes),
    scalenormtype_(TwoNorm),
    scalevalue_(1.0),
    status_(Undefined),
    curBlksz_(0),
    curLSNum_(0),
    numrhs_(0),
    firstcallCheckStatus_(true),
    firstcallDefineResForm_(true),
    firstcallDefineScaleForm_(true),
    lossDetected_(false)
{
  // 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>
StatusTestImpResNorm<ScalarType,MV,OP>::~StatusTestImpResNorm()
{}

template <class ScalarType, class MV, class OP>
void StatusTestImpResNorm<ScalarType,MV,OP>::reset()
{
  status_ = Undefined;
  curBlksz_ = 0;
  curLSNum_ = 0;
  curLSIdx_.resize(0);
  numrhs_ = 0;
  ind_.resize(0);
  currTolerance_ = tolerance_;
  firstcallCheckStatus_ = true;
  lossDetected_ = false;
  curSoln_ = Teuchos::null;
}

template <class ScalarType, class MV, class OP>
int StatusTestImpResNorm<ScalarType,MV,OP>::defineResForm( NormType TypeOfNorm )
{
  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_==false,StatusTestError,
        "StatusTestResNorm::defineResForm(): The residual form has already been defined.");
  firstcallDefineResForm_ = false;

  resnormtype_ = TypeOfNorm;

  return(0);
}

template <class ScalarType, class MV, class OP>
int StatusTestImpResNorm<ScalarType,MV,OP>::defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm,
                                                         MagnitudeType ScaleValue )
{
  TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_==false,StatusTestError,
        "StatusTestResNorm::defineScaleForm(): The scaling type has already been defined.");
  firstcallDefineScaleForm_ = false;

  scaletype_ = TypeOfScaling;
  scalenormtype_ = TypeOfNorm;
  scalevalue_ = ScaleValue;

  return(0);
}

template <class ScalarType, class MV, class OP>
StatusType StatusTestImpResNorm<ScalarType,MV,OP>::
checkStatus (Iteration<ScalarType,MV,OP>* iSolver)
{
  using Teuchos::as;
  using Teuchos::RCP;

  const MagnitudeType zero = STM::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_;
    }
  }

  // mfh 23 Apr 2012: I don't know exactly what this code does.  It
  // has something to do with picking the block of right-hand sides
  // which we're currently checking.
  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_;
    }
  }

  //
  // Get the "native" residual norms from the solver for this block of
  // right-hand sides.  If the solver's getNativeResiduals() method
  // actually returns a multivector, compute the norms of the columns
  // of the multivector explicitly.  Otherwise, we assume that
  // resvector_ contains the norms.
  //
  // Note that "compute the norms explicitly" doesn't necessarily mean
  // the "explicit" residual norms (in the sense discussed in this
  // class' documentation).  These are just some vectors returned by
  // the solver.  Some Krylov methods, like CG, compute a residual
  // vector "recursively."  This is an "implicit residual" in the
  // sense of this class' documentation.  It equals the explicit
  // residual in exact arithmetic, but due to rounding error, it is
  // usually different than the explicit residual.
  //
  // FIXME (mfh 23 Apr 2012) This method does _not_ respect the
  // OrthoManager used by the solver.
  //
  std::vector<MagnitudeType> tmp_resvector( curBlksz_ );
  RCP<const MV> residMV = iSolver->getNativeResiduals (&tmp_resvector);
  if (! residMV.is_null ()) {
    // We got a multivector back.  Compute the norms explicitly.
    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];
      }
    }
  }
  //
  // Scale the unscaled residual norms we computed or obtained above.
  //
  if (scalevector_.size () > 0) {
    // There are per-vector scaling factors to apply.
    typename std::vector<int>::iterator p = curLSIdx_.begin();
    for (; p<curLSIdx_.end(); ++p) {
      // Check if this index is valid
      if (*p != -1) {
        // Scale the 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 { // There are no per-vector scaling factors.
    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_;
      }
    }
  }

  // Count how many scaled residual norms (in testvector_) pass, using
  // the current tolerance (currTolerance_) rather than the original
  // tolerance (tolerance_).  If at least quorum_ of them pass, we
  // have a quorum for the whole test to pass.
  //
  // We also check here whether any of the scaled residual norms is
  // NaN, and throw an exception in that case.
  int have = 0;
  ind_.resize( curLSIdx_.size() );
  std::vector<int> lclInd( curLSIdx_.size() );
  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) {
      if (testvector_[ *p ] > currTolerance_) {
        // The current residual norm doesn't pass.  Do nothing.
      } else if (testvector_[ *p ] <= currTolerance_) {
        ind_[have] = *p;
        lclInd[have] = i;
        have++; // Yay, the current residual norm passes!
      } else {
        // Throw an std::exception if the current residual norm is
        // NaN.  We know that it's NaN because it is not less than,
        // equal to, or greater than the current tolerance.  This is
        // only possible if either the residual norm or the current
        // tolerance is NaN; we assume the former.  We also mark the
        // test as failed, in case you want to catch the exception.
        status_ = Failed;
        TEUCHOS_TEST_FOR_EXCEPTION(true, StatusTestError, "Belos::"
          "StatusTestImpResNorm::checkStatus(): One or more of the current "
          "implicit residual norms is NaN.");
      }
    }
  }
  // "have" is the number of residual norms that passed.
  ind_.resize(have);
  lclInd.resize(have);

  // Now check the exact residuals
  if (have) { // At least one residual norm has converged.
    //
    // Compute the explicit residual norm(s) from the current solution update.
    //
    RCP<MV> cur_update = iSolver->getCurrentUpdate ();
    curSoln_ = lp.updateSolution (cur_update);
    RCP<MV> cur_res = MVT::Clone (*curSoln_, MVT::GetNumberVecs (*curSoln_));
    lp.computeCurrResVec (&*cur_res, &*curSoln_);
    tmp_resvector.resize (MVT::GetNumberVecs (*cur_res));
    std::vector<MagnitudeType> tmp_testvector (have);
    MVT::MvNorm (*cur_res, tmp_resvector, resnormtype_);

    // Scale the explicit residual norm(s), just like the implicit norm(s).
    if ( scalevector_.size() > 0 ) {
      for (int i=0; i<have; ++i) {
        // Scale the vector accordingly
        if ( scalevector_[ ind_[i] ] != zero ) {
          // Don't intentionally divide by zero.
          tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevector_[ ind_[i] ] / scalevalue_;
        } else {
          tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
        }
      }
    }
    else {
      for (int i=0; i<have; ++i) {
        tmp_testvector[ i ] = tmp_resvector[ lclInd[i] ] / scalevalue_;
      }
    }

    //
    // Check whether the explicit residual norms also pass the
    // convergence test.  If not, check whether we want to try
    // iterating a little more to force both implicit and explicit
    // residual norms to pass.
    //
    int have2 = 0;
    for (int i = 0; i < have; ++i) {
      // testvector_ contains the implicit (i.e., recursive, computed
      // by the algorithm) (possibly scaled) residuals.  All of these
      // pass the convergence test.
      //
      // tmp_testvector contains the explicit (i.e., ||B-AX||)
      // (possibly scaled) residuals.  We're checking whether these
      // pass as well.  The explicit residual norms only have to meet
      // the _original_ tolerance (tolerance_), not the current
      // tolerance (currTolerance_).
      if (tmp_testvector[i] <= tolerance_) {
        ind_[have2] = ind_[i];
        have2++; // This right-hand side has converged.
      }
      else {
        // Absolute difference between the current explicit and
        // implicit residual norm.
        const MagnitudeType diff = STM::magnitude (testvector_[ind_[i]] - tmp_testvector[i]);
        if (diff > currTolerance_) {
          // If the above difference is bigger than the current
          // convergence tolerance, report a loss of accuracy, but
          // mark this right-hand side converged.  (The latter tells
          // users not to iterate further on this right-hand side,
          // since it probably won't do much good.)  Note that the
          // current tolerance may have been changed by the branch
          // below in a previous call to this method.
          lossDetected_ = true;
          ind_[have2] = ind_[i];
          have2++; // Count this right-hand side as converged.
        }
        else {
          // Otherwise, the explicit and implicit residuals are pretty
          // close together, and the implicit residual has converged,
          // but the explicit residual hasn't converged.  Reduce the
          // convergence tolerance by some formula related to the
          // difference, and keep iterating.
          //
          // mfh 23 Apr 2012: I have no idea why the currTolerance_
          // update formula is done the way it's done below.  It
          // doesn't make sense to me.  It definitely makes the
          // current tolerance smaller, though, which is what we want.

          // We define these constants in this way, rather than simply
          // writing 1.5 resp. 0.1, to ensure no rounding error from
          // translating from float to MagnitudeType.  Remember that
          // 0.1 doesn't have an exact representation in binary finite
          // floating-point arithmetic.
          const MagnitudeType onePointFive = as<MagnitudeType>(3) / as<MagnitudeType> (2);
          const MagnitudeType oneTenth = STM::one () / as<MagnitudeType> (10);

          currTolerance_ = currTolerance_ - onePointFive * diff;
          while (currTolerance_ < STM::zero ()) {
            currTolerance_ += oneTenth * diff;
          }
        }
      }
    }
    have = have2;
    ind_.resize(have);
  }

  // Check whether we've met the quorum of vectors necessary for the
  // whole test to pass.
  int need = (quorum_ == -1) ? curNumRHS_: quorum_;
  status_ = (have >= need) ? Passed : Failed;

  // Return the current status
  return status_;
}

template <class ScalarType, class MV, class OP>
void StatusTestImpResNorm<ScalarType,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 ScalarType, class MV, class OP>
void StatusTestImpResNorm<ScalarType,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:
    if (lossDetected_)
      os << "Unconverged (LoA)";
    else
      os << "Unconverged";
    break;
  case  Undefined:
  default:
    os << "**";
    break;
  }
  os << std::left << std::setfill(' ');
    return;
}

template <class ScalarType, class MV, class OP>
StatusType StatusTestImpResNorm<ScalarType,MV,OP>::
firstCallCheckStatusSetup (Iteration<ScalarType,MV,OP>* iSolver)
{
  int i;
  const MagnitudeType zero = STM::zero ();
  const MagnitudeType one = STM::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) {
      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) {
      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_IMPRESNORM_H */