Thyra Version of the Day
Thyra_DefaultFiniteDifferenceModelEvaluator_def.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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (bartlettra@ornl.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #ifndef THYRA_DEFAULT_FINITE_DIFFERENCE_MODEL_EVALUATOR_DEF_HPP
00043 #define THYRA_DEFAULT_FINITE_DIFFERENCE_MODEL_EVALUATOR_DEF_HPP
00044 
00045 #include "Thyra_DefaultFiniteDifferenceModelEvaluator_decl.hpp"
00046 
00047 
00048 namespace Thyra {
00049 
00050 
00051 // Constructors/initializers/accessors/utilities
00052 
00053 
00054 template<class Scalar>
00055 DefaultFiniteDifferenceModelEvaluator<Scalar>::DefaultFiniteDifferenceModelEvaluator()
00056 {}
00057 
00058 
00059 template<class Scalar>
00060 void DefaultFiniteDifferenceModelEvaluator<Scalar>::initialize(
00061   const RCP<ModelEvaluator<Scalar> > &thyraModel,
00062   const RCP<DirectionalFiniteDiffCalculator<Scalar> > &direcFiniteDiffCalculator_in
00063   )
00064 {
00065   this->ModelEvaluatorDelegatorBase<Scalar>::initialize(thyraModel);
00066   direcFiniteDiffCalculator_ = direcFiniteDiffCalculator_in;
00067 }
00068 
00069 
00070 // Public functions overridden from Teuchos::Describable
00071 
00072 
00073 template<class Scalar>
00074 std::string DefaultFiniteDifferenceModelEvaluator<Scalar>::description() const
00075 {
00076   const RCP<const ModelEvaluator<Scalar> >
00077     thyraModel = this->getUnderlyingModel();
00078   std::ostringstream oss;
00079   oss << "Thyra::DefaultFiniteDifferenceModelEvaluator{";
00080   oss << "thyraModel=";
00081   if(thyraModel.get())
00082     oss << "\'"<<thyraModel->description()<<"\'";
00083   else
00084     oss << "NULL";
00085   oss << "}";
00086   return oss.str();
00087 }
00088 
00089 
00090 // Private functions overridden from ModelEvaulatorDefaultBase
00091 
00092 
00093 template<class Scalar>
00094 ModelEvaluatorBase::OutArgs<Scalar>
00095 DefaultFiniteDifferenceModelEvaluator<Scalar>::createOutArgsImpl() const
00096 {
00097   typedef ModelEvaluatorBase MEB;
00098   typedef MEB::DerivativeSupport DS;
00099   const RCP<const ModelEvaluator<Scalar> >
00100     thyraModel = this->getUnderlyingModel();
00101   const MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
00102   const int l_Np = wrappedOutArgs.Np(), l_Ng = wrappedOutArgs.Ng();
00103   MEB::OutArgsSetup<Scalar> outArgs;
00104   outArgs.setModelEvalDescription(this->description());
00105   outArgs.set_Np_Ng(l_Np,l_Ng);
00106   outArgs.setSupports(wrappedOutArgs);
00107   if (wrappedOutArgs.supports(MEB::OUT_ARG_f)) {
00108     for( int l = 0; l < l_Np; ++l ) {
00109       outArgs.setSupports(MEB::OUT_ARG_DfDp,l,MEB::DERIV_MV_BY_COL);
00110     }
00111   }
00112   for( int j = 0; j < l_Ng; ++j ) {
00113     for( int l = 0; l < l_Np; ++l ) {
00114       outArgs.setSupports( MEB::OUT_ARG_DgDp , j, l, MEB::DERIV_MV_BY_COL);
00115     }
00116   }
00117   // ToDo: Add support for more derivatives as needed!
00118   return outArgs;
00119 }
00120 
00121 
00122 template<class Scalar>
00123 void DefaultFiniteDifferenceModelEvaluator<Scalar>::evalModelImpl(
00124   const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00125   const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00126   ) const
00127 {
00128   using Teuchos::rcp;
00129   using Teuchos::rcp_const_cast;
00130   using Teuchos::rcp_dynamic_cast;
00131   using Teuchos::OSTab;
00132   typedef Teuchos::ScalarTraits<Scalar> ST;
00133   //typedef typename ST::magnitudeType ScalarMag;
00134   typedef ModelEvaluatorBase MEB;
00135   namespace DFDCT = DirectionalFiniteDiffCalculatorTypes;
00136 
00137   typedef RCP<VectorBase<Scalar> > V_ptr;
00138   typedef RCP<const VectorBase<Scalar> > CV_ptr;
00139   typedef RCP<MultiVectorBase<Scalar> > MV_ptr;
00140 
00141   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_BEGIN(
00142     "Thyra::DefaultFiniteDifferenceModelEvaluator",inArgs,outArgs
00143     );
00144 
00145   //
00146   // Note: Just do derivatives DfDp(l) and DgDp(j,l) for now!
00147   //
00148 
00149   const RCP<const VectorSpaceBase<Scalar> >
00150     p_space = thyraModel->get_p_space(0),
00151     g_space = thyraModel->get_g_space(0);
00152 
00153   //
00154   // A) Compute the base point
00155   //
00156 
00157   if(out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00158     *out << "\nComputing the base point ...\n";
00159 
00160   const int l_Np = outArgs.Np();
00161   const int l_Ng = outArgs.Ng();
00162   MEB::InArgs<Scalar> wrappedInArgs = inArgs;
00163   MEB::OutArgs<Scalar> baseFunc = thyraModel->createOutArgs();
00164   if( outArgs.supports(MEB::OUT_ARG_f) && outArgs.get_f().get() )
00165     baseFunc.set_f(outArgs.get_f());
00166   for( int j = 0; j < l_Ng; ++j ) {
00167     V_ptr g_j;
00168     if( (g_j=outArgs.get_g(j)).get() )
00169       baseFunc.set_g(j,g_j);
00170   }
00171   // 2007/08/27: We really should really try to allow some derivatives to pass
00172   // through and some derivatives to be computed by finite differences. Right
00173   // now, if you use this class, all derivatives w.r.t. parameters are finite
00174   // differenced and that is not given the user enough control!
00175 
00176   thyraModel->evalModel(wrappedInArgs,baseFunc);
00177 
00178   bool failed = baseFunc.isFailed();
00179 
00180   //
00181   // B) Compute the derivatives
00182   //
00183  
00184   if(!failed) {
00185 
00186     // a) Determine what derivatives you need to support first
00187 
00188     Array<int> compute_DfDp;
00189     Array<Array<int> > compute_DgDp(l_Ng);
00190     DFDCT::SelectedDerivatives selectedDerivs;
00191 
00192     for ( int l = 0; l < l_Np; ++l ) {
00193 
00194       // DfDp(l)
00195       if(
00196         outArgs.supports(MEB::OUT_ARG_DfDp,l).none()==false
00197         &&
00198         outArgs.get_DfDp(l).isEmpty()==false
00199         )
00200       {
00201         selectedDerivs.supports(MEB::OUT_ARG_DfDp,l);
00202         compute_DfDp.push_back(true);
00203       }
00204       else
00205       {
00206         compute_DfDp.push_back(false);
00207       }
00208 
00209       // DgDp(j=0...,l)
00210       for ( int j = 0; j < l_Ng; ++j ) {
00211         if(
00212           outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none()==false
00213           &&
00214           outArgs.get_DgDp(j,l).isEmpty()==false
00215           )
00216         {
00217           selectedDerivs.supports(MEB::OUT_ARG_DgDp,j,l);
00218           compute_DgDp[j].push_back(true);
00219         }
00220         else
00221         {
00222           compute_DgDp[j].push_back(false);
00223         }
00224       }
00225     }
00226 
00227     // b) Create the deriv OutArgs and set the output objects that need to be
00228     // computed with finite differences
00229  
00230     MEB::OutArgs<Scalar>
00231       deriv = direcFiniteDiffCalculator_->createOutArgs(
00232         *thyraModel, selectedDerivs );
00233 
00234     for ( int l = 0; l < l_Np; ++l ) {
00235       if ( compute_DfDp[l] )
00236         deriv.set_DfDp(l,outArgs.get_DfDp(l));
00237       for ( int j = 0; j < l_Ng; ++j ) {
00238         if ( compute_DgDp[j][l] )
00239           deriv.set_DgDp(j,l,outArgs.get_DgDp(j,l));
00240       }
00241     }
00242 
00243     // c) Compute the missing functions with finite differences!
00244 
00245     direcFiniteDiffCalculator_->calcDerivatives(
00246       *thyraModel,inArgs,baseFunc,deriv
00247       );
00248 
00249   }
00250 
00251   if(failed) {
00252     if(out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00253       *out << "\nEvaluation failed, returning NaNs ...\n";
00254     outArgs.setFailed();
00255   }
00256 
00257   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00258  
00259 }
00260 
00261 
00262 } // namespace Thyra
00263 
00264 
00265 #endif // THYRA_DEFAULT_FINITE_DIFFERENCE_MODEL_EVALUATOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines