Thyra Version of the Day
Thyra_DefaultInverseLinearOp_def.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_DEFAULT_INVERSE_LINEAR_OP_DEF_HPP
00043 #define THYRA_DEFAULT_INVERSE_LINEAR_OP_DEF_HPP
00044 
00045 #include "Thyra_DefaultInverseLinearOp_decl.hpp"
00046 #include "Thyra_MultiVectorStdOps.hpp"
00047 #include "Thyra_AssertOp.hpp"
00048 #include "Teuchos_Utils.hpp"
00049 #include "Teuchos_TypeNameTraits.hpp"
00050 
00051 
00052 namespace Thyra {
00053 
00054 
00055 // Constructors/initializers/accessors
00056 
00057 
00058 template<class Scalar>
00059 DefaultInverseLinearOp<Scalar>::DefaultInverseLinearOp()
00060 {}
00061 
00062 
00063 template<class Scalar>
00064 DefaultInverseLinearOp<Scalar>::DefaultInverseLinearOp(
00065   const Teuchos::RCP<LinearOpWithSolveBase<Scalar> > &lows,
00066   const SolveCriteria<Scalar> *fwdSolveCriteria,
00067   const EThrowOnSolveFailure throwOnFwdSolveFailure,
00068   const SolveCriteria<Scalar> *adjSolveCriteria,
00069   const EThrowOnSolveFailure throwOnAdjSolveFailure
00070   )
00071 {
00072   initializeImpl(
00073     lows,fwdSolveCriteria,throwOnFwdSolveFailure
00074     ,adjSolveCriteria,throwOnAdjSolveFailure
00075     );
00076 }
00077 
00078 
00079 template<class Scalar>
00080 DefaultInverseLinearOp<Scalar>::DefaultInverseLinearOp(
00081   const Teuchos::RCP<const LinearOpWithSolveBase<Scalar> > &lows,
00082   const SolveCriteria<Scalar> *fwdSolveCriteria,
00083   const EThrowOnSolveFailure throwOnFwdSolveFailure,
00084   const SolveCriteria<Scalar> *adjSolveCriteria,
00085   const EThrowOnSolveFailure throwOnAdjSolveFailure
00086   )
00087 {
00088   initializeImpl(
00089     lows,fwdSolveCriteria,throwOnFwdSolveFailure
00090     ,adjSolveCriteria,throwOnAdjSolveFailure
00091     );
00092 }
00093 
00094 
00095 template<class Scalar>
00096 void DefaultInverseLinearOp<Scalar>::initialize(
00097   const Teuchos::RCP<LinearOpWithSolveBase<Scalar> > &lows,
00098   const SolveCriteria<Scalar> *fwdSolveCriteria,
00099   const EThrowOnSolveFailure throwOnFwdSolveFailure,
00100   const SolveCriteria<Scalar> *adjSolveCriteria,
00101   const EThrowOnSolveFailure throwOnAdjSolveFailure
00102   )
00103 {
00104   initializeImpl(
00105     lows,fwdSolveCriteria,throwOnFwdSolveFailure
00106     ,adjSolveCriteria,throwOnAdjSolveFailure
00107     );
00108 }
00109 
00110 
00111 template<class Scalar>
00112 void DefaultInverseLinearOp<Scalar>::initialize(
00113   const Teuchos::RCP<const LinearOpWithSolveBase<Scalar> > &lows,
00114   const SolveCriteria<Scalar> *fwdSolveCriteria,
00115   const EThrowOnSolveFailure throwOnFwdSolveFailure,
00116   const SolveCriteria<Scalar> *adjSolveCriteria,
00117   const EThrowOnSolveFailure throwOnAdjSolveFailure
00118   )
00119 {
00120   initializeImpl(
00121     lows,fwdSolveCriteria,throwOnFwdSolveFailure
00122     ,adjSolveCriteria,throwOnAdjSolveFailure
00123     );
00124 }
00125 
00126 
00127 template<class Scalar>
00128 void DefaultInverseLinearOp<Scalar>::uninitialize()
00129 {
00130   lows_.uninitialize();
00131   fwdSolveCriteria_ = Teuchos::null;
00132   adjSolveCriteria_ = Teuchos::null;
00133 }
00134 
00135 
00136 // Overridden form InverseLinearOpBase
00137 
00138 
00139 template<class Scalar>
00140 bool DefaultInverseLinearOp<Scalar>::isLowsConst() const
00141 {
00142   return lows_.isConst();
00143 }
00144 
00145 
00146 template<class Scalar>
00147 Teuchos::RCP<LinearOpWithSolveBase<Scalar> >
00148 DefaultInverseLinearOp<Scalar>::getNonconstLows()
00149 {
00150   return lows_.getNonconstObj();
00151 }
00152 
00153 
00154 template<class Scalar>
00155 Teuchos::RCP<const LinearOpWithSolveBase<Scalar> >
00156 DefaultInverseLinearOp<Scalar>::getLows() const
00157 {
00158   return lows_.getConstObj();
00159 }
00160 
00161 
00162 // Overridden from LinearOpBase
00163 
00164 
00165 template<class Scalar>
00166 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00167 DefaultInverseLinearOp<Scalar>::range() const
00168 {
00169   assertInitialized();
00170   return lows_.getConstObj()->domain();
00171 }
00172 
00173 
00174 template<class Scalar>
00175 Teuchos::RCP< const VectorSpaceBase<Scalar> >
00176 DefaultInverseLinearOp<Scalar>::domain() const
00177 {
00178   assertInitialized();
00179   return lows_.getConstObj()->range();
00180 }
00181 
00182 
00183 template<class Scalar>
00184 Teuchos::RCP<const LinearOpBase<Scalar> >
00185 DefaultInverseLinearOp<Scalar>::clone() const
00186 {
00187   return Teuchos::null; // Not supported yet but could be!
00188 }
00189 
00190 
00191 // Overridden from Teuchos::Describable
00192 
00193                                                 
00194 template<class Scalar>
00195 std::string DefaultInverseLinearOp<Scalar>::description() const
00196 {
00197   assertInitialized();
00198   std::ostringstream oss;
00199   oss
00200     << Teuchos::Describable::description() << "{"
00201     << "lows="<<lows_.getConstObj()->description()
00202     << ",fwdSolveCriteria="<<(fwdSolveCriteria_.get()?"...":"DEFAULT")
00203     << ",adjSolveCriteria="<<(adjSolveCriteria_.get()?"...":"DEFAULT")
00204     << "}";
00205   return oss.str();
00206 }
00207 
00208 
00209 template<class Scalar>
00210 void DefaultInverseLinearOp<Scalar>::describe(
00211   Teuchos::FancyOStream &out,
00212   const Teuchos::EVerbosityLevel verbLevel
00213   ) const
00214 {
00215   using Teuchos::RCP;
00216   using Teuchos::OSTab;
00217   assertInitialized();
00218   OSTab tab(out);
00219   switch(verbLevel) {
00220     case Teuchos::VERB_DEFAULT:
00221     case Teuchos::VERB_LOW:
00222       out << this->description() << std::endl;
00223       break;
00224     case Teuchos::VERB_MEDIUM:
00225     case Teuchos::VERB_HIGH:
00226     case Teuchos::VERB_EXTREME:
00227     {
00228       out
00229         << Teuchos::Describable::description() << "{"
00230         << "rangeDim=" << this->range()->dim()
00231         << ",domainDim=" << this->domain()->dim() << "}:\n";
00232       OSTab tab2(out);
00233       out <<  "lows = ";
00234       if(!lows_.getConstObj().get()) {
00235         out << " NULL\n";
00236       }
00237       else {
00238         out << Teuchos::describe(*lows_.getConstObj(),verbLevel);
00239       }
00240       break;
00241     }
00242     default:
00243       TEUCHOS_TEST_FOR_EXCEPT(true); // Should never be called!
00244   }
00245 }
00246 
00247 
00248 // protected
00249 
00250 
00251 // Overridden from LinearOpBase
00252 
00253 
00254 template<class Scalar>
00255 bool DefaultInverseLinearOp<Scalar>::opSupportedImpl(EOpTransp M_trans) const
00256 {
00257   if (nonnull(lows_)) {
00258     return solveSupports(*lows_.getConstObj(),M_trans);
00259   }
00260   return false;
00261 }
00262 
00263 
00264 template<class Scalar>
00265 void DefaultInverseLinearOp<Scalar>::applyImpl(
00266   const EOpTransp M_trans,
00267   const MultiVectorBase<Scalar> &X,
00268   const Ptr<MultiVectorBase<Scalar> > &Y,
00269   const Scalar alpha,
00270   const Scalar beta
00271   ) const
00272 {
00273   typedef Teuchos::ScalarTraits<Scalar> ST;
00274   assertInitialized();
00275   // ToDo: Put in hooks for propogating verbosity level
00276   //
00277   // Y = alpha*op(M)*X + beta*Y
00278   //
00279   //   =>
00280   //
00281   // Y = beta*Y
00282   // Y += alpha*inv(op(lows))*X
00283   //
00284   Teuchos::RCP<MultiVectorBase<Scalar> > T;
00285   if(beta==ST::zero()) {
00286     T = Teuchos::rcpFromPtr(Y);
00287   }
00288   else {
00289     T = createMembers(Y->range(),Y->domain()->dim());
00290     scale(beta, Y);
00291   }
00292   //
00293   const Ptr<const SolveCriteria<Scalar> > solveCriteria = 
00294     (
00295       real_trans(M_trans)==NOTRANS
00296       ? fwdSolveCriteria_.ptr()
00297       : adjSolveCriteria_.ptr()
00298       );
00299   assign(T.ptr(), ST::zero()); // Have to initialize before solve!
00300   SolveStatus<Scalar> solveStatus =
00301     Thyra::solve<Scalar>(*lows_.getConstObj(), M_trans, X, T.ptr(), solveCriteria);
00302 
00303   TEUCHOS_TEST_FOR_EXCEPTION(
00304     nonnull(solveCriteria) && solveStatus.solveStatus!=SOLVE_STATUS_CONVERGED
00305     && ( real_trans(M_trans)==NOTRANS
00306          ? throwOnFwdSolveFailure_==THROW_ON_SOLVE_FAILURE
00307          : throwOnAdjSolveFailure_==THROW_ON_SOLVE_FAILURE )
00308     ,CatastrophicSolveFailure
00309     ,"Error, the LOWS object " << lows_.getConstObj()->description() << " returned an unconverged"
00310     "status of " << toString(solveStatus.solveStatus) << " with the message "
00311     << solveStatus.message << "."
00312     );
00313   //
00314   if(beta==ST::zero()) {
00315     scale(alpha, Y);
00316   }
00317   else {
00318     update( alpha, *T, Y );
00319   }
00320 }
00321 
00322 
00323 // private
00324 
00325 
00326 template<class Scalar>
00327 template<class LOWS>
00328 void DefaultInverseLinearOp<Scalar>::initializeImpl(
00329   const Teuchos::RCP<LOWS> &lows,
00330   const SolveCriteria<Scalar> *fwdSolveCriteria,
00331   const EThrowOnSolveFailure throwOnFwdSolveFailure,
00332   const SolveCriteria<Scalar> *adjSolveCriteria,
00333   const EThrowOnSolveFailure throwOnAdjSolveFailure
00334   )
00335 {
00336   lows_.initialize(lows);
00337   if(fwdSolveCriteria)
00338     fwdSolveCriteria_ = Teuchos::rcp(new SolveCriteria<Scalar>(*fwdSolveCriteria));
00339   else
00340     fwdSolveCriteria_ = Teuchos::null;
00341   if(adjSolveCriteria)
00342     adjSolveCriteria_ = Teuchos::rcp(new SolveCriteria<Scalar>(*adjSolveCriteria));
00343   else
00344     adjSolveCriteria_ = Teuchos::null;
00345   throwOnFwdSolveFailure_ = throwOnFwdSolveFailure;
00346   throwOnAdjSolveFailure_ = throwOnAdjSolveFailure;
00347   const std::string lowsLabel = lows_.getConstObj()->getObjectLabel();
00348   if(lowsLabel.length())
00349     this->setObjectLabel( "inv("+lowsLabel+")" );
00350 }
00351 
00352 
00353 } // end namespace Thyra
00354 
00355 
00356 // Related non-member functions
00357 
00358 
00359 template<class Scalar>
00360 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00361 Thyra::nonconstInverse(
00362   const RCP<LinearOpWithSolveBase<Scalar> > &A,
00363   const Ptr<const SolveCriteria<Scalar> > &fwdSolveCriteria,
00364   const EThrowOnSolveFailure throwOnFwdSolveFailure,
00365   const Ptr<const SolveCriteria<Scalar> > &adjSolveCriteria,
00366   const EThrowOnSolveFailure throwOnAdjSolveFailure
00367   )
00368 {
00369   return Teuchos::rcp(
00370     new DefaultInverseLinearOp<Scalar>(
00371       A, fwdSolveCriteria.get(), throwOnFwdSolveFailure,
00372       adjSolveCriteria.get(), throwOnAdjSolveFailure
00373       )
00374     );
00375 }
00376 
00377 template<class Scalar>
00378 Teuchos::RCP<Thyra::LinearOpBase<Scalar> >
00379 Thyra::inverse(
00380   const RCP<const LinearOpWithSolveBase<Scalar> > &A,
00381   const Ptr<const SolveCriteria<Scalar> > &fwdSolveCriteria,
00382   const EThrowOnSolveFailure throwOnFwdSolveFailure,
00383   const Ptr<const SolveCriteria<Scalar> > &adjSolveCriteria,
00384   const EThrowOnSolveFailure throwOnAdjSolveFailure
00385   )
00386 {
00387   return Teuchos::rcp(
00388     new DefaultInverseLinearOp<Scalar>(
00389       A, fwdSolveCriteria.get(), throwOnFwdSolveFailure,
00390       adjSolveCriteria.get(), throwOnAdjSolveFailure
00391       )
00392     );
00393 }
00394 
00395 
00396 //
00397 // Explicit instantiation macro
00398 //
00399 // Must be expanded from within the Thyra namespace!
00400 //
00401 
00402 
00403 #define THYRA_DEFAULT_INVERSE_LINEAR_OP_INSTANT(SCALAR) \
00404   \
00405   template class DefaultInverseLinearOp<SCALAR >; \
00406    \
00407   template RCP<LinearOpBase<SCALAR > > \
00408   nonconstInverse( \
00409     const RCP<LinearOpWithSolveBase<SCALAR > > &A, \
00410     const Ptr<const SolveCriteria<SCALAR > > &fwdSolveCriteria, \
00411     const EThrowOnSolveFailure throwOnFwdSolveFailure, \
00412     const Ptr<const SolveCriteria<SCALAR > > &adjSolveCriteria, \
00413     const EThrowOnSolveFailure throwOnAdjSolveFailure \
00414     ); \
00415    \
00416   template RCP<LinearOpBase<SCALAR > > \
00417   inverse( \
00418     const RCP<const LinearOpWithSolveBase<SCALAR > > &A, \
00419     const Ptr<const SolveCriteria<SCALAR > > &fwdSolveCriteria, \
00420     const EThrowOnSolveFailure throwOnFwdSolveFailure, \
00421     const Ptr<const SolveCriteria<SCALAR > > &adjSolveCriteria, \
00422     const EThrowOnSolveFailure throwOnAdjSolveFailure \
00423     );
00424 
00425 
00426 #endif  // THYRA_DEFAULT_INVERSE_LINEAR_OP_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines