NLPInterfacePack_NLPDirectThyraModelEvaluator.cpp

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   if ( !is_null(model_outArgs.get_W()) || !is_null(model_outArgs.get_W_op()) ) {
00263     if (model_inArgs.supports(MEB::IN_ARG_alpha))
00264       model_inArgs.set_alpha(0.0);
00265     if (model_inArgs.supports(MEB::IN_ARG_beta))
00266       model_inArgs.set_beta(1.0);
00267   }
00268   //
00269   // Evaluate the functions
00270   //
00271   model_->evalModel(model_inArgs,model_outArgs);
00272   //
00273   // Postprocess the evaluation
00274   //
00275   postprocessBaseOutArgs(&model_outArgs,Gf,f,recalc_c?c:NULL);
00276   // Setup solve components
00277   const VectorSpaceThyra *space_c;
00278   const VectorSpaceThyra *space_xD;
00279   Teuchos::RCP<const Thyra::VectorBase<value_type> > thyra_c;
00280   Teuchos::RCP<Thyra::VectorBase<value_type> > thyra_py;
00281   RCP<MatrixOp> D_used = rcp(D,false);
00282   RCP<Thyra::MultiVectorBase<value_type> > thyra_D;
00283   if( py || ( ( rGf || D ) && conDirecFiniteDiffCalculator_.get() ) ) {
00284     space_c = &dyn_cast<const VectorSpaceThyra>(c->space()),
00285       space_xD = &dyn_cast<const VectorSpaceThyra>(py->space());
00286     get_thyra_vector(*space_c,*c,&thyra_c);
00287     get_thyra_vector(*space_xD,py,&thyra_py);
00288   }
00289   if( D || rGf ) {
00290     if(!D) D_used = this->factory_D()->create();
00291     thyra_D =
00292       rcp_const_cast<Thyra::MultiVectorBase<value_type> >(
00293         rcp_dynamic_cast<const Thyra::MultiVectorBase<value_type> >(
00294           dyn_cast<MultiVectorMutableThyra>(*D_used).thyra_multi_vec()
00295           )
00296         );
00297   }
00298   if(
00299     ( rGf || D )
00300     &&
00301     conDirecFiniteDiffCalculator_.get()
00302     &&
00303     ( new_thyra_N || DfDp_is_const() )
00304     )
00305   {
00306     if(out.get() && trace)
00307       *out << "\nComputing thyra_N using directional finite differences ...\n";
00308     // !
00309     if(thyra_c.get())
00310       model_outArgs.set_f(rcp_const_cast<Thyra::VectorBase<value_type> >(thyra_c)); // Okay, will not be changed!
00311     typedef Thyra::DirectionalFiniteDiffCalculator<value_type>::SelectedDerivatives SelectedDerivatives; 
00312     MEB::OutArgs<value_type>
00313       model_fdOutArgs = conDirecFiniteDiffCalculator_->createOutArgs(
00314         *model_,SelectedDerivatives().supports(MEB::OUT_ARG_DfDp,0)
00315         );
00316     model_fdOutArgs.set_DfDp(p_idx_,DerivMV(thyra_N_.assert_not_null(),MEB::DERIV_MV_BY_COL));
00317     conDirecFiniteDiffCalculator_->calcDerivatives(
00318       *model_
00319       ,model_inArgs // The base point to compute the derivatives at
00320       ,model_outArgs // The base function values
00321       ,model_fdOutArgs // The derivatives that we want to compute
00322       );
00323   }
00324   // Perform solve
00325   if( ( D || rGf ) && py ) {
00326     // Solve for py and D together
00327     if(out.get() && trace)
00328       *out << "\nSolving C*[py,D] = -[c,N] simultaneously ...\n";
00329     const int nind = thyra_N_->domain()->dim();
00330     RCP<Thyra::MultiVectorBase<value_type> > 
00331       thyra_cN = Thyra::createMembers(thyra_N_->range(),nind+1);
00332     Thyra::assign(&*thyra_cN->col(0),*thyra_c);
00333     Thyra::assign(&*thyra_cN->subView(Teuchos::Range1D(1,nind)),*thyra_N_);
00334     RCP<Thyra::MultiVectorBase<value_type> > 
00335       thyra_pyD = Thyra::createMembers(thyra_D->range(),nind+1);
00336     Thyra::assign(&*thyra_pyD,0.0);
00337     Thyra::SolveStatus<value_type>
00338       solveStatus = Thyra::solve(*thyra_C_,Thyra::NOTRANS,*thyra_cN,&*thyra_pyD);
00339     if(out.get() && trace) {
00340       *out
00341         << "\nsolve status:\n";
00342       OSTab(out).o() << solveStatus;
00343     }
00344     Thyra::scale(-1.0,&*thyra_pyD);
00345     Thyra::assign(&*thyra_py,*thyra_pyD->col(0));
00346     Thyra::assign(&*thyra_D,*thyra_pyD->subView(Teuchos::Range1D(1,nind)));
00347   }
00348   else {
00349     // Solve for py or D
00350     if( py ) {
00351       if(out.get() && trace)
00352         *out << "\nSolving C*py = -c ...\n";
00353       Thyra::assign(&*thyra_py,0.0);
00354       Thyra::SolveStatus<value_type>
00355         solveStatus = Thyra::solve(*thyra_C_,Thyra::NOTRANS,*thyra_c,&*thyra_py);
00356       if(out.get() && trace) {
00357         *out
00358           << "\nsolve status:\n";
00359         OSTab(out).o() << solveStatus;
00360       }
00361       Thyra::Vt_S(&*thyra_py,-1.0);
00362     }
00363     if( D || rGf ) {
00364       if(out.get() && trace)
00365         *out << "\nSolving C*D = -N ...\n";
00366       Thyra::assign(&*thyra_D,0.0);
00367       Thyra::SolveStatus<value_type>
00368         solveStatus = Thyra::solve(*thyra_C_,Thyra::NOTRANS,*thyra_N_,&*thyra_D);
00369       if(out.get() && trace) {
00370         *out
00371           << "\nsolve status:\n";
00372         OSTab(out).o() << solveStatus;
00373       }
00374       Thyra::scale(-1.0,&*thyra_D);
00375     }
00376   }
00377   if(thyra_py.get()) {
00378     free_thyra_vector(*space_c,*c,&thyra_c);
00379     commit_thyra_vector(*space_xD,py,&thyra_py);
00380   }
00381   // Compute reduced gradient
00382   if(rGf) {
00383     // rGf = D' * Gf_xD + Gf_xI
00384     const Range1D
00385       var_dep = basis_sys_->var_dep(),
00386       var_indep = basis_sys_->var_indep();
00387     if(objDirecFiniteDiffCalculator_.get()) {
00388       if(out.get() && trace)
00389         *out << "\nComputing rGf using directional finite differences ...\n";
00390       //
00391       // Compute each component of the reduced gradient as
00392       //
00393       // rGf(i) = fd_product( model.g, [ D(:,i); e(i) ] )
00394       //
00395       AbstractLinAlgPack::VectorDenseMutableEncap d_rGf(*rGf);
00396       TVecPtr e_i = createMember(model_->get_p_space(p_idx_));
00397       Thyra::assign(&*e_i,0.0);
00398       MEB::InArgs<value_type> dir = model_->createInArgs();
00399       dir.set_p(p_idx_,e_i);
00400       MEB::OutArgs<value_type> bfunc = model_->createOutArgs();
00401       if(f) {
00402         bfunc.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00403         Thyra::set_ele(0,*f,&*bfunc.get_g(g_idx_));
00404       }
00405       MEB::OutArgs<value_type> var = model_->createOutArgs();
00406       var.set_g(g_idx_,createMember(model_->get_g_space(g_idx_)));
00407       const int np = var_indep.size();
00408       for( int i = 0; i < np; ++i ) {
00409         set_ele(i,1.0,&*e_i);
00410         dir.set_x(thyra_D->col(i));
00411         objDirecFiniteDiffCalculator_->calcVariations(
00412           *model_,model_inArgs,dir,bfunc,var
00413           );
00414         dir.set_x(Teuchos::null);
00415         Thyra::set_ele(i,0.0,&*e_i);
00416         d_rGf()[i] = Thyra::get_ele(*var.get_g(g_idx_),0);
00417       }
00418     }
00419     else {
00420       LinAlgOpPack::V_MtV( rGf, *D_used, BLAS_Cpp::trans, *Gf->sub_view(var_dep) );
00421       LinAlgOpPack::Vp_V( rGf, *Gf->sub_view(var_indep) );
00422       // ToDo: Just compute the operators associated with Gf and not Gf directly!
00423     }
00424   }
00425   // * ToDo: Add specialized algorithm for computing D using an inexact Jacobian
00426   // * ToDo: Add in logic for inexact solves
00427   if(out.get() && trace)
00428     *out << "\nLeaving MoochoPack::NLPDirectThyraModelEvaluator::calc_point(...) ...\n";
00429 }
00430 
00431 void NLPDirectThyraModelEvaluator::calc_semi_newton_step(
00432   const Vector &x
00433   ,VectorMutable *c
00434   ,bool recalc_c
00435   ,VectorMutable *py
00436   ) const
00437 {
00438   if(recalc_c) {
00439     //
00440     // Recompute c
00441     //
00442     TEST_FOR_EXCEPT(true);
00443   }
00444   // Compute py = - inv(C)*c
00445   TEST_FOR_EXCEPT(true);
00446 }
00447 
00448 } // end namespace NLPInterfacePack

Generated on Wed May 12 21:51:33 2010 for MOOCHO/Thyra Adapter Software by  doxygen 1.4.7