BelosXpetraStatusTestGenResSubNorm.hpp

Go to the documentation of this file.

// @HEADER
// *****************************************************************************
//        MueLu: A package for multigrid based preconditioning
//
// Copyright 2012 NTESS and the MueLu contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER

#ifndef BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_HPP
#define BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_HPP

#include "Xpetra_ConfigDefs.hpp"

#include "Xpetra_BlockedCrsMatrix.hpp"

#include "MueLu_Exceptions.hpp"

#include <BelosConfigDefs.hpp>
#include <BelosTypes.hpp>
#include <BelosOperatorT.hpp>
#include <BelosXpetraAdapterOperator.hpp>
#include <BelosStatusTestGenResSubNorm.hpp>

namespace Belos {

template <class Scalar, class LocalOrdinal, class GlobalOrdinal, class Node>
class StatusTestGenResSubNorm<Scalar, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > >
  : public StatusTestResNorm<Scalar, Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node>, Belos::OperatorT<Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> > > {
 public:
  // Convenience typedefs
  typedef Xpetra::MultiVector<Scalar, LocalOrdinal, GlobalOrdinal, Node> MV;
  typedef Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> BCRS;
  typedef Xpetra::MapExtractor<Scalar, LocalOrdinal, GlobalOrdinal, Node> ME;
  typedef Belos::OperatorT<MV> OP;

  typedef Teuchos::ScalarTraits<Scalar> SCT;
  typedef typename SCT::magnitudeType MagnitudeType;
  typedef MultiVecTraits<Scalar, MV> MVT;
  typedef OperatorTraits<Scalar, MV, OP> OT;



  StatusTestGenResSubNorm(MagnitudeType Tolerance, size_t subIdx, int quorum = -1, bool showMaxResNormOnly = false)
    : tolerance_(Tolerance)
    , subIdx_(subIdx)
    , quorum_(quorum)
    , showMaxResNormOnly_(showMaxResNormOnly)
    , 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)
    , mapExtractor_(Teuchos::null) {}

  virtual ~StatusTestGenResSubNorm(){};




  int defineResForm(NormType TypeOfNorm) {
    TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineResForm_ == false, StatusTestError,
                               "StatusTestGenResSubNorm::defineResForm(): The residual form has already been defined.");
    firstcallDefineResForm_ = false;

    resnormtype_ = TypeOfNorm;

