BelosLSQRStatusTest.hpp

Go to the documentation of this file.

/*
// @HEADER
// *****************************************************************************
//                 Belos: Block Linear Solvers Package
//
// Copyright 2004-2016 NTESS and the Belos contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
*/

#ifndef BELOS_LSQR_STATUS_TEST_HPP
#define BELOS_LSQR_STATUS_TEST_HPP

#include "BelosStatusTest.hpp"
#include "BelosLSQRIter.hpp"

namespace Belos {


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

public:

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




  LSQRStatusTest( MagnitudeType condMax = 0.0,
                  int term_iter_max = 1,
                  MagnitudeType rel_rhs_err = 0.0,
                  MagnitudeType rel_mat_err = 0.0 );

  virtual ~LSQRStatusTest();




  Belos::StatusType checkStatus(Belos::Iteration<ScalarType,MV,OP> *iSolver );

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




  void reset();

  int setCondLim(MagnitudeType condMax) {
    condMax_ = condMax;
    rcondMin_ = (condMax > 0) ? (Teuchos::ScalarTraits< MagnitudeType >::one() / condMax) : Teuchos::ScalarTraits< MagnitudeType >::eps();
    return(0);}

  int setTermIterMax(int term_iter_max) {
    term_iter_max_ = term_iter_max;
    if (term_iter_max_ < 1)
      term_iter_max_ = 1;
    return(0);}

  int setRelRhsErr(MagnitudeType rel_rhs_err) {
    rel_rhs_err_ = rel_rhs_err;
    return(0);}

  int setRelMatErr(MagnitudeType rel_mat_err) {
    rel_mat_err_ = rel_mat_err;
    return(0);}




  MagnitudeType getCondMaxLim() const {return(condMax_);}

  int getTermIterMax() const {return(term_iter_max_);}

  MagnitudeType getRelRhsErr() const {return(rel_rhs_err_);}

  MagnitudeType getMatErr() const {return(rel_mat_err_);}

  MagnitudeType getMatCondNum() const {return(matCondNum_);}

  MagnitudeType getMatNorm() const {return(matNorm_);}

  int getTermIter() const { return term_iter_; }

  MagnitudeType getResidNorm() const {return(resNorm_);}

  MagnitudeType getLSResidNorm() const {return(matResNorm_);}




  void print(std::ostream& os, int indent = 0) const;

  void printStatus(std::ostream& os, Belos::StatusType type) const;



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


  std::string description() const
  {
    std::ostringstream oss;
    oss << "LSQRStatusTest<>: [ limit of condition number = " << condMax_ << " ]";
    return oss.str();
  }

private:



  MagnitudeType condMax_;

  int term_iter_max_;

  MagnitudeType rel_rhs_err_;

  MagnitudeType rel_mat_err_;

  MagnitudeType rcondMin_;

  Belos::StatusType status_;

  // term_iter_ records the number of consecutive "successful" iterations.
  // convergence requires that term_iter_max consecutive iterates satisfy the other convergence tests
  int term_iter_;

  // condition number of the operator
  MagnitudeType matCondNum_;

  // Frobenius norm of the operator
  MagnitudeType matNorm_;

  // residual norm for the linear system
  MagnitudeType resNorm_;

  // least squares residual, operator^Transpose * residual
  MagnitudeType matResNorm_;


};

template <class ScalarType, class MV, class OP>
LSQRStatusTest<ScalarType,MV,OP>::
LSQRStatusTest (MagnitudeType condMax /* = 0 */,
                int term_iter_max /* = 1 */,
                MagnitudeType rel_rhs_err /* = 0 */,
                MagnitudeType rel_mat_err /* = 0 */)
  : condMax_(condMax),
    term_iter_max_ (term_iter_max),
    rel_rhs_err_ (rel_rhs_err),
    rel_mat_err_ (rel_mat_err),
    rcondMin_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
    status_ (Belos::Undefined),
    term_iter_ (0),
    matCondNum_ ( Teuchos::ScalarTraits<MagnitudeType>::one() ),
    matNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
    resNorm_  ( Teuchos::ScalarTraits<MagnitudeType>::zero() ),
    matResNorm_ ( Teuchos::ScalarTraits<MagnitudeType>::zero() )
{}

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

