Thyra Version of the Day
Thyra_SolveSupportTypes.hpp
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_SOLVE_SUPPORT_TYPES_HPP
00030 #define THYRA_SOLVE_SUPPORT_TYPES_HPP
00031 
00032 #include "Thyra_OperatorVectorTypes.hpp"
00033 #include "Teuchos_ParameterList.hpp"
00034 #include "Teuchos_FancyOStream.hpp"
00035 #include "Teuchos_Describable.hpp"
00036 
00037 
00038 namespace Thyra {
00039 
00040 
00047 enum ESolveMeasureNormType {
00049   SOLVE_MEASURE_ONE,
00051   SOLVE_MEASURE_NORM_RESIDUAL,
00053   SOLVE_MEASURE_NORM_SOLUTION,
00055   SOLVE_MEASURE_NORM_INIT_RESIDUAL,
00057   SOLVE_MEASURE_NORM_RHS
00058 };
00059 
00060 
00065 inline
00066 const std::string toString(const ESolveMeasureNormType solveMeasureNormType)
00067 {
00068   switch(solveMeasureNormType) {
00069     case SOLVE_MEASURE_ONE:
00070       return "SOLVE_MEASURE_ONE";
00071     case SOLVE_MEASURE_NORM_RESIDUAL:
00072       return "SOLVE_MEASURE_NORM_RESIDUAL";
00073     case SOLVE_MEASURE_NORM_SOLUTION:
00074       return "SOLVE_MEASURE_NORM_SOLUTION";
00075     case SOLVE_MEASURE_NORM_INIT_RESIDUAL:
00076       return "SOLVE_MEASURE_NORM_INIT_RESIDUAL";
00077     case SOLVE_MEASURE_NORM_RHS:
00078       return "SOLVE_MEASURE_NORM_RHS";
00079     default:
00080       TEST_FOR_EXCEPT(true);
00081   }
00082   return NULL; // Never be called!
00083 }
00084 
00085 
00099 struct SolveMeasureType {
00101   ESolveMeasureNormType  numerator;
00103   ESolveMeasureNormType  denominator;
00105   SolveMeasureType()
00106     :numerator(SOLVE_MEASURE_ONE),denominator(SOLVE_MEASURE_ONE)
00107     {}
00109   SolveMeasureType(ESolveMeasureNormType _numerator, ESolveMeasureNormType _denominator)
00110     :numerator(_numerator),denominator(_denominator)
00111     {}
00113   void set(ESolveMeasureNormType _numerator, ESolveMeasureNormType _denominator)
00114     { numerator = _numerator; denominator = _denominator; }
00118   bool useDefault() const
00119     { return ( numerator==SOLVE_MEASURE_ONE && denominator==SOLVE_MEASURE_ONE ); }
00121   bool operator()(ESolveMeasureNormType numerator_in,
00122     ESolveMeasureNormType denominator_in
00123     ) const
00124     { return ( numerator==numerator_in && denominator==denominator_in ); }
00126   bool contains(ESolveMeasureNormType measure) const
00127     { return ( numerator==measure || denominator==measure ); }
00128 };
00129 
00130 
00135 inline
00136 std::ostream& operator<<(std::ostream &out, const SolveMeasureType &solveMeasureType)
00137 {
00138   out << "("<<toString(solveMeasureType.numerator)
00139       << "/"<<toString(solveMeasureType.denominator)<<")";
00140   return out;
00141 }
00142 
00143 
00147 template<class Scalar>
00148 class ReductionFunctional : public Teuchos::Describable {
00149 public:
00150 
00153 
00162   typename ScalarTraits<Scalar>::magnitudeType
00163   reduce( const VectorBase<Scalar> &v ) const
00164     {
00165 #ifdef THYRA_DEBUG
00166       TEST_FOR_EXCEPTION(!isCompatible(v), Exceptions::IncompatibleVectorSpaces,
00167         "Error, the vector v="<<v.description()<<" is not compatiable with"
00168         " *this="<<this->description()<<"!");
00169 #endif
00170       return reduceImpl(v);
00171     }      
00172 
00176   bool isCompatible( const VectorBase<Scalar> &v ) const
00177     { return isCompatibleImpl(v); }
00178 
00180 
00181 protected:
00182 
00185 
00187   virtual typename ScalarTraits<Scalar>::magnitudeType
00188   reduceImpl( const VectorBase<Scalar> &v ) const = 0;
00189 
00191   virtual bool isCompatibleImpl( const VectorBase<Scalar> &v ) const = 0;
00192 
00194 
00195 };
00196 
00197 
00298 template <class Scalar>
00299 struct SolveCriteria {
00301   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00303   static ScalarMag unspecifiedTolerance() { return ScalarMag(-1.0); }
00306   SolveMeasureType solveMeasureType;
00309   ScalarMag requestedTol;
00315   RCP<ParameterList> extraParameters;
00318   RCP<const ReductionFunctional<Scalar> > numeratorReductionFunc;
00321   RCP<const ReductionFunctional<Scalar> > denominatorReductionFunc;
00323   SolveCriteria()
00324     : requestedTol(unspecifiedTolerance())
00325     {}
00327   SolveCriteria(
00328     SolveMeasureType solveMeasureType_in,
00329     ScalarMag requestedTol_in,
00330     const RCP<ParameterList> &extraParameters_in = Teuchos::null,
00331     const RCP<ReductionFunctional<Scalar> > &numeratorReductionFunc_in = Teuchos::null,
00332     const RCP<ReductionFunctional<Scalar> > &denominatorReductionFunc_in = Teuchos::null
00333     )
00334     : solveMeasureType(solveMeasureType_in),
00335       requestedTol(requestedTol_in), 
00336       extraParameters(extraParameters_in),
00337       numeratorReductionFunc(numeratorReductionFunc_in),
00338       denominatorReductionFunc(denominatorReductionFunc_in)
00339     {}
00340 };
00341 
00342 
00347 template<class Scalar>
00348 std::ostream& operator<<(std::ostream &out, const SolveCriteria<Scalar> &solveCriteria)
00349 {
00350   out << typeName(solveCriteria) << "{";
00351   out << "solveMeasureType="<<solveCriteria.solveMeasureType;
00352   out << ", requestedTol="<<solveCriteria.requestedTol;
00353   if (nonnull(solveCriteria.extraParameters)) {
00354     out << ", extraParameters="<<solveCriteria.extraParameters;
00355   }
00356   if (nonnull(solveCriteria.numeratorReductionFunc)) {
00357     out << ", numeratorReductionFunc="<<solveCriteria.numeratorReductionFunc->description();
00358   }
00359   if (nonnull(solveCriteria.denominatorReductionFunc)) {
00360     out << ", denominatorReductionFunc="<<solveCriteria.denominatorReductionFunc->description();
00361   }
00362   out << "}";
00363   return out;
00364 }
00365 
00366 
00371 template <class Scalar>
00372 struct THYRA_DEPRECATED BlockSolveCriteria {
00374   SolveCriteria<Scalar> solveCriteria;
00376   int                     numRhs;
00378   BlockSolveCriteria()
00379     : solveCriteria(), numRhs(1)
00380     {}
00382   BlockSolveCriteria( const SolveCriteria<Scalar> &_solveCriteria, int _numRhs )
00383     : solveCriteria(_solveCriteria), numRhs(_numRhs)
00384     {}
00385 };
00386 
00387 
00392 class CatastrophicSolveFailure : public std::runtime_error
00393 {public: CatastrophicSolveFailure(const std::string& what_arg) : std::runtime_error(what_arg) {}};
00394 
00395 
00400 enum ESolveStatus {
00401   SOLVE_STATUS_CONVERGED        
00402   ,SOLVE_STATUS_UNCONVERGED     
00403   ,SOLVE_STATUS_UNKNOWN         
00404 };
00405 
00406 
00411 inline
00412 const std::string toString(const ESolveStatus solveStatus)
00413 {
00414   switch(solveStatus) {
00415     case SOLVE_STATUS_CONVERGED:    return "SOLVE_STATUS_CONVERGED";
00416     case SOLVE_STATUS_UNCONVERGED:  return "SOLVE_STATUS_UNCONVERGED";
00417     case SOLVE_STATUS_UNKNOWN:      return "SOLVE_STATUS_UNKNOWN";
00418     default: TEST_FOR_EXCEPT(true);
00419   }
00420   return ""; // Never be called!
00421 }
00422 
00423 
00430 template <class Scalar>
00431 struct SolveStatus {
00433   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00435   static ScalarMag unknownTolerance() { return ScalarMag(-1); }
00437   ESolveStatus solveStatus;
00441   ScalarMag achievedTol;
00443   std::string message;
00446   RCP<ParameterList> extraParameters;
00448   SolveStatus()
00449     :solveStatus(SOLVE_STATUS_UNKNOWN), achievedTol(unknownTolerance())
00450     {}
00453   static std::string achievedTolToString( const ScalarMag &achievedTol )
00454     {
00455       if(achievedTol==unknownTolerance()) return "unknownTolerance()";
00456       std::ostringstream oss; oss << achievedTol; return oss.str();
00457     }
00458 };
00459 
00460 
00465 template <class Scalar>
00466 std::ostream& operator<<( std::ostream& out_arg, const SolveStatus<Scalar> &solveStatus )
00467 {
00468   RCP<Teuchos::FancyOStream>
00469     out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
00470   Teuchos::OSTab tab(out);
00471   *out
00472     << "solveStatus = " << toString(solveStatus.solveStatus) << std::endl
00473     << "achievedTol = " << SolveStatus<Scalar>::achievedTolToString(solveStatus.achievedTol) << std::endl;
00474   *out << "message:";
00475   if (solveStatus.message.length()) {
00476     Teuchos::OSTab tab2(out);
00477     *out << "\n" << solveStatus.message << "\n";
00478   }
00479   *out << "extraParameters:";
00480   if(solveStatus.extraParameters.get()) {
00481     *out << "\n";
00482     Teuchos::OSTab tab3(out);
00483     solveStatus.extraParameters->print(*out, 1000, true);
00484   }
00485   else {
00486     *out << " NONE\n";
00487   }
00488   return out_arg;
00489 }
00490 
00491 
00497 enum ESupportSolveUse {
00498   SUPPORT_SOLVE_UNSPECIFIED  
00499   ,SUPPORT_SOLVE_FORWARD_ONLY  
00500   ,SUPPORT_SOLVE_TRANSPOSE_ONLY  
00501   ,SUPPORT_SOLVE_FORWARD_AND_TRANSPOSE  
00502 };
00503 
00504 
00509 enum EPreconditionerInputType {
00510   PRECONDITIONER_INPUT_TYPE_AS_OPERATOR  
00511   ,PRECONDITIONER_INPUT_TYPE_AS_MATRIX   
00512 };
00513 
00514 
00519 template <class Scalar>
00520 void accumulateSolveStatusInit(
00521   const Ptr<SolveStatus<Scalar> > &overallSolveStatus
00522   )
00523 {
00524   overallSolveStatus->solveStatus = SOLVE_STATUS_CONVERGED;
00525 }
00526 
00527 
00544 template <class Scalar>
00545 void accumulateSolveStatus(
00546   const SolveCriteria<Scalar>, // ToDo: Never used, need to take this out!
00547   const SolveStatus<Scalar> &solveStatus,
00548   const Ptr<SolveStatus<Scalar> > &overallSolveStatus
00549   )
00550 {
00551   switch(solveStatus.solveStatus) {
00552     case SOLVE_STATUS_UNCONVERGED:
00553     {
00554       // First, if we see any unconverged solve status, then the entire block is
00555       // unconverged!
00556       overallSolveStatus->solveStatus = SOLVE_STATUS_UNCONVERGED;
00557       overallSolveStatus->message = solveStatus.message;
00558       overallSolveStatus->extraParameters = solveStatus.extraParameters;
00559       break;
00560     }
00561     case SOLVE_STATUS_UNKNOWN:
00562     {
00563       // Next, if any solve status is unknown, then if the overall solve
00564       // status says converged, then we have to mark it as unknown.  Note that
00565       // unknown could mean that the system is actually converged!
00566       switch(overallSolveStatus->solveStatus) {
00567         case SOLVE_STATUS_CONVERGED:
00568           overallSolveStatus->solveStatus = SOLVE_STATUS_UNKNOWN;
00569           break;
00570         case SOLVE_STATUS_UNCONVERGED:
00571         case SOLVE_STATUS_UNKNOWN:
00572           // If we get here then the overall solve status is either unknown
00573           // already or says unconverged and this will not change here!
00574           overallSolveStatus->message = solveStatus.message;
00575           overallSolveStatus->extraParameters = solveStatus.extraParameters;
00576           break;
00577         default:
00578           TEST_FOR_EXCEPT(true); // Corrupted enum?
00579       }
00580       break;
00581     }
00582     case SOLVE_STATUS_CONVERGED:
00583     {
00584       // If we get here then the overall solve status is either unknown,
00585       // unconverged, or converged and this will not change here!
00586       if(overallSolveStatus->message == "")
00587         overallSolveStatus->message = solveStatus.message;
00588       break;
00589     }
00590     default:
00591       TEST_FOR_EXCEPT(true); // Corrupted enum?
00592   }
00593   // Update the achieved tolerence to the maximum returned
00594   if( solveStatus.achievedTol > overallSolveStatus->achievedTol ) {
00595     overallSolveStatus->achievedTol = solveStatus.achievedTol;
00596   }
00597   // Set a message if none is set
00598   if(overallSolveStatus->message == "")
00599     overallSolveStatus->message = solveStatus.message;
00600   // Set the extra parameters if none is set
00601   if(overallSolveStatus->extraParameters.get()==NULL)
00602     overallSolveStatus->extraParameters = solveStatus.extraParameters;
00603 }
00604 
00605 
00610 template <class Scalar>
00611 THYRA_DEPRECATED
00612 void accumulateSolveStatus(
00613   const SolveCriteria<Scalar>, // ToDo: Never used, need to take this out!
00614   const SolveStatus<Scalar> &solveStatus,
00615   SolveStatus<Scalar> *overallSolveStatus
00616   )
00617 {
00618   accumulateSolveStatus(
00619     SolveCriteria<Scalar>(),
00620     solveStatus, Teuchos::ptr(overallSolveStatus)
00621     );
00622 }
00623 
00624 
00625 } // namespace Thyra
00626 
00627 
00628 #endif // THYRA_SOLVE_SUPPORT_TYPES_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines