Thyra_SolveSupportTypes.hpp

// @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"

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, ESolveMeasureNormType _denominator) const
    { return ( numerator==_numerator && denominator==_denominator ); }
};

inline
std::string toString(const SolveMeasureType& solveMeasureType)
{
  std::ostringstream oss;
  oss << "("<<toString(solveMeasureType.numerator)<<")/("<<solveMeasureType.denominator<<")";
  return oss.str();
}

template <class Scalar>
struct SolveCriteria {
  typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
  static ScalarMag unspecifiedTolerance() { return ScalarMag(-1); }
  SolveMeasureType solveMeasureType;
  ScalarMag requestedTol;
  Teuchos::RCP<Teuchos::ParameterList> extraParameters;
  SolveCriteria()
    :solveMeasureType()
    ,requestedTol(unspecifiedTolerance())
    {}
  SolveCriteria(
    SolveMeasureType _solveMeasureType, ScalarMag _requestedTol
    ,const Teuchos::RCP<Teuchos::ParameterList> &_extraParameters = Teuchos::null
    )
    :solveMeasureType(_solveMeasureType),requestedTol(_requestedTol),extraParameters(_extraParameters)
    {}
};

template <class Scalar>
struct 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;
  Teuchos::RCP<Teuchos::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 )
{
  Teuchos::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 tab(out);
    *out << "\n" << solveStatus.message << "\n";
  }
  *out << "extraParameters:";
  if(solveStatus.extraParameters.get()) {
    *out << "\n";
    solveStatus.extraParameters->print(Teuchos::OSTab(out).o(),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 accumulateSolveStatus(
  const SolveCriteria<Scalar>    &overallSolveCriteria
  ,const SolveStatus<Scalar>     &solveStatus
  ,SolveStatus<Scalar>           *overallSolveStatus
  )
{
#ifdef TEUCHOS_DEBUG
  TEST_FOR_EXCEPT(overallSolveStatus==NULL);
#endif
  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;
}

} // namespace Thyra

#endif // THYRA_SOLVE_SUPPORT_TYPES_HPP

---