template <class ScalarType, class MV, class OP>
void LSQRStatusTest<ScalarType,MV,OP>::reset()
{
  status_ = Belos::Undefined;
}

template <class ScalarType, class MV, class OP>
Belos::StatusType LSQRStatusTest<ScalarType,MV,OP>::checkStatus( Belos::Iteration<ScalarType,MV,OP>* iSolver)
{
  const MagnitudeType MTzero = Teuchos::ScalarTraits<MagnitudeType>::zero();
  const MagnitudeType MTone = Teuchos::ScalarTraits<MagnitudeType>::one();
  if (condMax_ > MTzero )
    {
        rcondMin_ = MTone / condMax_;
    }
  else
    {
        rcondMin_ = Teuchos::ScalarTraits< MagnitudeType >::eps();
    }

  bool termIterFlag = false;
  LSQRIter<ScalarType,MV,OP>* solver = dynamic_cast< LSQRIter<ScalarType,MV,OP>* > (iSolver);
  TEUCHOS_ASSERT(solver != NULL);
  LSQRIterationState< ScalarType, MV > state = solver->getState();
  //
  //   LSQR solves a least squares problem.  A converged preconditioned residual norm
  // suffices for convergence, but is not necessary.  LSQR sometimes returns a larger
  // relative residual norm than what would have been returned by a linear solver.
  // This section evaluates three stopping criteria.  In the Solver Manager, this test
  // is combined with a generic number of iteration test.
  //   If the linear system includes a preconditioner, then the least squares problem
  // is solved for the preconditioned linear system.  Preconditioning changes the least
  // squares problem (in the sense of changing the norms), and the solution depends
  // on the preconditioner in this sense.
  //   In the context of Linear Least Squares problems, preconditioning refers
  // to the regularization matrix.  Here the regularization matrix is always a scalar
  // multiple of the identity (standard form least squres).
  //   The "loss of accuracy" concept is not yet implemented here, becuase it is unclear
  // what this means for linear least squares.  LSQR solves an inconsistent system
  // in a least-squares sense.  "Loss of accuracy" would correspond to
  // the difference between the preconditioned residual and the unpreconditioned residual.
  //

  std::cout << " X " << state.sol_norm
            << "  b-AX " << state.resid_norm
            << "  Atr  " << state.mat_resid_norm
            << "  A " << state.frob_mat_norm
            << "  cond  " << state.mat_cond_num
            << "  relResNorm " << state.resid_norm/state.bnorm
            << "  LS " << state.mat_resid_norm /( state.resid_norm * state.frob_mat_norm )
            << std::endl;

  const MagnitudeType zero = Teuchos::ScalarTraits<MagnitudeType>::zero();
  const ScalarType one = Teuchos::ScalarTraits<ScalarType>::one();
  ScalarType stop_crit_1 = zero; // b = 0, done
  if( state.bnorm  > zero )
    {
      stop_crit_1 = state.resid_norm / state.bnorm;
    }
  ScalarType stop_crit_2 = zero;
  if( state.frob_mat_norm  > zero  && state.resid_norm > zero )
    {
      stop_crit_2 = (state.resid_norm > zero) ? state.mat_resid_norm / (state.frob_mat_norm * state.resid_norm) : zero;
    }
  else
    {
     if( state.resid_norm == zero )
       {
         stop_crit_2 = zero;
       }
     else
       {
         stop_crit_2 = one; // Initial mat_norm always vanishes
       }
    }
  ScalarType stop_crit_3 = one / state.mat_cond_num;
  ScalarType resid_tol = rel_rhs_err_ + rel_mat_err_ * state.frob_mat_norm * state.sol_norm / state.bnorm;
  ScalarType resid_tol_mach = Teuchos::ScalarTraits< MagnitudeType >::eps() + Teuchos::ScalarTraits< MagnitudeType >::eps() * state.frob_mat_norm * state.sol_norm / state.bnorm;

  // The expected use case for our users is that the linear system will almost
  // always be compatible, but occasionally may not be.  However, some users
  // may use LSQR for more general cases.  This is why we include the full
  // suite of tests, for both compatible and incompatible systems.
  //
  // Users will have to be educated that sometimes they will get an answer X
  // that does _not_ satisfy the linear system AX=B, but _does_ satisfy the
  // corresponding least-squares problem.  Perhaps the solution manager should
  // provide them with a way to find out.

  // stop_crit_1 is for compatible linear systems.
  // stop_crit_2 is for incompatible linear systems.
  // stop_crit_3 is for either compatible or incompatible linear systems.

  // Have we met any of the stopping criteria?
  if (stop_crit_1 <= resid_tol || stop_crit_2 <= rel_mat_err_ || stop_crit_3 <= rcondMin_ || stop_crit_1 <= resid_tol_mach || stop_crit_2 <= Teuchos::ScalarTraits< MagnitudeType >::eps() || stop_crit_3 <= Teuchos::ScalarTraits< MagnitudeType >::eps()) {
    termIterFlag = true;

    if (stop_crit_1 <= resid_tol )
      std::cout << "Conv: stop_crit_1  " << stop_crit_1  << " resid_tol " << resid_tol << std::endl;

    if (stop_crit_1 <=  resid_tol_mach )
      std::cout << "Conv: stop_crit_1  " << stop_crit_1  << " resid_tol_mach " << resid_tol_mach << std::endl;

    if (stop_crit_2 <= rel_mat_err_ )
      std::cout << "Conv: stop_crit_2  " << stop_crit_2  << " rel_mat_err " << rel_mat_err_ << std::endl;

    if (stop_crit_2 <=   Teuchos::ScalarTraits< MagnitudeType >::eps() )
      std::cout << "Conv: stop_crit_2  " << stop_crit_2  << " eps " <<   Teuchos::ScalarTraits< MagnitudeType >::eps()   << std::endl;

    if (stop_crit_3 <= rcondMin_ )
      std::cout << "Conv: stop_crit_3  " << stop_crit_3  << " rcondMin_ " << rcondMin_ << std::endl;

    if (stop_crit_3 <=   Teuchos::ScalarTraits< MagnitudeType >::eps() )
      std::cout << "Conv: stop_crit_3  " << stop_crit_3  << " eps " <<   Teuchos::ScalarTraits< MagnitudeType >::eps()   << std::endl;
  }

  // update number of consecutive successful iterations
  if (!termIterFlag) {
    term_iter_ = 0;
  } else {
    term_iter_++;
  }
  status_ = (term_iter_ < term_iter_max_) ? Belos::Failed : Belos::Passed;

  matCondNum_ = state.mat_cond_num; // information that defined convergence
  matNorm_ = state.frob_mat_norm;   // in accessible variables
  resNorm_  = state.resid_norm;
  matResNorm_ = state.mat_resid_norm;

  return status_;
}

template <class ScalarType, class MV, class OP>
void LSQRStatusTest<ScalarType,MV,OP>::print(std::ostream& os, int indent) const
{
  for (int j = 0; j < indent; j++)
    os << ' ';
  printStatus(os, status_);
  os << "limit of condition number = " << condMax_ << std::endl;
  os << "limit of condition number = " << condMax_ << std::endl;
}

template <class ScalarType, class MV, class OP>
void LSQRStatusTest<ScalarType,MV,OP>::printStatus(std::ostream&os, Belos::StatusType type) const
{
  os << std::left << std::setw(13) << std::setfill('.');
  switch (type) {
  case Belos::Passed:
    os << "Passed";
    break;
  case Belos::Failed:
    os << "Failed";
    break;
  case Belos::Undefined:
  default:
    os << "Undefined";
    break;
  }
  os << std::left << std::setfill(' ');
  return;
}

} // end Belos namespace


#endif /* BELOS_LSQR_STATUS_TEST_HPP */