NLPInterfacePack_NLPDirectThyraModelEvaluator.cpp

Go to the documentation of this file.
00001 // @HEADER
00002 // ***********************************************************************
00003 // 
00004 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
00005 //                  Copyright (2003) 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 Roscoe A. Bartlett (rabartl@sandia.gov) 
00025 // 
00026 // ***********************************************************************
00027 // @HEADER
00028 
00029 #include "NLPInterfacePack_NLPDirectThyraModelEvaluator.hpp"
00030 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00031 #include "AbstractLinAlgPack_VectorOut.hpp"
00032 #include "AbstractLinAlgPack_ThyraAccessors.hpp"
00033 #include "AbstractLinAlgPack_VectorSpaceThyra.hpp"
00034 #include "AbstractLinAlgPack_VectorMutableThyra.hpp"
00035 #include "AbstractLinAlgPack_MultiVectorMutableThyra.hpp"
00036 #include "AbstractLinAlgPack_MatrixOpNonsingThyra.hpp"
00037 #include "AbstractLinAlgPack_VectorSpaceBlocked.hpp"
00038 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00039 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
00040 #include "AbstractLinAlgPack_BasisSystem.hpp"
00041 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00042 #include "Thyra_ModelEvaluatorHelpers.hpp"
00043 #include "Thyra_DetachedVectorView.hpp"
00044 #include "Teuchos_AbstractFactoryStd.hpp"
00045 #include "Teuchos_TestForException.hpp"
00046 #include "Teuchos_dyn_cast.hpp"
00047 
00048 namespace NLPInterfacePack {
00049 
00050 NLPDirectThyraModelEvaluator::NLPDirectThyraModelEvaluator()
00051 {}
00052 
00053 NLPDirectThyraModelEvaluator::NLPDirectThyraModelEvaluator(
00054   const Teuchos::RCP<Thyra::ModelEvaluator<value_type> >  &model   
00055   ,const int                                                      p_idx 
00056   ,const int                                                      g_idx 
00057   ,const objDirecFiniteDiffCalculator_ptr_t                       objDirecFiniteDiffCalculator
00058   ,const conDirecFiniteDiffCalculator_ptr_t                       conDirecFiniteDiffCalculator
00059   )
00060 {
00061   initialize(
00062     model,p_idx,g_idx
00063     ,objDirecFiniteDiffCalculator,conDirecFiniteDiffCalculator
00064     );
00065 }
00066 
00067 void NLPDirectThyraModelEvaluator::initialize(
00068   const Teuchos::RCP<Thyra::ModelEvaluator<value_type> >  &model
00069   ,const int                                                      p_idx
00070   ,const int                                                      g_idx
00071   ,const objDirecFiniteDiffCalculator_ptr_t                       objDirecFiniteDiffCalculator
00072   ,const conDirecFiniteDiffCalculator_ptr_t                       conDirecFiniteDiffCalculator
00073   )
00074 {
00075   typedef Thyra::ModelEvaluatorBase MEB;
00076   if(objDirecFiniteDiffCalculator.get())
00077     this->set_objDirecFiniteDiffCalculator(objDirecFiniteDiffCalculator);
00078   if(conDirecFiniteDiffCalculator.get())
00079     this->set_conDirecFiniteDiffCalculator(conDirecFiniteDiffCalculator);
00080   initializeBase(model,p_idx,g_idx);
00081   Thyra::ModelEvaluatorBase::OutArgs<double> model_outArgs = model->createOutArgs();
00082   MEB::DerivativeProperties model_W_properties = model_outArgs.get_W_properties();
00083   if( p_idx >= 0 ) {
00084     TEST_FOR_EXCEPTION(
00085       (conDirecFiniteDiffCalculator_.get()==0  && !model_outArgs.supports(MEB::OUT_ARG_DfDp,p_idx).supports(MEB::DERIV_MV_BY_COL))
00086       ,std::invalid_argument
00087       ,"Error, model must support computing DfDp("<<p_idx<<") as a"
00088       " column-oriented multi-vector if not using finite differences!"
00089       );
00090   }
00091 }
00092 
00093 // Overridden public members from NLP
00094 
00095 void NLPDirectThyraModelEvaluator::initialize(bool test_setup)
00096 {
00097   if(initialized_) {
00098     NLPDirect::initialize(test_setup);
00099     return;
00100   }
00101   NLPThyraModelEvaluatorBase::initialize(test_setup);
00102   NLPDirect::initialize(test_setup);
00103 }
00104 
00105 void NLPDirectThyraModelEvaluator::unset_quantities()
00106 {
00107   NLPDirect::unset_quantities();
00108 }
00109 
00110 // Overridden public members from NLPObjGrad
00111 
00112 bool NLPDirectThyraModelEvaluator::supports_Gf() const
00113 {
00114   if(objDirecFiniteDiffCalculator_.get())
00115     return false;
00116   return true;
00117   // ToDo: Change this to false when only the operator for model_DgDx is
00118   // supported!
00119 }
00120 
00121 bool NLPDirectThyraModelEvaluator::supports_Gf_prod() const
00122 {
00123   return ( objDirecFiniteDiffCalculator_.get() != NULL );
00124 }
00125 
00126 value_type NLPDirectThyraModelEvaluator::calc_Gf_prod(
00127   const Vector& x, const Vector& d, bool newx
00128   ) const
00129 {
00130   TEST_FOR_EXCEPT(objDirecFiniteDiffCalculator_.get()==0);
00131   typedef Thyra::ModelEvaluatorBase MEB;
00132   MEB::InArgs<value_type>  basePoint   = model_->createInArgs();
00133   MEB::InArgs<value_type>  directions  = model_->createInArgs();
00134   MEB::OutArgs<value_type> baseFunc    = model_->createOutArgs();
00135   MEB::OutArgs<value_type> variations  = model_->createOutArgs();
00136   set_x(x,&basePoint);
00137   set_x(d,&directions);
00138   variations.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00139   objDirecFiniteDiffCalculator_->calcVariations(
00140     *model_,basePoint,directions,baseFunc,variations
00141     );
00142   return Thyra::get_ele(*variations.get_g(g_idx_),0);
00143 }
00144 
00145 // Overridden public members from NLPDirect
00146 
00147 Range1D NLPDirectThyraModelEvaluator::var_dep() const
00148 {
00149   return basis_sys_->var_dep();
00150 }
00151 
00152 Range1D NLPDirectThyraModelEvaluator::var_indep() const
00153 {
00154   return basis_sys_->var_indep();
00155 }
00156 
00157 const NLPDirect::mat_fcty_ptr_t
00158 NLPDirectThyraModelEvaluator::factory_D() const
00159 {
00160   return basis_sys_->factory_D();
00161 }
00162 
00163 const NLPDirect::mat_sym_nonsing_fcty_ptr_t
00164 NLPDirectThyraModelEvaluator::factory_S() const
00165 {
00166   return basis_sys_->factory_S();
00167 }
00168 
00169 void NLPDirectThyraModelEvaluator::calc_point(
00170   const Vector     &x
00171   ,value_type      *f
00172   ,VectorMutable   *c
00173   ,bool            recalc_c
00174   ,VectorMutable   *Gf
00175   ,VectorMutable   *py
00176   ,VectorMutable   *rGf
00177   ,MatrixOp        *GcU
00178   ,MatrixOp        *D
00179   ,MatrixOp        *Uz
00180   ) const
00181 {
00182   using Teuchos::FancyOStream;
00183   using Teuchos::OSTab;
00184   using Teuchos::dyn_cast;
00185   using Teuchos::RCP;
00186   using Teuchos::rcp;
00187   using Teuchos::rcp_const_cast;
00188   using Teuchos::rcp_dynamic_cast;
00189   using AbstractLinAlgPack::VectorSpaceThyra;
00190   using AbstractLinAlgPack::VectorMutableThyra;
00191   using AbstractLinAlgPack::MultiVectorMutableThyra;
00192   using AbstractLinAlgPack::MatrixOpThyra;
00193   using AbstractLinAlgPack::MatrixOpNonsingThyra;
00194   typedef Thyra::ModelEvaluatorBase MEB;
00195   typedef Teuchos::RCP<Thyra::VectorBase<value_type> > TVecPtr;
00196   typedef MEB::DerivativeMultiVector<value_type> DerivMV;
00197   typedef MEB::Derivative<value_type> Deriv;
00198   //
00199   // Get output and verbosity
00200   //
00201   const Teuchos::RCP<Teuchos::FancyOStream>
00202     out = this->getOStream();
00203   const Teuchos::EVerbosityLevel
00204     verbLevel = ( showModelEvaluatorTrace() ? this->getVerbLevel() : Teuchos::VERB_NONE );
00205   Teuchos::OSTab tab(out);
00206   typedef Teuchos::VerboseObjectTempState<MEB> VOTSME;
00207   VOTSME modelOutputTempState(model_,out,verbLevel);
00208   typedef Teuchos::VerboseObjectTempState<Thyra::DirectionalFiniteDiffCalculator<value_type> > VOTSDFDC;
00209   VOTSDFDC objDirecFiniteDiffCalculatorOutputTempState(objDirecFiniteDiffCalculator_,out,verbLevel);
00210   VOTSDFDC conDirecFiniteDiffCalculatorOutputTempState(conDirecFiniteDiffCalculator_,out,verbLevel);
00211   const bool trace = static_cast<int>(verbLevel) >= static_cast<int>(Teuchos::VERB_LOW);
00212   if(out.get() && trace)
00213     *out << "\nEntering MoochoPack::NLPDirectThyraModelEvaluator::calc_point(...) ...\n";
00214   //
00215   // Validate input
00216   //
00217   TEST_FOR_EXCEPT(GcU!=NULL);  // Can't handle these yet!
00218   TEST_FOR_EXCEPT(Uz!=NULL);
00219   TEST_FOR_EXCEPTION(
00220     objDirecFiniteDiffCalculator_.get()!=NULL && Gf!=NULL, std::logic_error
00221     ,"Error, can not compute full gradient vector Gf when using directional finite differences!"
00222     );
00223   //
00224   // Compute using derivative objects that come from the underlying model.
00225   //
00226   // Set the input and output arguments
00227   //
00228   // ToDo: Disallow computation Gf and instead just compute
00229   // the operators for this!
00230   //
00231   MEB::InArgs<value_type>  model_inArgs  = model_->createInArgs();
00232   MEB::OutArgs<value_type> model_outArgs = model_->createOutArgs();
00233   NLPObjGrad::ObjGradInfo obj_grad_info;
00234   obj_grad_info.Gf = Gf;
00235   obj_grad_info.f = f;
00236   if(recalc_c) obj_grad_info.c = c;
00237   preprocessBaseInOutArgs(
00238     x,true,NULL,&obj_grad_info,NULL
00239     ,&model_inArgs,&model_outArgs,NULL,NULL,NULL,NULL
00240     );
00241   if( py || rGf || D ) {
00242     if(thyra_C_.get()==NULL)
00243       thyra_C_ = model_->create_W();
00244     model_outArgs.set_W(thyra_C_.assert_not_null());
00245   }
00246   bool new_thyra_N = false;
00247   if( rGf || D ) {
00248     if(thyra_N_.get()==NULL) {
00249       thyra_N_ = Thyra::create_DfDp_mv(*model_,p_idx_,MEB::DERIV_MV_BY_COL).getMultiVector();
00250       new_thyra_N = true;
00251     }
00252     if(
00253       !conDirecFiniteDiffCalculator_.get()
00254       &&
00255       ( new_thyra_N
00256         || model_outArgs.get_DfDp_properties(p_idx_).linearity!=MEB::DERIV_LINEARITY_CONST )
00257       )
00258     {
00259       model_outArgs.set_DfDp(p_idx_,DerivMV(thyra_N_.assert_not_null(),MEB::DERIV_MV_BY_COL));
00260     }
00261   }
00262   //
00263   // Evaluate the functions
00264   //
00265   model_->evalModel(model_inArgs,model_outArgs);
00266   //
00267   // Postprocess the evaluation
00268   //
00269   postprocessBaseOutArgs(&model_outArgs,Gf,f,recalc_c?c:NULL);
00270   // Setup solve components
00271   const VectorSpaceThyra                                       *space_c;
00272   const VectorSpaceThyra                                       *space_xD;
00273   Teuchos::RCP<const Thyra::VectorBase<value_type> >   thyra_c;
00274   Teuchos::RCP<Thyra::VectorBase<value_type> >         thyra_py;
00275   RCP<MatrixOp>                                        D_used = rcp(D,false);
00276   RCP<Thyra::MultiVectorBase<value_type> >             thyra_D;
00277   if( py || ( ( rGf || D ) && conDirecFiniteDiffCalculator_.get() ) ) {
00278     space_c  = &dyn_cast<const VectorSpaceThyra>(c->space()),
00279       space_xD = &dyn_cast<const VectorSpaceThyra>(py->space());
00280     get_thyra_vector(*space_c,*c,&thyra_c);
00281     get_thyra_vector(*space_xD,py,&thyra_py);
00282   }
00283   if( D || rGf ) {
00284     if(!D) D_used = this->factory_D()->create();
00285     thyra_D =
00286       rcp_const_cast<Thyra::MultiVectorBase<value_type> >(
00287         rcp_dynamic_cast<const Thyra::MultiVectorBase<value_type> >(
00288           dyn_cast<MultiVectorMutableThyra>(*D_used).thyra_multi_vec()
00289           )
00290         );
00291   }
00292   if(
00293     ( rGf || D )
00294     &&
00295     conDirecFiniteDiffCalculator_.get()
00296     &&
00297     ( new_thyra_N || DfDp_is_const() )
00298     )
00299   {
00300     if(out.get() && trace)
00301       *out << "\nComputing thyra_N using directional finite differences ...\n";
00302     // !
00303     if(thyra_c.get())
00304       model_outArgs.set_f(rcp_const_cast<Thyra::VectorBase<value_type> >(thyra_c)); // Okay, will not be changed!
00305     typedef Thyra::DirectionalFiniteDiffCalculator<value_type>::SelectedDerivatives SelectedDerivatives; 
00306     MEB::OutArgs<value_type>
00307       model_fdOutArgs = conDirecFiniteDiffCalculator_->createOutArgs(
00308         *model_,SelectedDerivatives().supports(MEB::OUT_ARG_DfDp,0)
00309         );
00310     model_fdOutArgs.set_DfDp(p_idx_,DerivMV(thyra_N_.assert_not_null(),MEB::DERIV_MV_BY_COL));
00311     conDirecFiniteDiffCalculator_->calcDerivatives(
00312       *model_
00313       ,model_inArgs     // The base point to compute the derivatives at
00314       ,model_outArgs    // The base function values
00315       ,model_fdOutArgs  // The derivatives that we want to compute
00316       );
00317   }
00318   // Perform solve
00319   if( ( D || rGf ) && py ) {
00320     // Solve for py and D together
00321     if(out.get() && trace)
00322       *out << "\nSolving C*[py,D] = -[c,N] simultaneously ...\n";
00323     const int nind = thyra_N_->domain()->dim();
00324     RCP<Thyra::MultiVectorBase<value_type> > 
00325       thyra_cN = Thyra::createMembers(thyra_N_->range(),nind+1);
00326     Thyra::assign(&*thyra_cN->col(0),*thyra_c);
00327     Thyra::assign(&*thyra_cN->subView(Teuchos::Range1D(1,nind)),*thyra_N_);
00328     RCP<Thyra::MultiVectorBase<value_type> > 
00329       thyra_pyD = Thyra::createMembers(thyra_D->range(),nind+1);
00330     Thyra::assign(&*thyra_pyD,0.0);
00331     Thyra::SolveStatus<value_type>
00332       solveStatus = Thyra::solve(*thyra_C_,Thyra::NOTRANS,*thyra_cN,&*thyra_pyD);
00333     if(out.get() && trace) {
00334       *out
00335         << "\nsolve status:\n";
00336       OSTab(out).o() << solveStatus;
00337     }
00338     Thyra::scale(-1.0,&*thyra_pyD);
00339     Thyra::assign(&*thyra_py,*thyra_pyD->col(0));
00340     Thyra::assign(&*thyra_D,*thyra_pyD->subView(Teuchos::Range1D(1,nind)));
00341   }
00342   else {
00343     // Solve for py or D
00344     if( py ) {
00345       if(out.get() && trace)
00346         *out << "\nSolving C*py = -c ...\n";
00347       Thyra::assign(&*thyra_py,0.0);
00348       Thyra::SolveStatus<value_type>
00349         solveStatus = Thyra::solve(*thyra_C_,Thyra::NOTRANS,*thyra_c,&*thyra_py);
00350       if(out.get() && trace) {
00351         *out
00352           << "\nsolve status:\n";
00353         OSTab(out).o() << solveStatus;
00354       }
00355       Thyra::Vt_S(&*thyra_py,-1.0);
00356     }
00357     if( D || rGf ) {
00358       if(out.get() && trace)
00359         *out << "\nSolving C*D = -N ...\n";
00360       Thyra::assign(&*thyra_D,0.0);
00361       Thyra::SolveStatus<value_type>
00362         solveStatus = Thyra::solve(*thyra_C_,Thyra::NOTRANS,*thyra_N_,&*thyra_D);
00363       if(out.get() && trace) {
00364         *out
00365           << "\nsolve status:\n";
00366         OSTab(out).o() << solveStatus;
00367       }
00368       Thyra::scale(-1.0,&*thyra_D);
00369     }
00370   }
00371   if(thyra_py.get()) {
00372     free_thyra_vector(*space_c,*c,&thyra_c);
00373     commit_thyra_vector(*space_xD,py,&thyra_py);
00374   }
00375   // Compute reduced gradient
00376   if(rGf) {
00377     // rGf = D' * Gf_xD + Gf_xI
00378     const Range1D
00379       var_dep   = basis_sys_->var_dep(),
00380       var_indep = basis_sys_->var_indep();
00381     if(objDirecFiniteDiffCalculator_.get()) {
00382       if(out.get() && trace)
00383         *out << "\nComputing rGf using directional finite differences ...\n";
00384       //
00385       // Compute each component of the reduced gradient as
00386       //
00387       //    rGf(i) = fd_product( model.g, [ D(:,i); e(i) ] )
00388       //
00389       AbstractLinAlgPack::VectorDenseMutableEncap d_rGf(*rGf);
00390       TVecPtr e_i = createMember(model_->get_p_space(p_idx_));
00391       Thyra::assign(&*e_i,0.0);
00392       MEB::InArgs<value_type> dir = model_->createInArgs();
00393       dir.set_p(p_idx_,e_i);
00394       MEB::OutArgs<value_type> bfunc = model_->createOutArgs();
00395       if(f) {
00396         bfunc.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00397         Thyra::set_ele(0,*f,&*bfunc.get_g(g_idx_));
00398       }
00399       MEB::OutArgs<value_type> var = model_->createOutArgs();
00400       var.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00401       const int np = var_indep.size();
00402       for( int i = 0; i < np; ++i ) {
00403         set_ele(i,1.0,&*e_i);
00404         dir.set_x(thyra_D->col(i));
00405         objDirecFiniteDiffCalculator_->calcVariations(
00406           *model_,model_inArgs,dir,bfunc,var
00407           );
00408         dir.set_x(Teuchos::null);
00409         Thyra::set_ele(i,0.0,&*e_i);
00410         d_rGf()[i] = Thyra::get_ele(*var.get_g(g_idx_),0);
00411       }
00412     }
00413     else {
00414       LinAlgOpPack::V_MtV( rGf, *D_used, BLAS_Cpp::trans, *Gf->sub_view(var_dep) );
00415       LinAlgOpPack::Vp_V( rGf, *Gf->sub_view(var_indep) );
00416       // ToDo: Just compute the operators associated with Gf and not Gf directly!
00417     }
00418   }
00419   // * ToDo: Add specialized algorithm for computing D using an inexact Jacobian
00420   // * ToDo: Add in logic for inexact solves
00421   if(out.get() && trace)
00422     *out << "\nLeaving MoochoPack::NLPDirectThyraModelEvaluator::calc_point(...) ...\n";
00423 }
00424 
00425 void NLPDirectThyraModelEvaluator::calc_semi_newton_step(
00426   const Vector    &x
00427   ,VectorMutable  *c
00428   ,bool           recalc_c
00429   ,VectorMutable  *py
00430   ) const
00431 {
00432   if(recalc_c) {
00433     //
00434     // Recompute c
00435     //
00436     TEST_FOR_EXCEPT(true);
00437   }
00438   // Compute py = - inv(C)*c
00439   TEST_FOR_EXCEPT(true);
00440 }
00441 
00442 } // end namespace NLPInterfacePack

Generated on Tue Oct 20 12:51:48 2009 for MOOCHO (Single Doxygen Collection) by doxygen 1.4.7