Thyra_SolveSupportTypes.hpp

// @HEADER
// *****************************************************************************
//    Thyra: Interfaces and Support for Abstract Numerical Algorithms
//
// Copyright 2004 NTESS and the Thyra contributors.
// SPDX-License-Identifier: BSD-3-Clause
// *****************************************************************************
// @HEADER
 
#ifndef THYRA_SOLVE_SUPPORT_TYPES_HPP
#define THYRA_SOLVE_SUPPORT_TYPES_HPP
 
#include "Thyra_OperatorVectorTypes.hpp"
#include "Teuchos_ParameterList.hpp"
#include "Teuchos_FancyOStream.hpp"
#include "Teuchos_Describable.hpp"
 
 
namespace Thyra {
 
 
enum ESolveMeasureNormType {
  SOLVE_MEASURE_ONE,
  SOLVE_MEASURE_NORM_RESIDUAL,
  SOLVE_MEASURE_NORM_SOLUTION,
  SOLVE_MEASURE_NORM_INIT_RESIDUAL,
  SOLVE_MEASURE_NORM_RHS
};
 
 
inline
const std::string toString(const ESolveMeasureNormType solveMeasureNormType)
{
  switch(solveMeasureNormType) {
    case SOLVE_MEASURE_ONE:
      return "SOLVE_MEASURE_ONE";
    case SOLVE_MEASURE_NORM_RESIDUAL:
      return "SOLVE_MEASURE_NORM_RESIDUAL";
    case SOLVE_MEASURE_NORM_SOLUTION:
      return "SOLVE_MEASURE_NORM_SOLUTION";
    case SOLVE_MEASURE_NORM_INIT_RESIDUAL:
      return "SOLVE_MEASURE_NORM_INIT_RESIDUAL";
    case SOLVE_MEASURE_NORM_RHS:
      return "SOLVE_MEASURE_NORM_RHS";
    default:
      TEUCHOS_TEST_FOR_EXCEPT(true);
  }
  TEUCHOS_UNREACHABLE_RETURN("");
}
 
 
struct SolveMeasureType {
  ESolveMeasureNormType  numerator;
  ESolveMeasureNormType  denominator;
  SolveMeasureType()
    :numerator(SOLVE_MEASURE_ONE),denominator(SOLVE_MEASURE_ONE)
    {}
  SolveMeasureType(ESolveMeasureNormType _numerator, ESolveMeasureNormType _denominator)
    :numerator(_numerator),denominator(_denominator)
    {}
  void set(ESolveMeasureNormType _numerator, ESolveMeasureNormType _denominator)
    { numerator = _numerator; denominator = _denominator; }
  bool useDefault() const
    { return ( numerator==SOLVE_MEASURE_ONE && denominator==SOLVE_MEASURE_ONE ); }
  bool operator()(ESolveMeasureNormType numerator_in,
    ESolveMeasureNormType denominator_in
    ) const
    { return ( numerator==numerator_in && denominator==denominator_in ); }
  bool contains(ESolveMeasureNormType measure) const
    { return ( numerator==measure || denominator==measure ); }
};
 
 
inline
std::ostream& operator<<(std::ostream &out, const SolveMeasureType &solveMeasureType)
{
  out << "("<<toString(solveMeasureType.numerator)
      << "/"<<toString(solveMeasureType.denominator)<<")";
  return out;
}
 
 
template<class Scalar>
class ReductionFunctional : public Teuchos::Describable {
public:
 
 
  typename ScalarTraits<Scalar>::magnitudeType
  reduce( const VectorBase<Scalar> &v ) const
    {
#ifdef THYRA_DEBUG
      TEUCHOS_TEST_FOR_EXCEPTION(!isCompatible(v), Exceptions::IncompatibleVectorSpaces,
        "Error, the vector v="<<v.description()<<" is not compatiable with"
        " *this="<<this->description()<<"!");
#endif
      return reduceImpl(v);
    }      
 
  bool isCompatible( const VectorBase<Scalar> &v ) const
    { return isCompatibleImpl(v); }
 
 
protected:
 
 
  virtual typename ScalarTraits<Scalar>::magnitudeType
  reduceImpl( const VectorBase<Scalar> &v ) const = 0;
 
