Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_ModelEvaluatorDefaultBase.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_MODEL_EVALUATOR_DEFAULT_BASE_HPP
00030 #define THYRA_MODEL_EVALUATOR_DEFAULT_BASE_HPP
00031 
00032 #include "Thyra_VectorBase.hpp"
00033 
00034 #include "Thyra_ModelEvaluator.hpp"
00035 #include "Thyra_LinearOpWithSolveFactoryHelpers.hpp"
00036 
00037 
00038 #ifdef HAVE_THYRA_ME_POLYNOMIAL
00039 
00040 
00041 // Define the polynomial traits class specializtaion
00042 // Teuchos::PolynomialTraits<VectorBase > before there is any chance of an
00043 // instantiation requiring the definition of this class.  The Intel C++ 9.1
00044 // compiler is having trouble compiling Thyra_EpetraModelEvaluator.cpp because
00045 // Thyra_ModelEvaluatorBase_decl.hpp is only including
00046 // "Teuchos_Polynomial.hpp" which does not contain the needed specialization.
00047 // By including it here, all client code is likely to compile just fine.  I am
00048 // going through all of the in order to avoid having to put
00049 // Thyra_PolynomialVectorTraits.hpp in the interfaces directory.  The problem
00050 // with putting Thyra_PolynomialVectorTraits.hpp the interfaces directory is
00051 // that it contains calls to code in the support directory and creates an
00052 // unacceptable circular dependency that will cause problems once we move to
00053 // subpackages in the CMake build system.
00054 #include "Thyra_PolynomialVectorTraits.hpp"
00055 
00056 
00057 #endif // HAVE_THYRA_ME_POLYNOMIAL
00058 
00059 
00060 namespace Thyra {
00061 
00062 
00063 //
00064 // Undocumentated (from the user's perspective) types that are used to
00065 // implement ModelEvaluatorDefaultBase.  These types are defined outside of
00066 // the templated class since they do nt depend on the scalar type.
00067 //
00068 
00069 
00070 namespace ModelEvaluatorDefaultBaseTypes {
00071 
00072 
00073 // Type used to determine if the ModelEvaluatorDefaultBase implementation will
00074 // provide a LinearOpBase wrapper for derivative object given in
00075 // MultiVectorBase form.
00076 class DefaultDerivLinearOpSupport {
00077 public:
00078   DefaultDerivLinearOpSupport()
00079     :provideDefaultLinearOp_(false),
00080      mvImplOrientation_(ModelEvaluatorBase::DERIV_MV_BY_COL)
00081     {}
00082   DefaultDerivLinearOpSupport(
00083     const ModelEvaluatorBase::EDerivativeMultiVectorOrientation mvImplOrientation_in
00084     )
00085     :provideDefaultLinearOp_(true),
00086      mvImplOrientation_(mvImplOrientation_in)
00087     {}
00088   bool provideDefaultLinearOp() const
00089     { return provideDefaultLinearOp_; }
00090   ModelEvaluatorBase::EDerivativeMultiVectorOrientation mvImplOrientation() const
00091     { return mvImplOrientation_; }
00092 private:
00093   bool provideDefaultLinearOp_;
00094   ModelEvaluatorBase::EDerivativeMultiVectorOrientation mvImplOrientation_;
00095 };
00096 
00097 
00098 // Type used to determine if the ModelEvaluatorDefaultBase implementation will
00099 // provide an explicit transpose copy of a derivative object given in
00100 // MultiVectorBase form.
00101 class DefaultDerivMvAdjointSupport {
00102 public:
00103   DefaultDerivMvAdjointSupport()
00104     :provideDefaultAdjoint_(false),
00105      mvAdjointCopyOrientation_(ModelEvaluatorBase::DERIV_MV_BY_COL)
00106     {}
00107   DefaultDerivMvAdjointSupport(
00108     const ModelEvaluatorBase::EDerivativeMultiVectorOrientation mvAdjointCopyOrientation_in
00109     )
00110     :provideDefaultAdjoint_(true),
00111      mvAdjointCopyOrientation_(mvAdjointCopyOrientation_in)
00112     {}
00113   bool provideDefaultAdjoint() const
00114     { return provideDefaultAdjoint_; }
00115   ModelEvaluatorBase::EDerivativeMultiVectorOrientation mvAdjointCopyOrientation() const
00116     { return mvAdjointCopyOrientation_; }
00117 private:
00118   bool provideDefaultAdjoint_;
00119   ModelEvaluatorBase::EDerivativeMultiVectorOrientation mvAdjointCopyOrientation_;
00120 };
00121 
00122 
00123 // Type used to remember a pair of transposed multi-vectors to implement a
00124 // adjoint copy.
00125 template<class Scalar>
00126 struct MultiVectorAdjointPair {
00127   MultiVectorAdjointPair()
00128     {}
00129   MultiVectorAdjointPair(
00130     const RCP<MultiVectorBase<Scalar> > &in_mvOuter,
00131     const RCP<const MultiVectorBase<Scalar> > &in_mvImplAdjoint
00132     )
00133     : mvOuter(in_mvOuter),
00134       mvImplAdjoint(in_mvImplAdjoint)
00135     {}
00136   RCP<MultiVectorBase<Scalar> > mvOuter;
00137   RCP<const MultiVectorBase<Scalar> > mvImplAdjoint;
00138 private:
00139 };
00140 
00141 
00142 
00143 
00144 } // namespace ModelEvaluatorDefaultBaseTypes
00145 
00146 
00174 template<class Scalar>
00175 class ModelEvaluatorDefaultBase : virtual public ModelEvaluator<Scalar>
00176 {
00177 public:
00178 
00181 
00183   int Np() const;
00185   int Ng() const;
00187   RCP<LinearOpBase<Scalar> > create_DfDp_op(int l) const;
00189   RCP<LinearOpBase<Scalar> > create_DgDx_dot_op(int j) const;
00191   RCP<LinearOpBase<Scalar> > create_DgDx_op(int j) const;
00193   RCP<LinearOpBase<Scalar> > create_DgDp_op(int j, int l) const;
00195   RCP<LinearOpWithSolveBase<Scalar> > create_W() const;
00197   ModelEvaluatorBase::OutArgs<Scalar> createOutArgs() const;
00199   void evalModel(
00200     const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00201     const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00202     ) const;
00203 
00205 
00206 protected:
00207 
00210 
00219   void initializeDefaultBase();
00220 
00222 
00223 private:
00224 
00227 
00229   virtual RCP<LinearOpBase<Scalar> > create_DfDp_op_impl(int l) const;
00230 
00232   virtual RCP<LinearOpBase<Scalar> > create_DgDx_dot_op_impl(int j) const;
00233 
00235   virtual RCP<LinearOpBase<Scalar> > create_DgDx_op_impl(int j) const;
00236 
00238   virtual RCP<LinearOpBase<Scalar> > create_DgDp_op_impl(int j, int l) const;
00239   
00241 
00244 
00246   virtual ModelEvaluatorBase::OutArgs<Scalar> createOutArgsImpl() const = 0;
00247 
00249   virtual void evalModelImpl(
00250     const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00251     const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00252     ) const = 0;
00253 
00255 
00256 protected:
00257 
00259   ModelEvaluatorDefaultBase();
00260 
00261 private:
00262   
00263   // //////////////////////////////
00264   // Private tpyes
00265 
00266   typedef ModelEvaluatorDefaultBaseTypes::DefaultDerivLinearOpSupport
00267   DefaultDerivLinearOpSupport;
00268 
00269   typedef ModelEvaluatorDefaultBaseTypes::DefaultDerivMvAdjointSupport
00270   DefaultDerivMvAdjointSupport;
00271 
00272   typedef ModelEvaluatorDefaultBaseTypes::MultiVectorAdjointPair<Scalar>
00273   MultiVectorAdjointPair;
00274   
00275   // //////////////////////////////
00276   // Private data members
00277 
00278   bool isInitialized_;
00279 
00280   Array<DefaultDerivLinearOpSupport> DfDp_default_op_support_;
00281 
00282   Array<DefaultDerivLinearOpSupport> DgDx_dot_default_op_support_;
00283 
00284   Array<DefaultDerivLinearOpSupport> DgDx_default_op_support_;
00285 
00286   Array<Array<DefaultDerivLinearOpSupport> > DgDp_default_op_support_;
00287   Array<Array<DefaultDerivMvAdjointSupport> > DgDp_default_mv_support_;
00288 
00289   bool default_W_support_;
00290 
00291   ModelEvaluatorBase::OutArgs<Scalar> prototypeOutArgs_;
00292 
00293   // //////////////////////////////
00294   // Private member functions
00295 
00296   void lazyInitializeDefaultBase() const;
00297 
00298   void assert_l(const int l) const;
00299 
00300   void assert_j(const int j) const;
00301 
00302   // //////////////////////////////
00303   // Private static functions
00304 
00305   static DefaultDerivLinearOpSupport
00306   determineDefaultDerivLinearOpSupport(
00307     const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl
00308     );
00309 
00310   static RCP<LinearOpBase<Scalar> >
00311   createDefaultLinearOp(
00312     const DefaultDerivLinearOpSupport &defaultLinearOpSupport,
00313     const RCP<const VectorSpaceBase<Scalar> > &fnc_space,
00314     const RCP<const VectorSpaceBase<Scalar> > &var_space
00315     );
00316 
00317   static ModelEvaluatorBase::DerivativeSupport
00318   updateDefaultLinearOpSupport(
00319     const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
00320     const DefaultDerivLinearOpSupport &defaultLinearOpSupport
00321     );
00322 
00323   static ModelEvaluatorBase::Derivative<Scalar>
00324   getOutArgImplForDefaultLinearOpSupport(
00325     const ModelEvaluatorBase::Derivative<Scalar> &deriv,
00326     const DefaultDerivLinearOpSupport &defaultLinearOpSupport
00327     );
00328 
00329   static DefaultDerivMvAdjointSupport
00330   determineDefaultDerivMvAdjointSupport(
00331     const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
00332     const VectorSpaceBase<Scalar> &fnc_space,
00333     const VectorSpaceBase<Scalar> &var_space
00334     );
00335 
00336   static ModelEvaluatorBase::DerivativeSupport
00337   updateDefaultDerivMvAdjointSupport(
00338     const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
00339     const DefaultDerivMvAdjointSupport &defaultMvAdjointSupport
00340     );
00341   
00342 };
00343 
00344 
00345 } // namespace Thyra
00346 
00347 
00348 //
00349 // Implementations
00350 //
00351 
00352 
00353 #include "Thyra_ModelEvaluatorHelpers.hpp"
00354 #include "Thyra_DefaultScaledAdjointLinearOp.hpp"
00355 #include "Thyra_DetachedMultiVectorView.hpp"
00356 #include "Teuchos_Assert.hpp"
00357 
00358 
00359 namespace Thyra {
00360 
00361 
00362 // Overridden from ModelEvaluator
00363 
00364 
00365 template<class Scalar>
00366 int ModelEvaluatorDefaultBase<Scalar>::Np() const
00367 {
00368   lazyInitializeDefaultBase();
00369   return prototypeOutArgs_.Np();
00370 }
00371 
00372 
00373 template<class Scalar>
00374 int ModelEvaluatorDefaultBase<Scalar>::Ng() const
00375 {
00376   lazyInitializeDefaultBase();
00377   return prototypeOutArgs_.Ng();
00378 }
00379 
00380 
00381 template<class Scalar>
00382 RCP<LinearOpWithSolveBase<Scalar> >
00383 ModelEvaluatorDefaultBase<Scalar>::create_W() const
00384 {
00385   lazyInitializeDefaultBase();
00386   if (default_W_support_)
00387     return this->get_W_factory()->createOp();
00388   return Teuchos::null;
00389 }
00390 
00391 
00392 template<class Scalar>
00393 RCP<LinearOpBase<Scalar> >
00394 ModelEvaluatorDefaultBase<Scalar>::create_DfDp_op(int l) const
00395 {
00396   lazyInitializeDefaultBase();
00397 #ifdef TEUCHOS_DEBUG
00398   assert_l(l);
00399 #endif
00400   const DefaultDerivLinearOpSupport
00401     defaultLinearOpSupport = DfDp_default_op_support_[l];
00402   if (defaultLinearOpSupport.provideDefaultLinearOp()) {
00403     return createDefaultLinearOp(
00404       defaultLinearOpSupport,
00405       this->get_f_space(),
00406       this->get_p_space(l)
00407       );
00408   }
00409   return this->create_DfDp_op_impl(l);
00410 }
00411 
00412 
00413 template<class Scalar>
00414 RCP<LinearOpBase<Scalar> >
00415 ModelEvaluatorDefaultBase<Scalar>::create_DgDx_dot_op(int j) const
00416 {
00417   lazyInitializeDefaultBase();
00418 #ifdef TEUCHOS_DEBUG
00419   assert_j(j);
00420 #endif
00421   const DefaultDerivLinearOpSupport
00422     defaultLinearOpSupport = DgDx_dot_default_op_support_[j];
00423   if (defaultLinearOpSupport.provideDefaultLinearOp()) {
00424     return createDefaultLinearOp(
00425       defaultLinearOpSupport,
00426       this->get_g_space(j),
00427       this->get_x_space()
00428       );
00429   }
00430   return this->create_DgDx_dot_op_impl(j);
00431 }
00432 
00433 
00434 template<class Scalar>
00435 RCP<LinearOpBase<Scalar> >
00436 ModelEvaluatorDefaultBase<Scalar>::create_DgDx_op(int j) const
00437 {
00438   lazyInitializeDefaultBase();
00439 #ifdef TEUCHOS_DEBUG
00440   assert_j(j);
00441 #endif
00442   const DefaultDerivLinearOpSupport
00443     defaultLinearOpSupport = DgDx_default_op_support_[j];
00444   if (defaultLinearOpSupport.provideDefaultLinearOp()) {
00445     return createDefaultLinearOp(
00446       defaultLinearOpSupport,
00447       this->get_g_space(j),
00448       this->get_x_space()
00449       );
00450   }
00451   return this->create_DgDx_op_impl(j);
00452 }
00453 
00454 
00455 template<class Scalar>
00456 RCP<LinearOpBase<Scalar> >
00457 ModelEvaluatorDefaultBase<Scalar>::create_DgDp_op(int j, int l) const
00458 {
00459   lazyInitializeDefaultBase();
00460 #ifdef TEUCHOS_DEBUG
00461   assert_j(j);
00462   assert_l(l);
00463 #endif
00464   const DefaultDerivLinearOpSupport
00465     defaultLinearOpSupport = DgDp_default_op_support_[j][l];
00466   if (defaultLinearOpSupport.provideDefaultLinearOp()) {
00467     return createDefaultLinearOp(
00468       defaultLinearOpSupport,
00469       this->get_g_space(j),
00470       this->get_p_space(l)
00471       );
00472   }
00473   return this->create_DgDp_op_impl(j,l);
00474 }
00475 
00476 
00477 template<class Scalar>
00478 ModelEvaluatorBase::OutArgs<Scalar>
00479 ModelEvaluatorDefaultBase<Scalar>::createOutArgs() const
00480 {
00481   lazyInitializeDefaultBase();
00482   return prototypeOutArgs_;
00483 }
00484 
00485 
00486 template<class Scalar>
00487 void ModelEvaluatorDefaultBase<Scalar>::evalModel(
00488   const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00489   const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00490   ) const
00491 {
00492 
00493   using Teuchos::outArg;
00494   typedef ModelEvaluatorBase MEB;
00495 
00496   lazyInitializeDefaultBase();
00497 
00498   const int l_Np = outArgs.Np();
00499   const int l_Ng = outArgs.Ng();
00500 
00501   //
00502   // A) Assert that the inArgs and outArgs object match this class!
00503   //
00504 
00505 #ifdef TEUCHOS_DEBUG
00506   assertInArgsEvalObjects(*this,inArgs);
00507   assertOutArgsEvalObjects(*this,outArgs,&inArgs);
00508 #endif  
00509   
00510   //
00511   // B) Setup the OutArgs object for the underlying implementation's
00512   // evalModelImpl(...) function
00513   //
00514 
00515   MEB::OutArgs<Scalar> outArgsImpl = this->createOutArgsImpl();
00516   Array<MultiVectorAdjointPair> DgDp_temp_adjoint_copies;
00517 
00518   {
00519 
00520     outArgsImpl.setArgs(outArgs,true);
00521 
00522     // DfDp(l)
00523     if (outArgsImpl.supports(MEB::OUT_ARG_f)) {
00524       for ( int l = 0; l < l_Np; ++l ) {
00525         const DefaultDerivLinearOpSupport defaultLinearOpSupport =
00526           DfDp_default_op_support_[l];
00527         if (defaultLinearOpSupport.provideDefaultLinearOp()) {
00528           outArgsImpl.set_DfDp( l,
00529             getOutArgImplForDefaultLinearOpSupport(
00530               outArgs.get_DfDp(l), defaultLinearOpSupport
00531               )
00532             );
00533         }
00534         else {
00535           // DfDp(l) already set by outArgsImpl.setArgs(...)!
00536         }
00537       }
00538     }
00539 
00540     // DgDx_dot(j)
00541     for ( int j = 0; j < l_Ng; ++j ) {
00542       const DefaultDerivLinearOpSupport defaultLinearOpSupport =
00543         DgDx_dot_default_op_support_[j];
00544       if (defaultLinearOpSupport.provideDefaultLinearOp()) {
00545         outArgsImpl.set_DgDx_dot( j,
00546           getOutArgImplForDefaultLinearOpSupport(
00547             outArgs.get_DgDx_dot(j), defaultLinearOpSupport
00548             )
00549           );
00550       }
00551       else {
00552         // DgDx_dot(j) already set by outArgsImpl.setArgs(...)!
00553       }
00554     }
00555 
00556     // DgDx(j)
00557     for ( int j = 0; j < l_Ng; ++j ) {
00558       const DefaultDerivLinearOpSupport defaultLinearOpSupport =
00559         DgDx_default_op_support_[j];
00560       if (defaultLinearOpSupport.provideDefaultLinearOp()) {
00561         outArgsImpl.set_DgDx( j,
00562           getOutArgImplForDefaultLinearOpSupport(
00563             outArgs.get_DgDx(j), defaultLinearOpSupport
00564             )
00565           );
00566       }
00567       else {
00568         // DgDx(j) already set by outArgsImpl.setArgs(...)!
00569       }
00570     }
00571 
00572     // DgDp(j,l)
00573     for ( int j = 0; j < l_Ng; ++j ) {
00574       const Array<DefaultDerivLinearOpSupport> &DgDp_default_op_support_j =
00575         DgDp_default_op_support_[j];
00576       const Array<DefaultDerivMvAdjointSupport> &DgDp_default_mv_support_j =
00577         DgDp_default_mv_support_[j];
00578       for ( int l = 0; l < l_Np; ++l ) {
00579         const DefaultDerivLinearOpSupport defaultLinearOpSupport =
00580           DgDp_default_op_support_j[l];
00581         const DefaultDerivMvAdjointSupport defaultMvAdjointSupport =
00582           DgDp_default_mv_support_j[l];
00583         MEB::Derivative<Scalar> DgDp_j_l;
00584         if (!outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none())
00585           DgDp_j_l = outArgs.get_DgDp(j,l);
00586         if (
00587           defaultLinearOpSupport.provideDefaultLinearOp()
00588           && !is_null(DgDp_j_l.getLinearOp())
00589           )
00590         {
00591           outArgsImpl.set_DgDp( j, l,
00592             getOutArgImplForDefaultLinearOpSupport(
00593               DgDp_j_l, defaultLinearOpSupport
00594               )
00595             );
00596         }
00597         else if (
00598           defaultMvAdjointSupport.provideDefaultAdjoint()
00599           && !is_null(DgDp_j_l.getMultiVector())
00600           )
00601         {
00602           const RCP<MultiVectorBase<Scalar> > DgDp_j_l_mv = 
00603             DgDp_j_l.getMultiVector();
00604           if (
00605             defaultMvAdjointSupport.mvAdjointCopyOrientation()
00606             ==
00607             DgDp_j_l.getMultiVectorOrientation()
00608             )
00609           {
00610             // The orientation of the multi-vector is different so we need to
00611             // create a temporary copy to pass to evalModelImpl(...) and then
00612             // copy it back again!
00613             const RCP<MultiVectorBase<Scalar> > DgDp_j_l_mv_adj =
00614               createMembers(DgDp_j_l_mv->domain(), DgDp_j_l_mv->range()->dim());
00615             outArgsImpl.set_DgDp( j, l,
00616               MEB::Derivative<Scalar>(
00617                 DgDp_j_l_mv_adj,
00618                 getOtherDerivativeMultiVectorOrientation(
00619                   defaultMvAdjointSupport.mvAdjointCopyOrientation()
00620                   )
00621                 )
00622               );
00623             // Remember these multi-vectors so that we can do the transpose copy
00624             // back after the evaluation!
00625             DgDp_temp_adjoint_copies.push_back(
00626               MultiVectorAdjointPair(DgDp_j_l_mv, DgDp_j_l_mv_adj)
00627               );
00628           }
00629           else {
00630             // The form of the multi-vector is supported by evalModelImpl(..)
00631             // and is already set on the outArgsImpl object.
00632           }
00633         }
00634         else {
00635           // DgDp(j,l) already set by outArgsImpl.setArgs(...)!
00636         }
00637       }
00638     }
00639 
00640     // W
00641     {
00642       RCP<LinearOpWithSolveBase<Scalar> > W;
00643       if ( default_W_support_ && !is_null(W=outArgs.get_W()) ) {
00644         const RCP<const LinearOpWithSolveFactoryBase<Scalar> >
00645           W_factory = this->get_W_factory();
00646         // Extract the underlying W_op object (if it exists)
00647         RCP<const LinearOpBase<Scalar> > W_op_const;
00648         uninitializeOp<Scalar>(*W_factory, W.ptr(), outArg(W_op_const));
00649         RCP<LinearOpBase<Scalar> > W_op;
00650         if (!is_null(W_op_const)) {
00651           // Here we remove the const.  This is perfectly safe since
00652           // w.r.t. this class, we put a non-const object in there and we can
00653           // expect to change that object after the fact.  That is our
00654           // prerogative.
00655           W_op = Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(W_op_const);
00656         }
00657         else {
00658           // The W_op object has not been initialized yet so create it.  The
00659           // next time through, we should not have to do this!
00660           W_op = this->create_W_op();
00661         }
00662         outArgsImpl.set_W_op(W_op);
00663       }
00664     }
00665     
00666   }
00667 
00668   //
00669   // C) Evaluate the underlying model implementation!
00670   //
00671 
00672   this->evalModelImpl( inArgs, outArgsImpl );
00673 
00674   //
00675   // D) Post-process the output arguments
00676   //
00677 
00678   // Do explicit transposes for DgDp(j,l) if needed
00679   const int numMvAdjointCopies = DgDp_temp_adjoint_copies.size();
00680   for ( int adj_copy_i = 0; adj_copy_i < numMvAdjointCopies; ++adj_copy_i ) {
00681     const MultiVectorAdjointPair adjPair =
00682       DgDp_temp_adjoint_copies[adj_copy_i];
00683     doExplicitMultiVectorAdjoint( *adjPair.mvImplAdjoint, &*adjPair.mvOuter );
00684   }
00685   
00686   // Update W given W_op and W_factory
00687   {
00688     RCP<LinearOpWithSolveBase<Scalar> > W;
00689     if ( default_W_support_ && !is_null(W=outArgs.get_W()) ) {
00690       const RCP<const LinearOpWithSolveFactoryBase<Scalar> >
00691         W_factory = this->get_W_factory();
00692       W_factory->setOStream(this->getOStream());
00693       W_factory->setVerbLevel(this->getVerbLevel());
00694       initializeOp<Scalar>(*W_factory, outArgsImpl.get_W_op().getConst(), W.ptr());
00695     }
00696   }
00697   
00698 }
00699 
00700 
00701 // protected
00702 
00703 
00704 // Setup functions called by subclasses
00705 
00706 template<class Scalar>
00707 void ModelEvaluatorDefaultBase<Scalar>::initializeDefaultBase()
00708 {
00709 
00710   typedef ModelEvaluatorBase MEB;
00711 
00712   // In case we throw half way thorugh, set to uninitialized
00713   isInitialized_ = false;
00714   default_W_support_ = false;
00715 
00716   //
00717   // A) Get the InArgs and OutArgs from the subclass
00718   //
00719   
00720   const MEB::InArgs<Scalar> inArgs = this->createInArgs();
00721   const MEB::OutArgs<Scalar> outArgsImpl = this->createOutArgsImpl();
00722 
00723   //
00724   // B) Validate the subclasses InArgs and OutArgs
00725   //
00726 
00727 #ifdef TEUCHOS_DEBUG
00728   assertInArgsOutArgsSetup( this->description(), inArgs, outArgsImpl );
00729 #endif // TEUCHOS_DEBUG
00730 
00731   //
00732   // C) Set up support for default derivative objects and prototype OutArgs
00733   //
00734 
00735   const int l_Ng = outArgsImpl.Ng();
00736   const int l_Np = outArgsImpl.Np();
00737 
00738   // Set support for all outputs supported in the underly implementation
00739   MEB::OutArgsSetup<Scalar> outArgs;
00740   outArgs.setModelEvalDescription(this->description());
00741   outArgs.set_Np_Ng(l_Np,l_Ng);
00742   outArgs.setSupports(outArgsImpl);
00743 
00744   // DfDp
00745   DfDp_default_op_support_.clear();
00746   if (outArgs.supports(MEB::OUT_ARG_f)) {
00747     for ( int l = 0; l < l_Np; ++l ) {
00748       const MEB::DerivativeSupport DfDp_l_impl_support =
00749         outArgsImpl.supports(MEB::OUT_ARG_DfDp,l);
00750       const DefaultDerivLinearOpSupport DfDp_l_op_support =
00751         determineDefaultDerivLinearOpSupport(DfDp_l_impl_support);
00752       DfDp_default_op_support_.push_back(DfDp_l_op_support);
00753       outArgs.setSupports(
00754         MEB::OUT_ARG_DfDp, l,
00755         updateDefaultLinearOpSupport(
00756           DfDp_l_impl_support, DfDp_l_op_support
00757           )
00758         );
00759     }
00760   }
00761 
00762   // DgDx_dot
00763   DgDx_dot_default_op_support_.clear();
00764   for ( int j = 0; j < l_Ng; ++j ) {
00765     const MEB::DerivativeSupport DgDx_dot_j_impl_support =
00766       outArgsImpl.supports(MEB::OUT_ARG_DgDx_dot,j);
00767     const DefaultDerivLinearOpSupport DgDx_dot_j_op_support =
00768       determineDefaultDerivLinearOpSupport(DgDx_dot_j_impl_support);
00769     DgDx_dot_default_op_support_.push_back(DgDx_dot_j_op_support);
00770     outArgs.setSupports(
00771       MEB::OUT_ARG_DgDx_dot, j,
00772       updateDefaultLinearOpSupport(
00773         DgDx_dot_j_impl_support, DgDx_dot_j_op_support
00774         )
00775       );
00776   }
00777 
00778   // DgDx
00779   DgDx_default_op_support_.clear();
00780   for ( int j = 0; j < l_Ng; ++j ) {
00781     const MEB::DerivativeSupport DgDx_j_impl_support =
00782       outArgsImpl.supports(MEB::OUT_ARG_DgDx,j);
00783     const DefaultDerivLinearOpSupport DgDx_j_op_support =
00784       determineDefaultDerivLinearOpSupport(DgDx_j_impl_support);
00785     DgDx_default_op_support_.push_back(DgDx_j_op_support);
00786     outArgs.setSupports(
00787       MEB::OUT_ARG_DgDx, j,
00788       updateDefaultLinearOpSupport(
00789         DgDx_j_impl_support, DgDx_j_op_support
00790         )
00791       );
00792   }
00793 
00794   // DgDp
00795   DgDp_default_op_support_.clear();
00796   DgDp_default_mv_support_.clear();
00797   for ( int j = 0; j < l_Ng; ++j ) {
00798     DgDp_default_op_support_.push_back(Array<DefaultDerivLinearOpSupport>());
00799     DgDp_default_mv_support_.push_back(Array<DefaultDerivMvAdjointSupport>());
00800     for ( int l = 0; l < l_Np; ++l ) {
00801       const MEB::DerivativeSupport DgDp_j_l_impl_support =
00802         outArgsImpl.supports(MEB::OUT_ARG_DgDp,j,l);
00803       // LinearOpBase support
00804       const DefaultDerivLinearOpSupport DgDp_j_l_op_support =
00805         determineDefaultDerivLinearOpSupport(DgDp_j_l_impl_support);
00806       DgDp_default_op_support_[j].push_back(DgDp_j_l_op_support);
00807       outArgs.setSupports(
00808         MEB::OUT_ARG_DgDp, j, l,
00809         updateDefaultLinearOpSupport(
00810           DgDp_j_l_impl_support, DgDp_j_l_op_support
00811           )
00812         );
00813       // MultiVectorBase
00814       const DefaultDerivMvAdjointSupport DgDp_j_l_mv_support = 
00815         determineDefaultDerivMvAdjointSupport(
00816           DgDp_j_l_impl_support, *this->get_g_space(j), *this->get_p_space(l)
00817           );
00818       DgDp_default_mv_support_[j].push_back(DgDp_j_l_mv_support);
00819       outArgs.setSupports(
00820         MEB::OUT_ARG_DgDp, j, l,
00821         updateDefaultDerivMvAdjointSupport(
00822           outArgs.supports(MEB::OUT_ARG_DgDp, j, l),
00823           DgDp_j_l_mv_support
00824           )
00825         );
00826     }
00827   }
00828   // 2007/09/09: rabart: ToDo: Move the above code into a private helper
00829   // function!
00830   
00831   // W (given W_op and W_factory)
00832   default_W_support_ = false;
00833   if ( outArgsImpl.supports(MEB::OUT_ARG_W_op) && !is_null(this->get_W_factory())
00834     && !outArgsImpl.supports(MEB::OUT_ARG_W) )
00835   {
00836     default_W_support_ = true;
00837     outArgs.setSupports(MEB::OUT_ARG_W);
00838     outArgs.set_W_properties(outArgsImpl.get_W_properties());
00839   }
00840   
00841   //
00842   // D) All done!
00843   //
00844 
00845   prototypeOutArgs_ = outArgs;
00846   isInitialized_ = true;
00847 
00848 }
00849 
00850 
00851 // Private functions with default implementaton to be overridden by subclasses
00852 
00853 
00854 template<class Scalar>
00855 RCP<LinearOpBase<Scalar> >
00856 ModelEvaluatorDefaultBase<Scalar>::create_DfDp_op_impl(int l) const
00857 {
00858   typedef ModelEvaluatorBase MEB;
00859   MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
00860   TEST_FOR_EXCEPTION(
00861     outArgs.supports(MEB::OUT_ARG_DfDp,l).supports(MEB::DERIV_LINEAR_OP),
00862     std::logic_error,
00863     "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
00864     " supports the LinearOpBase form of DfDp("<<l<<") (as determined from its"
00865     " OutArgs object created by createOutArgsImpl())"
00866     " but this function create_DfDp_op_impl(...) has not been overriden"
00867     " to create such an object!"
00868     );
00869   return Teuchos::null;
00870 }
00871 
00872 
00873 template<class Scalar>
00874 RCP<LinearOpBase<Scalar> >
00875 ModelEvaluatorDefaultBase<Scalar>::create_DgDx_dot_op_impl(int j) const
00876 {
00877   typedef ModelEvaluatorBase MEB;
00878   MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
00879   TEST_FOR_EXCEPTION(
00880     outArgs.supports(MEB::OUT_ARG_DgDx_dot,j).supports(MEB::DERIV_LINEAR_OP),
00881     std::logic_error,
00882     "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
00883     " supports the LinearOpBase form of DgDx_dot("<<j<<") (as determined from"
00884     " its OutArgs object created by createOutArgsImpl())"
00885     " but this function create_DgDx_dot_op_impl(...) has not been overriden"
00886     " to create such an object!"
00887     );
00888   return Teuchos::null;
00889 }
00890 
00891 
00892 template<class Scalar>
00893 RCP<LinearOpBase<Scalar> >
00894 ModelEvaluatorDefaultBase<Scalar>::create_DgDx_op_impl(int j) const
00895 {
00896   typedef ModelEvaluatorBase MEB;
00897   MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
00898   TEST_FOR_EXCEPTION(
00899     outArgs.supports(MEB::OUT_ARG_DgDx,j).supports(MEB::DERIV_LINEAR_OP),
00900     std::logic_error,
00901     "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
00902     " supports the LinearOpBase form of DgDx("<<j<<") (as determined from"
00903     " its OutArgs object created by createOutArgsImpl())"
00904     " but this function create_DgDx_op_impl(...) has not been overriden"
00905     " to create such an object!"
00906     );
00907   return Teuchos::null;
00908 }
00909 
00910 
00911 template<class Scalar>
00912 RCP<LinearOpBase<Scalar> >
00913 ModelEvaluatorDefaultBase<Scalar>::create_DgDp_op_impl(int j, int l) const
00914 {
00915   typedef ModelEvaluatorBase MEB;
00916   MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
00917   TEST_FOR_EXCEPTION(
00918     outArgs.supports(MEB::OUT_ARG_DgDp,j,l).supports(MEB::DERIV_LINEAR_OP),
00919     std::logic_error,
00920     "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
00921     " supports the LinearOpBase form of DgDp("<<j<<","<<l<<")"
00922     " (as determined from its OutArgs object created by createOutArgsImpl())"
00923     " but this function create_DgDp_op_impl(...) has not been overriden"
00924     " to create such an object!"
00925     );
00926   return Teuchos::null;
00927 }
00928 
00929 
00930 // protected
00931 
00932 
00933 template<class Scalar>
00934 ModelEvaluatorDefaultBase<Scalar>::ModelEvaluatorDefaultBase()
00935   :isInitialized_(false), default_W_support_(false)
00936 {}
00937 
00938 
00939 // private
00940 
00941 
00942 template<class Scalar>
00943 void ModelEvaluatorDefaultBase<Scalar>::lazyInitializeDefaultBase() const
00944 {
00945   if (!isInitialized_)
00946     const_cast<ModelEvaluatorDefaultBase<Scalar>*>(this)->initializeDefaultBase();
00947 }
00948 
00949 
00950 template<class Scalar>
00951 void ModelEvaluatorDefaultBase<Scalar>::assert_l(const int l) const
00952 {
00953   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(l,0,this->Np());
00954 }
00955 
00956 
00957 template<class Scalar>
00958 void ModelEvaluatorDefaultBase<Scalar>::assert_j(const int j) const
00959 {
00960   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(j,0,this->Ng());
00961 }
00962 
00963 
00964 // Private static functions
00965 
00966 
00967 template<class Scalar>
00968 ModelEvaluatorDefaultBaseTypes::DefaultDerivLinearOpSupport
00969 ModelEvaluatorDefaultBase<Scalar>::determineDefaultDerivLinearOpSupport(
00970   const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl
00971   )
00972 {
00973   typedef ModelEvaluatorBase MEB;
00974   if (
00975     (
00976       derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
00977       ||
00978       derivSupportImpl.supports(MEB::DERIV_TRANS_MV_BY_ROW)
00979       )
00980     &&
00981     !derivSupportImpl.supports(MEB::DERIV_LINEAR_OP)
00982     )
00983   {
00984     return DefaultDerivLinearOpSupport(
00985       derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
00986       ? MEB::DERIV_MV_BY_COL
00987       : MEB::DERIV_TRANS_MV_BY_ROW
00988       );
00989   }
00990   return DefaultDerivLinearOpSupport();
00991 }
00992 
00993 
00994 template<class Scalar>
00995 RCP<LinearOpBase<Scalar> >
00996 ModelEvaluatorDefaultBase<Scalar>::createDefaultLinearOp(
00997   const DefaultDerivLinearOpSupport &defaultLinearOpSupport,
00998   const RCP<const VectorSpaceBase<Scalar> > &fnc_space,
00999   const RCP<const VectorSpaceBase<Scalar> > &var_space
01000   )
01001 {
01002   using Teuchos::rcp_implicit_cast;
01003   typedef LinearOpBase<Scalar> LOB;
01004   typedef ModelEvaluatorBase MEB;
01005   switch(defaultLinearOpSupport.mvImplOrientation()) {
01006     case MEB::DERIV_MV_BY_COL:
01007       // The MultiVector will do just fine as the LinearOpBase
01008       return createMembers(fnc_space, var_space->dim());
01009     case MEB::DERIV_TRANS_MV_BY_ROW:
01010       // We will have to implicitly transpose the underlying MultiVector
01011       return nonconstAdjoint<Scalar>(
01012         rcp_implicit_cast<LOB>(createMembers(var_space, fnc_space->dim()))
01013         );
01014 #ifdef TEUCHOS_DEBUG
01015     default:
01016       TEST_FOR_EXCEPT(true);
01017 #endif
01018   }
01019   return Teuchos::null; // Will never be called!
01020 }
01021 
01022 
01023 template<class Scalar>
01024 ModelEvaluatorBase::DerivativeSupport
01025 ModelEvaluatorDefaultBase<Scalar>::updateDefaultLinearOpSupport(
01026   const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
01027   const DefaultDerivLinearOpSupport &defaultLinearOpSupport
01028   )
01029 {
01030   typedef ModelEvaluatorBase MEB;
01031   MEB::DerivativeSupport derivSupport = derivSupportImpl;
01032   if (defaultLinearOpSupport.provideDefaultLinearOp())
01033     derivSupport.plus(MEB::DERIV_LINEAR_OP);
01034   return derivSupport;
01035 }
01036 
01037 
01038 template<class Scalar>
01039 ModelEvaluatorBase::Derivative<Scalar>
01040 ModelEvaluatorDefaultBase<Scalar>::getOutArgImplForDefaultLinearOpSupport(
01041   const ModelEvaluatorBase::Derivative<Scalar> &deriv,
01042   const DefaultDerivLinearOpSupport &defaultLinearOpSupport
01043   )
01044 {
01045 
01046   using Teuchos::rcp_dynamic_cast;
01047   typedef ModelEvaluatorBase MEB;
01048   typedef MultiVectorBase<Scalar> MVB;
01049   typedef ScaledAdjointLinearOpBase<Scalar> SALOB;
01050 
01051   // If the derivative is not a LinearOpBase then it it either empty or is a
01052   // MultiVectorBase object, and in either case we just return it since the
01053   // underlying evalModelImpl(...) functions should be able to handle it!
01054   if (is_null(deriv.getLinearOp()))
01055     return deriv;
01056 
01057   // The derivative is LinearOpBase so get out the underlying MultiVectorBase
01058   // object and return its derivative multi-vector form.
01059   switch(defaultLinearOpSupport.mvImplOrientation()) {
01060     case MEB::DERIV_MV_BY_COL: {
01061       return MEB::Derivative<Scalar>(
01062         rcp_dynamic_cast<MVB>(deriv.getLinearOp(),true),
01063         MEB::DERIV_MV_BY_COL
01064         );
01065     }
01066     case MEB::DERIV_TRANS_MV_BY_ROW: {
01067       return MEB::Derivative<Scalar>(
01068         rcp_dynamic_cast<MVB>(
01069           rcp_dynamic_cast<SALOB>(deriv.getLinearOp(),true)->getNonconstOrigOp()
01070           ),
01071         MEB::DERIV_TRANS_MV_BY_ROW
01072         );
01073     }
01074 #ifdef TEUCHOS_DEBUG
01075     default:
01076       TEST_FOR_EXCEPT(true);
01077 #endif
01078   }
01079 
01080   return ModelEvaluatorBase::Derivative<Scalar>(); // Should never get here!
01081 
01082 }
01083 
01084 
01085 template<class Scalar>
01086 ModelEvaluatorDefaultBaseTypes::DefaultDerivMvAdjointSupport
01087 ModelEvaluatorDefaultBase<Scalar>::determineDefaultDerivMvAdjointSupport(
01088   const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
01089   const VectorSpaceBase<Scalar> &fnc_space,
01090   const VectorSpaceBase<Scalar> &var_space
01091   )
01092 {
01093   typedef ModelEvaluatorBase MEB;
01094   // Here we will support the adjoint copy for of a multi-vector if both
01095   // spaces give rise to in-core vectors.
01096   const bool implSupportsMv =
01097     ( derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
01098       || derivSupportImpl.supports(MEB::DERIV_TRANS_MV_BY_ROW) );
01099   const bool bothSpacesHaveInCoreViews =
01100     ( fnc_space.hasInCoreView() && var_space.hasInCoreView() );
01101   if ( implSupportsMv && bothSpacesHaveInCoreViews ) {
01102     return DefaultDerivMvAdjointSupport(
01103       derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
01104       ? MEB::DERIV_TRANS_MV_BY_ROW
01105       : MEB::DERIV_MV_BY_COL
01106       );
01107   }
01108   // We can't provide an adjoint copy or such a copy is not needed!
01109   return DefaultDerivMvAdjointSupport();
01110 }
01111 
01112 
01113 template<class Scalar>
01114 ModelEvaluatorBase::DerivativeSupport
01115 ModelEvaluatorDefaultBase<Scalar>::updateDefaultDerivMvAdjointSupport(
01116   const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
01117   const DefaultDerivMvAdjointSupport &defaultMvAdjointSupport
01118   )
01119 {
01120   typedef ModelEvaluatorBase MEB;
01121   MEB::DerivativeSupport derivSupport = derivSupportImpl;
01122   if (defaultMvAdjointSupport.provideDefaultAdjoint())
01123     derivSupport.plus(defaultMvAdjointSupport.mvAdjointCopyOrientation());
01124   return derivSupport;
01125 }
01126 
01127 
01128 } // namespace Thyra
01129 
01130 
01131 #endif // THYRA_MODEL_EVALUATOR_DEFAULT_BASE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines