Thyra_SolveSupportTypes.hpp

Go to the documentation of this file.

// @HEADER
// ***********************************************************************
// 
//    Thyra: Interfaces and Support for Abstract Numerical Algorithms
//                 Copyright (2004) Sandia Corporation
// 
// Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
// license for use of this work by or on behalf of the U.S. Government.
// 
// This library is free software; you can redistribute it and/or modify
// it under the terms of the GNU Lesser General Public License as
// published by the Free Software Foundation; either version 2.1 of the
// License, or (at your option) any later version.
//  
// This library is distributed in the hope that it will be useful, but
// WITHOUT ANY WARRANTY; without even the implied warranty of
// MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
// Lesser General Public License for more details.
//  
// You should have received a copy of the GNU Lesser General Public
// License along with this library; if not, write to the Free Software
// Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
// USA
// Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
// 
// ***********************************************************************
// @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:
      TEST_FOR_EXCEPT(true);
  }
  return NULL; // Never be called!
}


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


template <class Scalar>
struct THYRA_DEPRECATED BlockSolveCriteria {
  SolveCriteria<Scalar> solveCriteria;
  int                     numRhs;
  BlockSolveCriteria()
    : solveCriteria(), numRhs(1)
    {}
  BlockSolveCriteria( const SolveCriteria<Scalar> &_solveCriteria, int _numRhs )
    : solveCriteria(_solveCriteria), numRhs(_numRhs)
    {}
};


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: TEST_FOR_EXCEPT(true);
  }
  return ""; // Never be called!
}


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, 1000, 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:
          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:
      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;
}


template <class Scalar>
THYRA_DEPRECATED
void accumulateSolveStatus(
  const SolveCriteria<Scalar>, // ToDo: Never used, need to take this out!
  const SolveStatus<Scalar> &solveStatus,
  SolveStatus<Scalar> *overallSolveStatus
  )
{
  accumulateSolveStatus(
    SolveCriteria<Scalar>(),
    solveStatus, Teuchos::ptr(overallSolveStatus)
    );
}


} // namespace Thyra


#endif // THYRA_SOLVE_SUPPORT_TYPES_HPP