  virtual bool isCompatibleImpl( const VectorBase<Scalar> &v ) const = 0;
 
 
};
 
 
template <class Scalar>
struct SolveCriteria {
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
  static ScalarMag unspecifiedTolerance() { return ScalarMag(-1.0); }
  SolveMeasureType solveMeasureType;
  ScalarMag requestedTol;
  RCP<ParameterList> extraParameters;
  RCP<const ReductionFunctional<Scalar> > numeratorReductionFunc;
  RCP<const ReductionFunctional<Scalar> > denominatorReductionFunc;
  SolveCriteria()
    : requestedTol(unspecifiedTolerance())
    {}
  SolveCriteria(
    SolveMeasureType solveMeasureType_in,
    ScalarMag requestedTol_in,
    const RCP<ParameterList> &extraParameters_in = Teuchos::null,
    const RCP<ReductionFunctional<Scalar> > &numeratorReductionFunc_in = Teuchos::null,
    const RCP<ReductionFunctional<Scalar> > &denominatorReductionFunc_in = Teuchos::null
    )
    : solveMeasureType(solveMeasureType_in),
      requestedTol(requestedTol_in), 
      extraParameters(extraParameters_in),
      numeratorReductionFunc(numeratorReductionFunc_in),
      denominatorReductionFunc(denominatorReductionFunc_in)
    {}
};
 
 
template<class Scalar>
std::ostream& operator<<(std::ostream &out, const SolveCriteria<Scalar> &solveCriteria)
{
  out << typeName(solveCriteria) << "{";
  out << "solveMeasureType="<<solveCriteria.solveMeasureType;
  out << ", requestedTol="<<solveCriteria.requestedTol;
  if (nonnull(solveCriteria.extraParameters)) {
    out << ", extraParameters="<<solveCriteria.extraParameters;
  }
  if (nonnull(solveCriteria.numeratorReductionFunc)) {
    out << ", numeratorReductionFunc="<<solveCriteria.numeratorReductionFunc->description();
  }
  if (nonnull(solveCriteria.denominatorReductionFunc)) {
    out << ", denominatorReductionFunc="<<solveCriteria.denominatorReductionFunc->description();
  }
  out << "}";
  return out;
}
 
 
class CatastrophicSolveFailure : public std::runtime_error
{public: CatastrophicSolveFailure(const std::string& what_arg) : std::runtime_error(what_arg) {}};
 
 
enum ESolveStatus {
  SOLVE_STATUS_CONVERGED        
  ,SOLVE_STATUS_UNCONVERGED     
  ,SOLVE_STATUS_UNKNOWN         
};
 
 
inline
const std::string toString(const ESolveStatus solveStatus)
{
  switch(solveStatus) {
    case SOLVE_STATUS_CONVERGED:    return "SOLVE_STATUS_CONVERGED";
    case SOLVE_STATUS_UNCONVERGED:  return "SOLVE_STATUS_UNCONVERGED";
    case SOLVE_STATUS_UNKNOWN:      return "SOLVE_STATUS_UNKNOWN";
    default: TEUCHOS_TEST_FOR_EXCEPT(true);
  }
  TEUCHOS_UNREACHABLE_RETURN("");
}
 
 
template <class Scalar>
struct SolveStatus {
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
  static ScalarMag unknownTolerance() { return ScalarMag(-1); }
  ESolveStatus solveStatus;
  ScalarMag achievedTol;
  std::string message;
  RCP<ParameterList> extraParameters;
  SolveStatus()
    :solveStatus(SOLVE_STATUS_UNKNOWN), achievedTol(unknownTolerance())
    {}
  static std::string achievedTolToString( const ScalarMag &achievedTol )
    {
      if(achievedTol==unknownTolerance()) return "unknownTolerance()";
      std::ostringstream oss; oss << achievedTol; return oss.str();
    }
};
 
 
template <class Scalar>
std::ostream& operator<<( std::ostream& out_arg, const SolveStatus<Scalar> &solveStatus )
{
  RCP<Teuchos::FancyOStream>
    out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
  Teuchos::OSTab tab(out);
  *out
    << "solveStatus = " << toString(solveStatus.solveStatus) << std::endl
    << "achievedTol = " << SolveStatus<Scalar>::achievedTolToString(solveStatus.achievedTol) << std::endl;
  *out << "message:";
  if (solveStatus.message.length()) {
    Teuchos::OSTab tab2(out);
    *out << "\n" << solveStatus.message << "\n";
  }
  *out << "extraParameters:";
  if(solveStatus.extraParameters.get()) {
    *out << "\n";
    Teuchos::OSTab tab3(out);
    solveStatus.extraParameters->print(*out, 10, true);
  }
  else {
    *out << " NONE\n";
  }
  return out_arg;
}
 
 
enum ESupportSolveUse {
  SUPPORT_SOLVE_UNSPECIFIED  
  ,SUPPORT_SOLVE_FORWARD_ONLY  
  ,SUPPORT_SOLVE_TRANSPOSE_ONLY  
  ,SUPPORT_SOLVE_FORWARD_AND_TRANSPOSE  
};
 
 
enum EPreconditionerInputType {
  PRECONDITIONER_INPUT_TYPE_AS_OPERATOR  
  ,PRECONDITIONER_INPUT_TYPE_AS_MATRIX   
};
 
 
template <class Scalar>
void accumulateSolveStatusInit(
  const Ptr<SolveStatus<Scalar> > &overallSolveStatus
  )
{
  overallSolveStatus->solveStatus = SOLVE_STATUS_CONVERGED;
}
 
 
template <class Scalar>
void accumulateSolveStatus(
  const SolveCriteria<Scalar>, // ToDo: Never used, need to take this out!
  const SolveStatus<Scalar> &solveStatus,
  const Ptr<SolveStatus<Scalar> > &overallSolveStatus
  )
{
  switch(solveStatus.solveStatus) {
    case SOLVE_STATUS_UNCONVERGED:
    {
      // First, if we see any unconverged solve status, then the entire block is
      // unconverged!
      overallSolveStatus->solveStatus = SOLVE_STATUS_UNCONVERGED;
      overallSolveStatus->message = solveStatus.message;
      overallSolveStatus->extraParameters = solveStatus.extraParameters;
      break;
    }
    case SOLVE_STATUS_UNKNOWN:
    {
      // Next, if any solve status is unknown, then if the overall solve
      // status says converged, then we have to mark it as unknown.  Note that
      // unknown could mean that the system is actually converged!
      switch(overallSolveStatus->solveStatus) {
        case SOLVE_STATUS_CONVERGED:
          overallSolveStatus->solveStatus = SOLVE_STATUS_UNKNOWN;
          break;
        case SOLVE_STATUS_UNCONVERGED:
        case SOLVE_STATUS_UNKNOWN:
          // If we get here then the overall solve status is either unknown
          // already or says unconverged and this will not change here!
          overallSolveStatus->message = solveStatus.message;
          overallSolveStatus->extraParameters = solveStatus.extraParameters;
          break;
        default:
          TEUCHOS_TEST_FOR_EXCEPT(true); // Corrupted enum?
      }
      break;
    }
    case SOLVE_STATUS_CONVERGED:
    {
      // If we get here then the overall solve status is either unknown,
      // unconverged, or converged and this will not change here!
      if(overallSolveStatus->message == "")
        overallSolveStatus->message = solveStatus.message;
      break;
    }
    default:
      TEUCHOS_TEST_FOR_EXCEPT(true); // Corrupted enum?
  }
  // Update the achieved tolerence to the maximum returned
  if( solveStatus.achievedTol > overallSolveStatus->achievedTol ) {
    overallSolveStatus->achievedTol = solveStatus.achievedTol;
  }
  // Set a message if none is set
  if(overallSolveStatus->message == "")
    overallSolveStatus->message = solveStatus.message;
  // Set the extra parameters if none is set
  if(overallSolveStatus->extraParameters.get()==NULL)
    overallSolveStatus->extraParameters = solveStatus.extraParameters;
}
 
 
} // namespace Thyra
 
 
#endif // THYRA_SOLVE_SUPPORT_TYPES_HPP