    return (0);
  }


  int defineScaleForm(ScaleType TypeOfScaling, NormType TypeOfNorm, MagnitudeType ScaleValue = Teuchos::ScalarTraits<MagnitudeType>::one()) {
    TEUCHOS_TEST_FOR_EXCEPTION(firstcallDefineScaleForm_ == false, StatusTestError,
                               "StatusTestGenResSubNorm::defineScaleForm(): The scaling type has already been defined.");
    firstcallDefineScaleForm_ = false;

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

    return (0);
  }


  int setTolerance(MagnitudeType tolerance) {
    tolerance_ = tolerance;
    return (0);
  }


  int setSubIdx(size_t subIdx) {
    subIdx_ = subIdx;
    return (0);
  }

  int setQuorum(int quorum) {
    quorum_ = quorum;
    return (0);
  }

  int setShowMaxResNormOnly(bool showMaxResNormOnly) {
    showMaxResNormOnly_ = showMaxResNormOnly;
    return (0);
  }




  StatusType checkStatus(Iteration<Scalar, MV, OP>* iSolver) {
    MagnitudeType zero                      = Teuchos::ScalarTraits<MagnitudeType>::zero();
    const LinearProblem<Scalar, 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_;
      }
    }

    //
    // 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));
    MvSubNorm(*cur_res, subIdx_, 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 pp = curLSIdx_.begin();
      for (; pp < curLSIdx_.end(); ++pp) {
        // Check if this index is valid
        if (*pp != -1) {
          // Scale the std::vector accordingly
          if (scalevector_[*pp] != zero) {
            // Don't intentionally divide by zero.
            testvector_[*pp] = resvector_[*pp] / scalevector_[*pp] / scalevalue_;
          } else {
            testvector_[*pp] = resvector_[*pp] / scalevalue_;
          }
        }
      }
    } else {
      typename std::vector<int>::iterator pp = curLSIdx_.begin();
      for (; pp < curLSIdx_.end(); ++pp) {
        // Check if this index is valid
        if (*pp != -1)
          testvector_[*pp] = resvector_[*pp] / 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 p2 = curLSIdx_.begin();
    for (; p2 < curLSIdx_.end(); ++p2) {
      // Check if this index is valid
      if (*p2 != -1) {
        // Check if any of the residuals are larger than the tolerance.
        if (testvector_[*p2] > tolerance_) {
          // do nothing.
        } else if (testvector_[*p2] == Teuchos::ScalarTraits<Scalar>::magnitude(Teuchos::ScalarTraits<Scalar>::zero())) {
          reset();
        } else if (testvector_[*p2] <= tolerance_) {
          ind_[have] = *p2;
          have++;
        } else {
          // Throw an std::exception if a NaN is found.
          status_ = Failed;
          TEUCHOS_TEST_FOR_EXCEPTION(true, StatusTestError, "StatusTestGenResSubNorm::checkStatus(): NaN has been detected.");
        }
      }
    }
    ind_.resize(have);
    int need = (quorum_ == -1) ? curNumRHS_ : quorum_;
    status_  = (have >= need) ? Passed : Failed;
    // Return the current status
    return status_;
  }

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



  void reset() {
    status_   = Undefined;
    curBlksz_ = 0;
    curLSNum_ = 0;
    curLSIdx_.resize(0);
    numrhs_ = 0;
    ind_.resize(0);
    firstcallCheckStatus_ = true;
    curSoln_              = Teuchos::null;
  }




  void print(std::ostream& os, int indent = 0) const {
    os.setf(std::ios_base::scientific);
    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;
  }

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



  Teuchos::RCP<MV> getSolution() { return curSoln_; }

  int getQuorum() const { return quorum_; }

  size_t getSubIdx() const { return subIdx_; }

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



  StatusType firstCallCheckStatusSetup(Iteration<Scalar, MV, OP>* iSolver) {
    int i;
    MagnitudeType zero                      = Teuchos::ScalarTraits<MagnitudeType>::zero();
    MagnitudeType one                       = Teuchos::ScalarTraits<MagnitudeType>::one();
    const LinearProblem<Scalar, 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;

      // try to access the underlying blocked operator
      Teuchos::RCP<const OP> Op = lp.getOperator();
      Teuchos::RCP<const Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xOp =
          Teuchos::rcp_dynamic_cast<const Belos::XpetraOp<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(Op);
      TEUCHOS_TEST_FOR_EXCEPTION(xOp.is_null(), MueLu::Exceptions::BadCast, "Bad cast from \'const Belos::OperatorT\' to \'const Belos::XpetraOp\'. The origin type is " << typeid(const OP).name() << ".");
      Teuchos::RCP<const Xpetra::Operator<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xIntOp =
          xOp->getOperator();
      TEUCHOS_TEST_FOR_EXCEPTION(xIntOp.is_null(), MueLu::Exceptions::BadCast, "Cannot access Xpetra::Operator stored in Belos::XpetraOperator.");
      Teuchos::RCP<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > xMat =
          Teuchos::rcp_dynamic_cast<const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xIntOp);
      TEUCHOS_TEST_FOR_EXCEPTION(xMat.is_null(), MueLu::Exceptions::RuntimeError, "Cannot access Xpetra::Matrix stored in Belos::XpetraOp. Error.");
      Teuchos::RCP<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> > bMat = Teuchos::rcp_dynamic_cast<const Xpetra::BlockedCrsMatrix<Scalar, LocalOrdinal, GlobalOrdinal, Node> >(xMat);
      TEUCHOS_TEST_FOR_EXCEPTION(bMat.is_null(), MueLu::Exceptions::BadCast, "Bad cast from \'const Xpetra::Matrix\' to \'const Xpetra::BlockedCrsMatrix\'. The origin type is " << typeid(const Xpetra::Matrix<Scalar, LocalOrdinal, GlobalOrdinal, Node>).name() << ". Note: you need a BlockedCrsMatrix object for the StatusTestGenResSubNorm to work!");
      mapExtractor_ = bMat->getRangeMapExtractor();
      TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor_.is_null(), MueLu::Exceptions::RuntimeError, "Could not extract map extractor from BlockedCrsMatrix. Error.");
      TEUCHOS_TEST_FOR_EXCEPTION(mapExtractor_->NumMaps() <= subIdx_, MueLu::Exceptions::RuntimeError, "The multivector is only split into " << mapExtractor_->NumMaps() << " sub parts. Cannot access sub-block " << subIdx_ << ".");

      // calculate initial norms
      if (scaletype_ == NormOfRHS) {
        Teuchos::RCP<const MV> rhs = lp.getRHS();
        numrhs_                    = MVT::GetNumberVecs(*rhs);
        scalevector_.resize(numrhs_);
        MvSubNorm(*rhs, subIdx_, scalevector_, scalenormtype_);
      } else if (scaletype_ == NormOfInitRes) {
        Teuchos::RCP<const MV> init_res = lp.getInitResVec();
        numrhs_                         = MVT::GetNumberVecs(*init_res);
        scalevector_.resize(numrhs_);
        MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
      } else if (scaletype_ == NormOfPrecInitRes) {
        Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
        numrhs_                         = MVT::GetNumberVecs(*init_res);
        scalevector_.resize(numrhs_);
        MvSubNorm(*init_res, subIdx_, scalevector_, scalenormtype_);
      } else if (scaletype_ == NormOfFullInitRes) {
        Teuchos::RCP<const MV> init_res = lp.getInitResVec();
        numrhs_                         = MVT::GetNumberVecs(*init_res);
        scalevector_.resize(numrhs_);
        MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
        scalevalue_ = one;
      } else if (scaletype_ == NormOfFullPrecInitRes) {
        Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
        numrhs_                         = MVT::GetNumberVecs(*init_res);
        scalevector_.resize(numrhs_);
        MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
        scalevalue_ = one;
      } else if (scaletype_ == NormOfFullScaledInitRes) {
        Teuchos::RCP<const MV> init_res = lp.getInitResVec();
        numrhs_                         = MVT::GetNumberVecs(*init_res);
        scalevector_.resize(numrhs_);
        MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
        MvScalingRatio(*init_res, subIdx_, scalevalue_);
      } else if (scaletype_ == NormOfFullScaledPrecInitRes) {
        Teuchos::RCP<const MV> init_res = lp.getInitPrecResVec();
        numrhs_                         = MVT::GetNumberVecs(*init_res);
        scalevector_.resize(numrhs_);
        MVT::MvNorm(*init_res, scalevector_, scalenormtype_);
        MvScalingRatio(*init_res, subIdx_, scalevalue_);
      } 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;
  }


  std::string description() const {
    std::ostringstream oss;
    oss << "Belos::StatusTestGenResSubNorm<>: " << 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 << " Exp";
    oss << " Res Vec [" << subIdx_ << "]) ";

    // 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 [" << subIdx_ << "]";
        else if (scaletype_ == NormOfPrecInitRes)
          oss << " Prec Res0 [" << subIdx_ << "]";
        else if (scaletype_ == NormOfFullInitRes)
          oss << " Full Res0 [" << subIdx_ << "]";
        else if (scaletype_ == NormOfFullPrecInitRes)
          oss << " Full Prec Res0 [" << subIdx_ << "]";
        else if (scaletype_ == NormOfFullScaledInitRes)
          oss << " scaled Full Res0 [" << subIdx_ << "]";
        else if (scaletype_ == NormOfFullScaledPrecInitRes)
          oss << " scaled Full Prec Res0 [" << subIdx_ << "]";
        else
          oss << " RHS [" << subIdx_ << "]";
        oss << ")";
      }
    }

    // TODO add a tagging name

    return oss.str();
  }




  // calculate norm of partial multivector
  void MvSubNorm(const MV& mv, size_t block, std::vector<typename Teuchos::ScalarTraits<Scalar>::magnitudeType>& normVec, NormType type = TwoNorm) {
    Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);

    Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);
    MVT::MvNorm(*SubVec, normVec, type);
  }

  // calculate ration of sub-vector length to full vector length (for scalevalue_)
  void MvScalingRatio(const MV& mv, size_t block, MagnitudeType& lengthRatio) {
    Teuchos::RCP<const MV> input = Teuchos::rcpFromRef(mv);

    Teuchos::RCP<const MV> SubVec = mapExtractor_->ExtractVector(input, block);

    lengthRatio = Teuchos::as<MagnitudeType>(SubVec->getGlobalLength()) / Teuchos::as<MagnitudeType>(input->getGlobalLength());
  }



  MagnitudeType tolerance_;

  size_t subIdx_;

  int quorum_;

  bool showMaxResNormOnly_;

  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_;

  Teuchos::RCP<const ME> mapExtractor_;
};

}  // namespace Belos

#endif /* BELOS_XPETRA_STATUS_TEST_GEN_RES_SUB_NORM_HPP */