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

Generated on Thu Sep 18 12:33:02 2008 for Thyra Package Browser (Single Doxygen Collection) by doxygen 1.3.9.1