Thyra_ModelEvaluatorDefaultBase.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 // 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   typedef ModelEvaluatorBase MEB;
00494 
00495   lazyInitializeDefaultBase();
00496 
00497   const int l_Np = outArgs.Np();
00498   const int l_Ng = outArgs.Ng();
00499 
00500   //
00501   // A) Assert that the inArgs and outArgs object match this class!
00502   //
00503 
00504 #ifdef TEUCHOS_DEBUG
00505   assertInArgsEvalObjects(*this,inArgs);
00506   assertOutArgsEvalObjects(*this,outArgs,&inArgs);
00507 #endif  
00508   
00509   //
00510   // B) Setup the OutArgs object for the underlying implementation's
00511   // evalModelImpl(...) function
00512   //
00513 
00514   MEB::OutArgs<Scalar> outArgsImpl = this->createOutArgsImpl();
00515   Array<MultiVectorAdjointPair> DgDp_temp_adjoint_copies;
00516 
00517   {
00518 
00519     outArgsImpl.setArgs(outArgs,true);
00520 
00521     // DfDp(l)
00522     if (outArgsImpl.supports(MEB::OUT_ARG_f)) {
00523       for ( int l = 0; l < l_Np; ++l ) {
00524         const DefaultDerivLinearOpSupport defaultLinearOpSupport =
00525           DfDp_default_op_support_[l];
00526         if (defaultLinearOpSupport.provideDefaultLinearOp()) {
00527           outArgsImpl.set_DfDp( l,
00528             getOutArgImplForDefaultLinearOpSupport(
00529               outArgs.get_DfDp(l), defaultLinearOpSupport
00530               )
00531             );
00532         }
00533         else {
00534           // DfDp(l) already set by outArgsImpl.setArgs(...)!
00535         }
00536       }
00537     }
00538 
00539     // DgDx_dot(j)
00540     for ( int j = 0; j < l_Ng; ++j ) {
00541       const DefaultDerivLinearOpSupport defaultLinearOpSupport =
00542         DgDx_dot_default_op_support_[j];
00543       if (defaultLinearOpSupport.provideDefaultLinearOp()) {
00544         outArgsImpl.set_DgDx_dot( j,
00545           getOutArgImplForDefaultLinearOpSupport(
00546             outArgs.get_DgDx_dot(j), defaultLinearOpSupport
00547             )
00548           );
00549       }
00550       else {
00551         // DgDx_dot(j) already set by outArgsImpl.setArgs(...)!
00552       }
00553     }
00554 
00555     // DgDx(j)
00556     for ( int j = 0; j < l_Ng; ++j ) {
00557       const DefaultDerivLinearOpSupport defaultLinearOpSupport =
00558         DgDx_default_op_support_[j];
00559       if (defaultLinearOpSupport.provideDefaultLinearOp()) {
00560         outArgsImpl.set_DgDx( j,
00561           getOutArgImplForDefaultLinearOpSupport(
00562             outArgs.get_DgDx(j), defaultLinearOpSupport
00563             )
00564           );
00565       }
00566       else {
00567         // DgDx(j) already set by outArgsImpl.setArgs(...)!
00568       }
00569     }
00570 
00571     // DgDp(j,l)
00572     for ( int j = 0; j < l_Ng; ++j ) {
00573       const Array<DefaultDerivLinearOpSupport> &DgDp_default_op_support_j =
00574         DgDp_default_op_support_[j];
00575       const Array<DefaultDerivMvAdjointSupport> &DgDp_default_mv_support_j =
00576         DgDp_default_mv_support_[j];
00577       for ( int l = 0; l < l_Np; ++l ) {
00578         const DefaultDerivLinearOpSupport defaultLinearOpSupport =
00579           DgDp_default_op_support_j[l];
00580         const DefaultDerivMvAdjointSupport defaultMvAdjointSupport =
00581           DgDp_default_mv_support_j[l];
00582         MEB::Derivative<Scalar> DgDp_j_l;
00583         if (!outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none())
00584           DgDp_j_l = outArgs.get_DgDp(j,l);
00585         if (
00586           defaultLinearOpSupport.provideDefaultLinearOp()
00587           && !is_null(DgDp_j_l.getLinearOp())
00588           )
00589         {
00590           outArgsImpl.set_DgDp( j, l,
00591             getOutArgImplForDefaultLinearOpSupport(
00592               DgDp_j_l, defaultLinearOpSupport
00593               )
00594             );
00595         }
00596         else if (
00597           defaultMvAdjointSupport.provideDefaultAdjoint()
00598           && !is_null(DgDp_j_l.getMultiVector())
00599           )
00600         {
00601           const RCP<MultiVectorBase<Scalar> > DgDp_j_l_mv = 
00602             DgDp_j_l.getMultiVector();
00603           if (
00604             defaultMvAdjointSupport.mvAdjointCopyOrientation()
00605             ==
00606             DgDp_j_l.getMultiVectorOrientation()
00607             )
00608           {
00609             // The orientation of the multi-vector is different so we need to
00610             // create a temporary copy to pass to evalModelImpl(...) and then
00611             // copy it back again!
00612             const RCP<MultiVectorBase<Scalar> > DgDp_j_l_mv_adj =
00613               createMembers(DgDp_j_l_mv->domain(), DgDp_j_l_mv->range()->dim());
00614             outArgsImpl.set_DgDp( j, l,
00615               MEB::Derivative<Scalar>(
00616                 DgDp_j_l_mv_adj,
00617                 getOtherDerivativeMultiVectorOrientation(
00618                   defaultMvAdjointSupport.mvAdjointCopyOrientation()
00619                   )
00620                 )
00621               );
00622             // Remember these multi-vectors so that we can do the transpose copy
00623             // back after the evaluation!
00624             DgDp_temp_adjoint_copies.push_back(
00625               MultiVectorAdjointPair(DgDp_j_l_mv, DgDp_j_l_mv_adj)
00626               );
00627           }
00628           else {
00629             // The form of the multi-vector is supported by evalModelImpl(..)
00630             // and is already set on the outArgsImpl object.
00631           }
00632         }
00633         else {
00634           // DgDp(j,l) already set by outArgsImpl.setArgs(...)!
00635         }
00636       }
00637     }
00638 
00639     // W
00640     {
00641       RCP<LinearOpWithSolveBase<Scalar> > W;
00642       if ( default_W_support_ && !is_null(W=outArgs.get_W()) ) {
00643         const RCP<const LinearOpWithSolveFactoryBase<Scalar> >
00644           W_factory = this->get_W_factory();
00645         // Extract the underlying W_op object (if it exists)
00646         RCP<const LinearOpBase<Scalar> > W_op_const;
00647         uninitializeOp<Scalar>( *W_factory, &*W, &W_op_const );
00648         RCP<LinearOpBase<Scalar> > W_op;
00649         if (!is_null(W_op_const)) {
00650           // Here we remove the const.  This is perfectly safe since
00651           // w.r.t. this class, we put a non-const object in there and we can
00652           // expect to change that object after the fact.  That is our
00653           // prerogative.
00654           W_op = Teuchos::rcp_const_cast<LinearOpBase<Scalar> >(W_op_const);
00655         }
00656         else {
00657           // The W_op object has not been initialized yet so create it.  The
00658           // next time through, we should not have to do this!
00659           W_op = this->create_W_op();
00660         }
00661         outArgsImpl.set_W_op(W_op);
00662       }
00663     }
00664     
00665   }
00666 
00667   //
00668   // C) Evaluate the underlying model implementation!
00669   //
00670 
00671   this->evalModelImpl( inArgs, outArgsImpl );
00672 
00673   //
00674   // D) Post-process the output arguments
00675   //
00676 
00677   // Do explicit transposes for DgDp(j,l) if needed
00678   const int numMvAdjointCopies = DgDp_temp_adjoint_copies.size();
00679   for ( int adj_copy_i = 0; adj_copy_i < numMvAdjointCopies; ++adj_copy_i ) {
00680     const MultiVectorAdjointPair adjPair =
00681       DgDp_temp_adjoint_copies[adj_copy_i];
00682     doExplicitMultiVectorAdjoint( *adjPair.mvImplAdjoint, &*adjPair.mvOuter );
00683   }
00684   
00685   // Update W given W_op and W_factory
00686   {
00687     RCP<LinearOpWithSolveBase<Scalar> > W;
00688     if ( default_W_support_ && !is_null(W=outArgs.get_W()) ) {
00689       const RCP<const LinearOpWithSolveFactoryBase<Scalar> >
00690         W_factory = this->get_W_factory();
00691       W_factory->setOStream(this->getOStream());
00692       W_factory->setVerbLevel(this->getVerbLevel());
00693       initializeOp<Scalar>(*W_factory, outArgsImpl.get_W_op().getConst(), &*W );
00694     }
00695   }
00696   
00697 }
00698 
00699 
00700 // protected
00701 
00702 
00703 // Setup functions called by subclasses
00704 
00705 template<class Scalar>
00706 void ModelEvaluatorDefaultBase<Scalar>::initializeDefaultBase()
00707 {
00708 
00709   typedef ModelEvaluatorBase MEB;
00710 
00711   // In case we throw half way thorugh, set to uninitialized
00712   isInitialized_ = false;
00713   default_W_support_ = false;
00714 
00715   //
00716   // A) Get the InArgs and OutArgs from the subclass
00717   //
00718   
00719   const MEB::InArgs<Scalar> inArgs = this->createInArgs();
00720   const MEB::OutArgs<Scalar> outArgsImpl = this->createOutArgsImpl();
00721 
00722   //
00723   // B) Validate the subclasses InArgs and OutArgs
00724   //
00725 
00726 #ifdef TEUCHOS_DEBUG
00727   assertInArgsOutArgsSetup( this->description(), inArgs, outArgsImpl );
00728 #endif // TEUCHOS_DEBUG
00729 
00730   //
00731   // C) Set up support for default derivative objects and prototype OutArgs
00732   //
00733 
00734   const int l_Ng = outArgsImpl.Ng();
00735   const int l_Np = outArgsImpl.Np();
00736 
00737   // Set support for all outputs supported in the underly implementation
00738   MEB::OutArgsSetup<Scalar> outArgs;
00739   outArgs.setModelEvalDescription(this->description());
00740   outArgs.set_Np_Ng(l_Np,l_Ng);
00741   outArgs.setSupports(outArgsImpl);
00742 
00743   // DfDp
00744   DfDp_default_op_support_.clear();
00745   if (outArgs.supports(MEB::OUT_ARG_f)) {
00746     for ( int l = 0; l < l_Np; ++l ) {
00747       const MEB::DerivativeSupport DfDp_l_impl_support =
00748         outArgsImpl.supports(MEB::OUT_ARG_DfDp,l);
00749       const DefaultDerivLinearOpSupport DfDp_l_op_support =
00750         determineDefaultDerivLinearOpSupport(DfDp_l_impl_support);
00751       DfDp_default_op_support_.push_back(DfDp_l_op_support);
00752       outArgs.setSupports(
00753         MEB::OUT_ARG_DfDp, l,
00754         updateDefaultLinearOpSupport(
00755           DfDp_l_impl_support, DfDp_l_op_support
00756           )
00757         );
00758     }
00759   }
00760 
00761   // DgDx_dot
00762   DgDx_dot_default_op_support_.clear();
00763   for ( int j = 0; j < l_Ng; ++j ) {
00764     const MEB::DerivativeSupport DgDx_dot_j_impl_support =
00765       outArgsImpl.supports(MEB::OUT_ARG_DgDx_dot,j);
00766     const DefaultDerivLinearOpSupport DgDx_dot_j_op_support =
00767       determineDefaultDerivLinearOpSupport(DgDx_dot_j_impl_support);
00768     DgDx_dot_default_op_support_.push_back(DgDx_dot_j_op_support);
00769     outArgs.setSupports(
00770       MEB::OUT_ARG_DgDx_dot, j,
00771       updateDefaultLinearOpSupport(
00772         DgDx_dot_j_impl_support, DgDx_dot_j_op_support
00773         )
00774       );
00775   }
00776 
00777   // DgDx
00778   DgDx_default_op_support_.clear();
00779   for ( int j = 0; j < l_Ng; ++j ) {
00780     const MEB::DerivativeSupport DgDx_j_impl_support =
00781       outArgsImpl.supports(MEB::OUT_ARG_DgDx,j);
00782     const DefaultDerivLinearOpSupport DgDx_j_op_support =
00783       determineDefaultDerivLinearOpSupport(DgDx_j_impl_support);
00784     DgDx_default_op_support_.push_back(DgDx_j_op_support);
00785     outArgs.setSupports(
00786       MEB::OUT_ARG_DgDx, j,
00787       updateDefaultLinearOpSupport(
00788         DgDx_j_impl_support, DgDx_j_op_support
00789         )
00790       );
00791   }
00792 
00793   // DgDp
00794   DgDp_default_op_support_.clear();
00795   DgDp_default_mv_support_.clear();
00796   for ( int j = 0; j < l_Ng; ++j ) {
00797     DgDp_default_op_support_.push_back(Array<DefaultDerivLinearOpSupport>());
00798     DgDp_default_mv_support_.push_back(Array<DefaultDerivMvAdjointSupport>());
00799     for ( int l = 0; l < l_Np; ++l ) {
00800       const MEB::DerivativeSupport DgDp_j_l_impl_support =
00801         outArgsImpl.supports(MEB::OUT_ARG_DgDp,j,l);
00802       // LinearOpBase support
00803       const DefaultDerivLinearOpSupport DgDp_j_l_op_support =
00804         determineDefaultDerivLinearOpSupport(DgDp_j_l_impl_support);
00805       DgDp_default_op_support_[j].push_back(DgDp_j_l_op_support);
00806       outArgs.setSupports(
00807         MEB::OUT_ARG_DgDp, j, l,
00808         updateDefaultLinearOpSupport(
00809           DgDp_j_l_impl_support, DgDp_j_l_op_support
00810           )
00811         );
00812       // MultiVectorBase
00813       const DefaultDerivMvAdjointSupport DgDp_j_l_mv_support = 
00814         determineDefaultDerivMvAdjointSupport(
00815           DgDp_j_l_impl_support, *this->get_g_space(j), *this->get_p_space(l)
00816           );
00817       DgDp_default_mv_support_[j].push_back(DgDp_j_l_mv_support);
00818       outArgs.setSupports(
00819         MEB::OUT_ARG_DgDp, j, l,
00820         updateDefaultDerivMvAdjointSupport(
00821           outArgs.supports(MEB::OUT_ARG_DgDp, j, l),
00822           DgDp_j_l_mv_support
00823           )
00824         );
00825     }
00826   }
00827   // 2007/09/09: rabart: ToDo: Move the above code into a private helper
00828   // function!
00829   
00830   // W (given W_op and W_factory)
00831   default_W_support_ = false;
00832   if ( outArgsImpl.supports(MEB::OUT_ARG_W_op) && !is_null(this->get_W_factory())
00833     && !outArgsImpl.supports(MEB::OUT_ARG_W) )
00834   {
00835     default_W_support_ = true;
00836     outArgs.setSupports(MEB::OUT_ARG_W);
00837     outArgs.set_W_properties(outArgsImpl.get_W_properties());
00838   }
00839   
00840   //
00841   // D) All done!
00842   //
00843 
00844   prototypeOutArgs_ = outArgs;
00845   isInitialized_ = true;
00846 
00847 }
00848 
00849 
00850 // Private functions with default implementaton to be overridden by subclasses
00851 
00852 
00853 template<class Scalar>
00854 RCP<LinearOpBase<Scalar> >
00855 ModelEvaluatorDefaultBase<Scalar>::create_DfDp_op_impl(int l) const
00856 {
00857   typedef ModelEvaluatorBase MEB;
00858   MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
00859   TEST_FOR_EXCEPTION(
00860     outArgs.supports(MEB::OUT_ARG_DfDp,l).supports(MEB::DERIV_LINEAR_OP),
00861     std::logic_error,
00862     "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
00863     " supports the LinearOpBase form of DfDp("<<l<<") (as determined from its"
00864     " OutArgs object created by createOutArgsImpl())"
00865     " but this function create_DfDp_op_impl(...) has not been overriden"
00866     " to create such an object!"
00867     );
00868   return Teuchos::null;
00869 }
00870 
00871 
00872 template<class Scalar>
00873 RCP<LinearOpBase<Scalar> >
00874 ModelEvaluatorDefaultBase<Scalar>::create_DgDx_dot_op_impl(int j) const
00875 {
00876   typedef ModelEvaluatorBase MEB;
00877   MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
00878   TEST_FOR_EXCEPTION(
00879     outArgs.supports(MEB::OUT_ARG_DgDx_dot,j).supports(MEB::DERIV_LINEAR_OP),
00880     std::logic_error,
00881     "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
00882     " supports the LinearOpBase form of DgDx_dot("<<j<<") (as determined from"
00883     " its OutArgs object created by createOutArgsImpl())"
00884     " but this function create_DgDx_dot_op_impl(...) has not been overriden"
00885     " to create such an object!"
00886     );
00887   return Teuchos::null;
00888 }
00889 
00890 
00891 template<class Scalar>
00892 RCP<LinearOpBase<Scalar> >
00893 ModelEvaluatorDefaultBase<Scalar>::create_DgDx_op_impl(int j) const
00894 {
00895   typedef ModelEvaluatorBase MEB;
00896   MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
00897   TEST_FOR_EXCEPTION(
00898     outArgs.supports(MEB::OUT_ARG_DgDx,j).supports(MEB::DERIV_LINEAR_OP),
00899     std::logic_error,
00900     "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
00901     " supports the LinearOpBase form of DgDx("<<j<<") (as determined from"
00902     " its OutArgs object created by createOutArgsImpl())"
00903     " but this function create_DgDx_op_impl(...) has not been overriden"
00904     " to create such an object!"
00905     );
00906   return Teuchos::null;
00907 }
00908 
00909 
00910 template<class Scalar>
00911 RCP<LinearOpBase<Scalar> >
00912 ModelEvaluatorDefaultBase<Scalar>::create_DgDp_op_impl(int j, int l) const
00913 {
00914   typedef ModelEvaluatorBase MEB;
00915   MEB::OutArgs<Scalar> outArgs = this->createOutArgsImpl();
00916   TEST_FOR_EXCEPTION(
00917     outArgs.supports(MEB::OUT_ARG_DgDp,j,l).supports(MEB::DERIV_LINEAR_OP),
00918     std::logic_error,
00919     "Error, The ModelEvaluator subclass "<<this->description()<<" says that it"
00920     " supports the LinearOpBase form of DgDp("<<j<<","<<l<<")"
00921     " (as determined from its OutArgs object created by createOutArgsImpl())"
00922     " but this function create_DgDp_op_impl(...) has not been overriden"
00923     " to create such an object!"
00924     );
00925   return Teuchos::null;
00926 }
00927 
00928 
00929 // protected
00930 
00931 
00932 template<class Scalar>
00933 ModelEvaluatorDefaultBase<Scalar>::ModelEvaluatorDefaultBase()
00934   :isInitialized_(false), default_W_support_(false)
00935 {}
00936 
00937 
00938 // private
00939 
00940 
00941 template<class Scalar>
00942 void ModelEvaluatorDefaultBase<Scalar>::lazyInitializeDefaultBase() const
00943 {
00944   if (!isInitialized_)
00945     const_cast<ModelEvaluatorDefaultBase<Scalar>*>(this)->initializeDefaultBase();
00946 }
00947 
00948 
00949 template<class Scalar>
00950 void ModelEvaluatorDefaultBase<Scalar>::assert_l(const int l) const
00951 {
00952   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(l,0,this->Np());
00953 }
00954 
00955 
00956 template<class Scalar>
00957 void ModelEvaluatorDefaultBase<Scalar>::assert_j(const int j) const
00958 {
00959   TEUCHOS_ASSERT_IN_RANGE_UPPER_EXCLUSIVE(j,0,this->Ng());
00960 }
00961 
00962 
00963 // Private static functions
00964 
00965 
00966 template<class Scalar>
00967 ModelEvaluatorDefaultBaseTypes::DefaultDerivLinearOpSupport
00968 ModelEvaluatorDefaultBase<Scalar>::determineDefaultDerivLinearOpSupport(
00969   const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl
00970   )
00971 {
00972   typedef ModelEvaluatorBase MEB;
00973   if (
00974     (
00975       derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
00976       ||
00977       derivSupportImpl.supports(MEB::DERIV_TRANS_MV_BY_ROW)
00978       )
00979     &&
00980     !derivSupportImpl.supports(MEB::DERIV_LINEAR_OP)
00981     )
00982   {
00983     return DefaultDerivLinearOpSupport(
00984       derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
00985       ? MEB::DERIV_MV_BY_COL
00986       : MEB::DERIV_TRANS_MV_BY_ROW
00987       );
00988   }
00989   return DefaultDerivLinearOpSupport();
00990 }
00991 
00992 
00993 template<class Scalar>
00994 RCP<LinearOpBase<Scalar> >
00995 ModelEvaluatorDefaultBase<Scalar>::createDefaultLinearOp(
00996   const DefaultDerivLinearOpSupport &defaultLinearOpSupport,
00997   const RCP<const VectorSpaceBase<Scalar> > &fnc_space,
00998   const RCP<const VectorSpaceBase<Scalar> > &var_space
00999   )
01000 {
01001   using Teuchos::rcp_implicit_cast;
01002   typedef LinearOpBase<Scalar> LOB;
01003   typedef ModelEvaluatorBase MEB;
01004   switch(defaultLinearOpSupport.mvImplOrientation()) {
01005     case MEB::DERIV_MV_BY_COL:
01006       // The MultiVector will do just fine as the LinearOpBase
01007       return createMembers(fnc_space, var_space->dim());
01008     case MEB::DERIV_TRANS_MV_BY_ROW:
01009       // We will have to implicitly transpose the underlying MultiVector
01010       return nonconstAdjoint<Scalar>(
01011         rcp_implicit_cast<LOB>(createMembers(var_space, fnc_space->dim()))
01012         );
01013 #ifdef TEUCHOS_DEBUG
01014     default:
01015       TEST_FOR_EXCEPT(true);
01016 #endif
01017   }
01018   return Teuchos::null; // Will never be called!
01019 }
01020 
01021 
01022 template<class Scalar>
01023 ModelEvaluatorBase::DerivativeSupport
01024 ModelEvaluatorDefaultBase<Scalar>::updateDefaultLinearOpSupport(
01025   const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
01026   const DefaultDerivLinearOpSupport &defaultLinearOpSupport
01027   )
01028 {
01029   typedef ModelEvaluatorBase MEB;
01030   MEB::DerivativeSupport derivSupport = derivSupportImpl;
01031   if (defaultLinearOpSupport.provideDefaultLinearOp())
01032     derivSupport.plus(MEB::DERIV_LINEAR_OP);
01033   return derivSupport;
01034 }
01035 
01036 
01037 template<class Scalar>
01038 ModelEvaluatorBase::Derivative<Scalar>
01039 ModelEvaluatorDefaultBase<Scalar>::getOutArgImplForDefaultLinearOpSupport(
01040   const ModelEvaluatorBase::Derivative<Scalar> &deriv,
01041   const DefaultDerivLinearOpSupport &defaultLinearOpSupport
01042   )
01043 {
01044 
01045   using Teuchos::rcp_dynamic_cast;
01046   typedef ModelEvaluatorBase MEB;
01047   typedef MultiVectorBase<Scalar> MVB;
01048   typedef ScaledAdjointLinearOpBase<Scalar> SALOB;
01049 
01050   // If the derivative is not a LinearOpBase then it it either empty or is a
01051   // MultiVectorBase object, and in either case we just return it since the
01052   // underlying evalModelImpl(...) functions should be able to handle it!
01053   if (is_null(deriv.getLinearOp()))
01054     return deriv;
01055 
01056   // The derivative is LinearOpBase so get out the underlying MultiVectorBase
01057   // object and return its derivative multi-vector form.
01058   switch(defaultLinearOpSupport.mvImplOrientation()) {
01059     case MEB::DERIV_MV_BY_COL: {
01060       return MEB::Derivative<Scalar>(
01061         rcp_dynamic_cast<MVB>(deriv.getLinearOp(),true),
01062         MEB::DERIV_MV_BY_COL
01063         );
01064     }
01065     case MEB::DERIV_TRANS_MV_BY_ROW: {
01066       return MEB::Derivative<Scalar>(
01067         rcp_dynamic_cast<MVB>(
01068           rcp_dynamic_cast<SALOB>(deriv.getLinearOp(),true)->getNonconstOrigOp()
01069           ),
01070         MEB::DERIV_TRANS_MV_BY_ROW
01071         );
01072     }
01073 #ifdef TEUCHOS_DEBUG
01074     default:
01075       TEST_FOR_EXCEPT(true);
01076 #endif
01077   }
01078 
01079   return ModelEvaluatorBase::Derivative<Scalar>(); // Should never get here!
01080 
01081 }
01082 
01083 
01084 template<class Scalar>
01085 ModelEvaluatorDefaultBaseTypes::DefaultDerivMvAdjointSupport
01086 ModelEvaluatorDefaultBase<Scalar>::determineDefaultDerivMvAdjointSupport(
01087   const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
01088   const VectorSpaceBase<Scalar> &fnc_space,
01089   const VectorSpaceBase<Scalar> &var_space
01090   )
01091 {
01092   typedef ModelEvaluatorBase MEB;
01093   // Here we will support the adjoint copy for of a multi-vector if both
01094   // spaces give rise to in-core vectors.
01095   const bool implSupportsMv =
01096     ( derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
01097       || derivSupportImpl.supports(MEB::DERIV_TRANS_MV_BY_ROW) );
01098   const bool bothSpacesHaveInCoreViews =
01099     ( fnc_space.hasInCoreView() && var_space.hasInCoreView() );
01100   if ( implSupportsMv && bothSpacesHaveInCoreViews ) {
01101     return DefaultDerivMvAdjointSupport(
01102       derivSupportImpl.supports(MEB::DERIV_MV_BY_COL)
01103       ? MEB::DERIV_TRANS_MV_BY_ROW
01104       : MEB::DERIV_MV_BY_COL
01105       );
01106   }
01107   // We can't provide an adjoint copy or such a copy is not needed!
01108   return DefaultDerivMvAdjointSupport();
01109 }
01110 
01111 
01112 template<class Scalar>
01113 ModelEvaluatorBase::DerivativeSupport
01114 ModelEvaluatorDefaultBase<Scalar>::updateDefaultDerivMvAdjointSupport(
01115   const ModelEvaluatorBase::DerivativeSupport &derivSupportImpl,
01116   const DefaultDerivMvAdjointSupport &defaultMvAdjointSupport
01117   )
01118 {
01119   typedef ModelEvaluatorBase MEB;
01120   MEB::DerivativeSupport derivSupport = derivSupportImpl;
01121   if (defaultMvAdjointSupport.provideDefaultAdjoint())
01122     derivSupport.plus(defaultMvAdjointSupport.mvAdjointCopyOrientation());
01123   return derivSupport;
01124 }
01125 
01126 
01127 } // namespace Thyra
01128 
01129 
01130 #endif // THYRA_MODEL_EVALUATOR_DEFAULT_BASE_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines
Generated on Wed Apr 13 10:00:17 2011 for Thyra by  doxygen 1.6.3