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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef THYRA_SOLVE_SUPPORT_TYPES_HPP
00043 #define THYRA_SOLVE_SUPPORT_TYPES_HPP
00044 
00045 #include "Thyra_OperatorVectorTypes.hpp"
00046 #include "Teuchos_ParameterList.hpp"
00047 #include "Teuchos_FancyOStream.hpp"
00048 #include "Teuchos_Describable.hpp"
00049 
00050 
00051 namespace Thyra {
00052 
00053 
00060 enum ESolveMeasureNormType {
00062   SOLVE_MEASURE_ONE,
00064   SOLVE_MEASURE_NORM_RESIDUAL,
00066   SOLVE_MEASURE_NORM_SOLUTION,
00068   SOLVE_MEASURE_NORM_INIT_RESIDUAL,
00070   SOLVE_MEASURE_NORM_RHS
00071 };
00072 
00073 
00078 inline
00079 const std::string toString(const ESolveMeasureNormType solveMeasureNormType)
00080 {
00081   switch(solveMeasureNormType) {
00082     case SOLVE_MEASURE_ONE:
00083       return "SOLVE_MEASURE_ONE";
00084     case SOLVE_MEASURE_NORM_RESIDUAL:
00085       return "SOLVE_MEASURE_NORM_RESIDUAL";
00086     case SOLVE_MEASURE_NORM_SOLUTION:
00087       return "SOLVE_MEASURE_NORM_SOLUTION";
00088     case SOLVE_MEASURE_NORM_INIT_RESIDUAL:
00089       return "SOLVE_MEASURE_NORM_INIT_RESIDUAL";
00090     case SOLVE_MEASURE_NORM_RHS:
00091       return "SOLVE_MEASURE_NORM_RHS";
00092     default:
00093       TEUCHOS_TEST_FOR_EXCEPT(true);
00094   }
00095   return NULL; // Never be called!
00096 }
00097 
00098 
00112 struct SolveMeasureType {
00114   ESolveMeasureNormType  numerator;
00116   ESolveMeasureNormType  denominator;
00118   SolveMeasureType()
00119     :numerator(SOLVE_MEASURE_ONE),denominator(SOLVE_MEASURE_ONE)
00120     {}
00122   SolveMeasureType(ESolveMeasureNormType _numerator, ESolveMeasureNormType _denominator)
00123     :numerator(_numerator),denominator(_denominator)
00124     {}
00126   void set(ESolveMeasureNormType _numerator, ESolveMeasureNormType _denominator)
00127     { numerator = _numerator; denominator = _denominator; }
00131   bool useDefault() const
00132     { return ( numerator==SOLVE_MEASURE_ONE && denominator==SOLVE_MEASURE_ONE ); }
00134   bool operator()(ESolveMeasureNormType numerator_in,
00135     ESolveMeasureNormType denominator_in
00136     ) const
00137     { return ( numerator==numerator_in && denominator==denominator_in ); }
00139   bool contains(ESolveMeasureNormType measure) const
00140     { return ( numerator==measure || denominator==measure ); }
00141 };
00142 
00143 
00148 inline
00149 std::ostream& operator<<(std::ostream &out, const SolveMeasureType &solveMeasureType)
00150 {
00151   out << "("<<toString(solveMeasureType.numerator)
00152       << "/"<<toString(solveMeasureType.denominator)<<")";
00153   return out;
00154 }
00155 
00156 
00160 template<class Scalar>
00161 class ReductionFunctional : public Teuchos::Describable {
00162 public:
00163 
00166 
00175   typename ScalarTraits<Scalar>::magnitudeType
00176   reduce( const VectorBase<Scalar> &v ) const
00177     {
00178 #ifdef THYRA_DEBUG
00179       TEUCHOS_TEST_FOR_EXCEPTION(!isCompatible(v), Exceptions::IncompatibleVectorSpaces,
00180         "Error, the vector v="<<v.description()<<" is not compatiable with"
00181         " *this="<<this->description()<<"!");
00182 #endif
00183       return reduceImpl(v);
00184     }      
00185 
00189   bool isCompatible( const VectorBase<Scalar> &v ) const
00190     { return isCompatibleImpl(v); }
00191 
00193 
00194 protected:
00195 
00198 
00200   virtual typename ScalarTraits<Scalar>::magnitudeType
00201   reduceImpl( const VectorBase<Scalar> &v ) const = 0;
00202 
00204   virtual bool isCompatibleImpl( const VectorBase<Scalar> &v ) const = 0;
00205 
00207 
00208 };
00209 
00210 
00311 template <class Scalar>
00312 struct SolveCriteria {
00314   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00316   static ScalarMag unspecifiedTolerance() { return ScalarMag(-1.0); }
00319   SolveMeasureType solveMeasureType;
00322   ScalarMag requestedTol;
00328   RCP<ParameterList> extraParameters;
00331   RCP<const ReductionFunctional<Scalar> > numeratorReductionFunc;
00334   RCP<const ReductionFunctional<Scalar> > denominatorReductionFunc;
00336   SolveCriteria()
00337     : requestedTol(unspecifiedTolerance())
00338     {}
00340   SolveCriteria(
00341     SolveMeasureType solveMeasureType_in,
00342     ScalarMag requestedTol_in,
00343     const RCP<ParameterList> &extraParameters_in = Teuchos::null,
00344     const RCP<ReductionFunctional<Scalar> > &numeratorReductionFunc_in = Teuchos::null,
00345     const RCP<ReductionFunctional<Scalar> > &denominatorReductionFunc_in = Teuchos::null
00346     )
00347     : solveMeasureType(solveMeasureType_in),
00348       requestedTol(requestedTol_in), 
00349       extraParameters(extraParameters_in),
00350       numeratorReductionFunc(numeratorReductionFunc_in),
00351       denominatorReductionFunc(denominatorReductionFunc_in)
00352     {}
00353 };
00354 
00355 
00360 template<class Scalar>
00361 std::ostream& operator<<(std::ostream &out, const SolveCriteria<Scalar> &solveCriteria)
00362 {
00363   out << typeName(solveCriteria) << "{";
00364   out << "solveMeasureType="<<solveCriteria.solveMeasureType;
00365   out << ", requestedTol="<<solveCriteria.requestedTol;
00366   if (nonnull(solveCriteria.extraParameters)) {
00367     out << ", extraParameters="<<solveCriteria.extraParameters;
00368   }
00369   if (nonnull(solveCriteria.numeratorReductionFunc)) {
00370     out << ", numeratorReductionFunc="<<solveCriteria.numeratorReductionFunc->description();
00371   }
00372   if (nonnull(solveCriteria.denominatorReductionFunc)) {
00373     out << ", denominatorReductionFunc="<<solveCriteria.denominatorReductionFunc->description();
00374   }
00375   out << "}";
00376   return out;
00377 }
00378 
00379 
00384 class CatastrophicSolveFailure : public std::runtime_error
00385 {public: CatastrophicSolveFailure(const std::string& what_arg) : std::runtime_error(what_arg) {}};
00386 
00387 
00392 enum ESolveStatus {
00393   SOLVE_STATUS_CONVERGED        
00394   ,SOLVE_STATUS_UNCONVERGED     
00395   ,SOLVE_STATUS_UNKNOWN         
00396 };
00397 
00398 
00403 inline
00404 const std::string toString(const ESolveStatus solveStatus)
00405 {
00406   switch(solveStatus) {
00407     case SOLVE_STATUS_CONVERGED:    return "SOLVE_STATUS_CONVERGED";
00408     case SOLVE_STATUS_UNCONVERGED:  return "SOLVE_STATUS_UNCONVERGED";
00409     case SOLVE_STATUS_UNKNOWN:      return "SOLVE_STATUS_UNKNOWN";
00410     default: TEUCHOS_TEST_FOR_EXCEPT(true);
00411   }
00412   return ""; // Never be called!
00413 }
00414 
00415 
00422 template <class Scalar>
00423 struct SolveStatus {
00425   typedef typename Teuchos::ScalarTraits<Scalar>::magnitudeType ScalarMag;
00427   static ScalarMag unknownTolerance() { return ScalarMag(-1); }
00429   ESolveStatus solveStatus;
00433   ScalarMag achievedTol;
00435   std::string message;
00438   RCP<ParameterList> extraParameters;
00440   SolveStatus()
00441     :solveStatus(SOLVE_STATUS_UNKNOWN), achievedTol(unknownTolerance())
00442     {}
00445   static std::string achievedTolToString( const ScalarMag &achievedTol )
00446     {
00447       if(achievedTol==unknownTolerance()) return "unknownTolerance()";
00448       std::ostringstream oss; oss << achievedTol; return oss.str();
00449     }
00450 };
00451 
00452 
00457 template <class Scalar>
00458 std::ostream& operator<<( std::ostream& out_arg, const SolveStatus<Scalar> &solveStatus )
00459 {
00460   RCP<Teuchos::FancyOStream>
00461     out = Teuchos::getFancyOStream(Teuchos::rcp(&out_arg,false));
00462   Teuchos::OSTab tab(out);
00463   *out
00464     << "solveStatus = " << toString(solveStatus.solveStatus) << std::endl
00465     << "achievedTol = " << SolveStatus<Scalar>::achievedTolToString(solveStatus.achievedTol) << std::endl;
00466   *out << "message:";
00467   if (solveStatus.message.length()) {
00468     Teuchos::OSTab tab2(out);
00469     *out << "\n" << solveStatus.message << "\n";
00470   }
00471   *out << "extraParameters:";
00472   if(solveStatus.extraParameters.get()) {
00473     *out << "\n";
00474     Teuchos::OSTab tab3(out);
00475     solveStatus.extraParameters->print(*out, 1000, true);
00476   }
00477   else {
00478     *out << " NONE\n";
00479   }
00480   return out_arg;
00481 }
00482 
00483 
00489 enum ESupportSolveUse {
00490   SUPPORT_SOLVE_UNSPECIFIED  
00491   ,SUPPORT_SOLVE_FORWARD_ONLY  
00492   ,SUPPORT_SOLVE_TRANSPOSE_ONLY  
00493   ,SUPPORT_SOLVE_FORWARD_AND_TRANSPOSE  
00494 };
00495 
00496 
00501 enum EPreconditionerInputType {
00502   PRECONDITIONER_INPUT_TYPE_AS_OPERATOR  
00503   ,PRECONDITIONER_INPUT_TYPE_AS_MATRIX   
00504 };
00505 
00506 
00511 template <class Scalar>
00512 void accumulateSolveStatusInit(
00513   const Ptr<SolveStatus<Scalar> > &overallSolveStatus
00514   )
00515 {
00516   overallSolveStatus->solveStatus = SOLVE_STATUS_CONVERGED;
00517 }
00518 
00519 
00536 template <class Scalar>
00537 void accumulateSolveStatus(
00538   const SolveCriteria<Scalar>, // ToDo: Never used, need to take this out!
00539   const SolveStatus<Scalar> &solveStatus,
00540   const Ptr<SolveStatus<Scalar> > &overallSolveStatus
00541   )
00542 {
00543   switch(solveStatus.solveStatus) {
00544     case SOLVE_STATUS_UNCONVERGED:
00545     {
00546       // First, if we see any unconverged solve status, then the entire block is
00547       // unconverged!
00548       overallSolveStatus->solveStatus = SOLVE_STATUS_UNCONVERGED;
00549       overallSolveStatus->message = solveStatus.message;
00550       overallSolveStatus->extraParameters = solveStatus.extraParameters;
00551       break;
00552     }
00553     case SOLVE_STATUS_UNKNOWN:
00554     {
00555       // Next, if any solve status is unknown, then if the overall solve
00556       // status says converged, then we have to mark it as unknown.  Note that
00557       // unknown could mean that the system is actually converged!
00558       switch(overallSolveStatus->solveStatus) {
00559         case SOLVE_STATUS_CONVERGED:
00560           overallSolveStatus->solveStatus = SOLVE_STATUS_UNKNOWN;
00561           break;
00562         case SOLVE_STATUS_UNCONVERGED:
00563         case SOLVE_STATUS_UNKNOWN:
00564           // If we get here then the overall solve status is either unknown
00565           // already or says unconverged and this will not change here!
00566           overallSolveStatus->message = solveStatus.message;
00567           overallSolveStatus->extraParameters = solveStatus.extraParameters;
00568           break;
00569         default:
00570           TEUCHOS_TEST_FOR_EXCEPT(true); // Corrupted enum?
00571       }
00572       break;
00573     }
00574     case SOLVE_STATUS_CONVERGED:
00575     {
00576       // If we get here then the overall solve status is either unknown,
00577       // unconverged, or converged and this will not change here!
00578       if(overallSolveStatus->message == "")
00579         overallSolveStatus->message = solveStatus.message;
00580       break;
00581     }
00582     default:
00583       TEUCHOS_TEST_FOR_EXCEPT(true); // Corrupted enum?
00584   }
00585   // Update the achieved tolerence to the maximum returned
00586   if( solveStatus.achievedTol > overallSolveStatus->achievedTol ) {
00587     overallSolveStatus->achievedTol = solveStatus.achievedTol;
00588   }
00589   // Set a message if none is set
00590   if(overallSolveStatus->message == "")
00591     overallSolveStatus->message = solveStatus.message;
00592   // Set the extra parameters if none is set
00593   if(overallSolveStatus->extraParameters.get()==NULL)
00594     overallSolveStatus->extraParameters = solveStatus.extraParameters;
00595 }
00596 
00597 
00598 } // namespace Thyra
00599 
00600 
00601 #endif // THYRA_SOLVE_SUPPORT_TYPES_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines