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

Generated on Tue Oct 20 12:47:13 2009 for Thyra Nonlinear Model Evaluator Support by doxygen 1.4.7