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