Thyra Package Browser (Single Doxygen Collection) Version of the Day
Thyra_DefaultFiniteDifferenceModelEvaluator_def.hpp
Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 //    Thyra: Interfaces and Support for Abstract Numerical Algorithms
00005 //                 Copyright (2004) Sandia Corporation
00006 // 
00007 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
00008 // license for use of this work by or on behalf of the U.S. Government.
00009 // 
00010 // This library is free software; you can redistribute it and/or modify
00011 // it under the terms of the GNU Lesser General Public License as
00012 // published by the Free Software Foundation; either version 2.1 of the
00013 // License, or (at your option) any later version.
00014 //  
00015 // This library is distributed in the hope that it will be useful, but
00016 // WITHOUT ANY WARRANTY; without even the implied warranty of
00017 // MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
00018 // Lesser General Public License for more details.
00019 //  
00020 // You should have received a copy of the GNU Lesser General Public
00021 // License along with this library; if not, write to the Free Software
00022 // Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307
00023 // USA
00024 // Questions? Contact Michael A. Heroux (maherou@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #ifndef THYRA_DEFAULT_FINITE_DIFFERENCE_MODEL_EVALUATOR_DEF_HPP
00030 #define THYRA_DEFAULT_FINITE_DIFFERENCE_MODEL_EVALUATOR_DEF_HPP
00031 
00032 #include "Thyra_DefaultFiniteDifferenceModelEvaluator_decl.hpp"
00033 
00034 
00035 namespace Thyra {
00036 
00037 
00038 // Constructors/initializers/accessors/utilities
00039 
00040 
00041 template<class Scalar>
00042 DefaultFiniteDifferenceModelEvaluator<Scalar>::DefaultFiniteDifferenceModelEvaluator()
00043 {}
00044 
00045 
00046 template<class Scalar>
00047 void DefaultFiniteDifferenceModelEvaluator<Scalar>::initialize(
00048   const RCP<ModelEvaluator<Scalar> > &thyraModel,
00049   const RCP<DirectionalFiniteDiffCalculator<Scalar> > &direcFiniteDiffCalculator_in
00050   )
00051 {
00052   this->ModelEvaluatorDelegatorBase<Scalar>::initialize(thyraModel);
00053   direcFiniteDiffCalculator_ = direcFiniteDiffCalculator_in;
00054 }
00055 
00056 
00057 // Public functions overridden from Teuchos::Describable
00058 
00059 
00060 template<class Scalar>
00061 std::string DefaultFiniteDifferenceModelEvaluator<Scalar>::description() const
00062 {
00063   const RCP<const ModelEvaluator<Scalar> >
00064     thyraModel = this->getUnderlyingModel();
00065   std::ostringstream oss;
00066   oss << "Thyra::DefaultFiniteDifferenceModelEvaluator{";
00067   oss << "thyraModel=";
00068   if(thyraModel.get())
00069     oss << "\'"<<thyraModel->description()<<"\'";
00070   else
00071     oss << "NULL";
00072   oss << "}";
00073   return oss.str();
00074 }
00075 
00076 
00077 // Private functions overridden from ModelEvaulatorDefaultBase
00078 
00079 
00080 template<class Scalar>
00081 ModelEvaluatorBase::OutArgs<Scalar>
00082 DefaultFiniteDifferenceModelEvaluator<Scalar>::createOutArgsImpl() const
00083 {
00084   typedef ModelEvaluatorBase MEB;
00085   typedef MEB::DerivativeSupport DS;
00086   const RCP<const ModelEvaluator<Scalar> >
00087     thyraModel = this->getUnderlyingModel();
00088   const MEB::OutArgs<Scalar> wrappedOutArgs = thyraModel->createOutArgs();
00089   const int l_Np = wrappedOutArgs.Np(), l_Ng = wrappedOutArgs.Ng();
00090   MEB::OutArgsSetup<Scalar> outArgs;
00091   outArgs.setModelEvalDescription(this->description());
00092   outArgs.set_Np_Ng(l_Np,l_Ng);
00093   outArgs.setSupports(wrappedOutArgs);
00094   if (wrappedOutArgs.supports(MEB::OUT_ARG_f)) {
00095     for( int l = 0; l < l_Np; ++l ) {
00096       outArgs.setSupports(MEB::OUT_ARG_DfDp,l,MEB::DERIV_MV_BY_COL);
00097     }
00098   }
00099   for( int j = 0; j < l_Ng; ++j ) {
00100     for( int l = 0; l < l_Np; ++l ) {
00101       outArgs.setSupports( MEB::OUT_ARG_DgDp , j, l, MEB::DERIV_MV_BY_COL);
00102     }
00103   }
00104   // ToDo: Add support for more derivatives as needed!
00105   return outArgs;
00106 }
00107 
00108 
00109 template<class Scalar>
00110 void DefaultFiniteDifferenceModelEvaluator<Scalar>::evalModelImpl(
00111   const ModelEvaluatorBase::InArgs<Scalar> &inArgs,
00112   const ModelEvaluatorBase::OutArgs<Scalar> &outArgs
00113   ) const
00114 {
00115   using Teuchos::rcp;
00116   using Teuchos::rcp_const_cast;
00117   using Teuchos::rcp_dynamic_cast;
00118   using Teuchos::OSTab;
00119   typedef Teuchos::ScalarTraits<Scalar> ST;
00120   //typedef typename ST::magnitudeType ScalarMag;
00121   typedef ModelEvaluatorBase MEB;
00122   namespace DFDCT = DirectionalFiniteDiffCalculatorTypes;
00123 
00124   typedef RCP<VectorBase<Scalar> > V_ptr;
00125   typedef RCP<const VectorBase<Scalar> > CV_ptr;
00126   typedef RCP<MultiVectorBase<Scalar> > MV_ptr;
00127 
00128   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_BEGIN(
00129     "Thyra::DefaultFiniteDifferenceModelEvaluator",inArgs,outArgs
00130     );
00131 
00132   //
00133   // Note: Just do derivatives DfDp(l) and DgDp(j,l) for now!
00134   //
00135 
00136   const RCP<const VectorSpaceBase<Scalar> >
00137     p_space = thyraModel->get_p_space(0),
00138     g_space = thyraModel->get_g_space(0);
00139 
00140   //
00141   // A) Compute the base point
00142   //
00143 
00144   if(out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00145     *out << "\nComputing the base point ...\n";
00146 
00147   const int l_Np = outArgs.Np();
00148   const int l_Ng = outArgs.Ng();
00149   MEB::InArgs<Scalar> wrappedInArgs = inArgs;
00150   MEB::OutArgs<Scalar> baseFunc = thyraModel->createOutArgs();
00151   if( outArgs.supports(MEB::OUT_ARG_f) && outArgs.get_f().get() )
00152     baseFunc.set_f(outArgs.get_f());
00153   for( int j = 0; j < l_Ng; ++j ) {
00154     V_ptr g_j;
00155     if( (g_j=outArgs.get_g(j)).get() )
00156       baseFunc.set_g(j,g_j);
00157   }
00158   // 2007/08/27: We really should really try to allow some derivatives to pass
00159   // through and some derivatives to be computed by finite differences. Right
00160   // now, if you use this class, all derivatives w.r.t. parameters are finite
00161   // differenced and that is not given the user enough control!
00162 
00163   thyraModel->evalModel(wrappedInArgs,baseFunc);
00164 
00165   bool failed = baseFunc.isFailed();
00166 
00167   //
00168   // B) Compute the derivatives
00169   //
00170  
00171   if(!failed) {
00172 
00173     // a) Determine what derivatives you need to support first
00174 
00175     Array<int> compute_DfDp;
00176     Array<Array<int> > compute_DgDp(l_Ng);
00177     DFDCT::SelectedDerivatives selectedDerivs;
00178 
00179     for ( int l = 0; l < l_Np; ++l ) {
00180 
00181       // DfDp(l)
00182       if(
00183         outArgs.supports(MEB::OUT_ARG_DfDp,l).none()==false
00184         &&
00185         outArgs.get_DfDp(l).isEmpty()==false
00186         )
00187       {
00188         selectedDerivs.supports(MEB::OUT_ARG_DfDp,l);
00189         compute_DfDp.push_back(true);
00190       }
00191       else
00192       {
00193         compute_DfDp.push_back(false);
00194       }
00195 
00196       // DgDp(j=0...,l)
00197       for ( int j = 0; j < l_Ng; ++j ) {
00198         if(
00199           outArgs.supports(MEB::OUT_ARG_DgDp,j,l).none()==false
00200           &&
00201           outArgs.get_DgDp(j,l).isEmpty()==false
00202           )
00203         {
00204           selectedDerivs.supports(MEB::OUT_ARG_DgDp,j,l);
00205           compute_DgDp[j].push_back(true);
00206         }
00207         else
00208         {
00209           compute_DgDp[j].push_back(false);
00210         }
00211       }
00212     }
00213 
00214     // b) Create the deriv OutArgs and set the output objects that need to be
00215     // computed with finite differences
00216  
00217     MEB::OutArgs<Scalar>
00218       deriv = direcFiniteDiffCalculator_->createOutArgs(
00219         *thyraModel, selectedDerivs );
00220 
00221     for ( int l = 0; l < l_Np; ++l ) {
00222       if ( compute_DfDp[l] )
00223         deriv.set_DfDp(l,outArgs.get_DfDp(l));
00224       for ( int j = 0; j < l_Ng; ++j ) {
00225         if ( compute_DgDp[j][l] )
00226           deriv.set_DgDp(j,l,outArgs.get_DgDp(j,l));
00227       }
00228     }
00229 
00230     // c) Compute the missing functions with finite differences!
00231 
00232     direcFiniteDiffCalculator_->calcDerivatives(
00233       *thyraModel,inArgs,baseFunc,deriv
00234       );
00235 
00236   }
00237 
00238   if(failed) {
00239     if(out.get() && includesVerbLevel(verbLevel,Teuchos::VERB_LOW))
00240       *out << "\nEvaluation failed, returning NaNs ...\n";
00241     outArgs.setFailed();
00242   }
00243 
00244   THYRA_MODEL_EVALUATOR_DECORATOR_EVAL_MODEL_END();
00245  
00246 }
00247 
00248 
00249 } // namespace Thyra
00250 
00251 
00252 #endif // THYRA_DEFAULT_FINITE_DIFFERENCE_MODEL_EVALUATOR_DEF_HPP
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines