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

Generated on Thu Sep 18 12:34:37 2008 for MOOCHO/Thyra Adapter Software by doxygen 1.3.9.1