MOOCHO (Single Doxygen Collection) Version of the Day
NLPInterfacePack_ExampleNLPDirect.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 <assert.h>
00030 
00031 #include <stdexcept>
00032 
00033 #include "NLPInterfacePack_ExampleNLPDirect.hpp"
00034 #include "ExampleNLPDirectRTOps.h"
00035 #include "AbstractLinAlgPack_BasisSystemComposite.hpp"
00036 #include "AbstractLinAlgPack_VectorMutable.hpp"
00037 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00038 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00039 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00040 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00041 #include "RTOpPack_RTOpC.hpp"
00042 #include "Teuchos_dyn_cast.hpp"
00043 #include "Teuchos_TestForException.hpp"
00044 #include "Teuchos_AbstractFactoryStd.hpp"
00045 
00046 namespace {
00047 
00048 static RTOpPack::RTOpC explnlp2_calc_py_D_op;
00049 
00050 class init_rtop_server_t {
00051 public:
00052   init_rtop_server_t() {
00053     TEST_FOR_EXCEPT(0!=RTOp_TOp_explnlp2_calc_py_D_construct(0,&explnlp2_calc_py_D_op.op()));
00054   }
00055 }; 
00056 init_rtop_server_t  init_rtop_server;
00057 
00058 } // end namespace
00059 
00060 namespace NLPInterfacePack {
00061 
00062 ExampleNLPDirect::ExampleNLPDirect(
00063   const VectorSpace::space_ptr_t&  vec_space
00064   ,value_type                      xo
00065   ,bool                            has_bounds
00066   ,bool                            dep_bounded
00067   )
00068   :ExampleNLPObjGrad(vec_space,xo,has_bounds,dep_bounded),initialized_(false)
00069 {
00070   namespace rcp = MemMngPack;
00071 
00072   // Create the factory object for D
00073   factory_D_ = Teuchos::rcp(new Teuchos::AbstractFactoryStd<MatrixOp,MatrixSymDiagStd>());
00074   NLPDirect::set_factories(
00075     Teuchos::rcp(
00076       new Teuchos::AbstractFactoryStd<MatrixSymOp,MatrixSymDiagStd>())               // D'*D
00077     ,Teuchos::rcp(
00078       new Teuchos::AbstractFactoryStd<MatrixSymOpNonsing,MatrixSymDiagStd>())    // S
00079     );
00080 }
00081 
00082 // Overridden public members from NLP
00083 
00084 void ExampleNLPDirect::initialize(bool test_setup)
00085 {
00086 
00087   if( initialized_ ) {
00088     NLPDirect::initialize(test_setup);
00089     ExampleNLPObjGrad::initialize(test_setup);
00090     return;
00091   }
00092 
00093   NLPDirect::initialize(test_setup);
00094   ExampleNLPObjGrad::initialize(test_setup);
00095 
00096   initialized_ = true;
00097 }
00098 
00099 bool ExampleNLPDirect::is_initialized() const
00100 {
00101   return initialized_;
00102 }
00103 
00104 // Overridden public members from NLPDirect
00105 
00106 Range1D ExampleNLPDirect::var_dep() const
00107 {
00108   return ExampleNLPObjGrad::var_dep();
00109 }
00110 
00111 Range1D ExampleNLPDirect::var_indep() const
00112 {
00113   return ExampleNLPObjGrad::var_indep();
00114 }
00115 
00116 const NLPDirect::mat_fcty_ptr_t
00117 ExampleNLPDirect::factory_D() const
00118 {
00119   return factory_D_;
00120 }
00121 
00122 void ExampleNLPDirect::calc_point(
00123   const Vector     &x
00124   ,value_type      *f
00125   ,VectorMutable   *c
00126   ,bool            recalc_c
00127   ,VectorMutable   *Gf
00128   ,VectorMutable   *py
00129   ,VectorMutable   *rGf
00130   ,MatrixOp        *GcU
00131   ,MatrixOp        *D
00132   ,MatrixOp        *Uz
00133   ) const
00134 {
00135   using Teuchos::dyn_cast;
00136   using LinAlgOpPack::Vp_MtV;
00137 
00138   assert_is_initialized();
00139 
00140   const size_type
00141     n = this->n();
00142 
00143   // Validate the input
00144 
00145 #ifdef TEUCHOS_DEBUG
00146   TEST_FOR_EXCEPTION(
00147     x.dim() != n, std::invalid_argument
00148     ,"ExampleNLPDirect::calc_point(...), Error x.dim() = " << x.dim()
00149     << " != n = " << n );
00150   TEST_FOR_EXCEPTION(
00151     c && !this->space_c()->is_compatible(c->space()), std::invalid_argument
00152     ,"ExampleNLPDirect::calc_point(...), Error c is not compatible" );
00153   TEST_FOR_EXCEPTION(
00154     Gf && !this->space_x()->is_compatible(Gf->space()), std::invalid_argument
00155     ,"ExampleNLPDirect::calc_point(...), Error, Gf is not compatible" );
00156   TEST_FOR_EXCEPTION(
00157     py && !this->space_x()->sub_space(this->var_dep())->is_compatible(py->space()), std::invalid_argument
00158     ,"ExampleNLPDirect::calc_point(...), Error, py is not compatible" );
00159   TEST_FOR_EXCEPTION(
00160     rGf && !this->space_x()->sub_space(this->var_dep())->is_compatible(rGf->space()), std::invalid_argument
00161     ,"ExampleNLPDirect::calc_point(...), Error, py is not compatible" );
00162   TEST_FOR_EXCEPTION(
00163     GcU, std::invalid_argument
00164     ,"ExampleNLPDirect::calc_point(...), Error, there are no undecomposed equalities" );
00165   TEST_FOR_EXCEPTION(
00166     D && !dynamic_cast<MatrixSymDiagStd*>(D), std::invalid_argument
00167     ,"ExampleNLPDirect::calc_point(...), Error, D is not compatible" );
00168   TEST_FOR_EXCEPTION(
00169     py!=NULL && c==NULL, std::invalid_argument
00170     ,"ExampleNLPDirect::calc_point(...) : "
00171     "Error, if py!=NULL then c!=NULL must also be true" );
00172 #endif
00173 
00174   // ///////////////////////////////////
00175   // Compute f(x), c(x) and Gf(x)
00176 
00177   typedef ExampleNLPDirect  this_t;
00178 
00179   // Make temp Gf if needed
00180   VectorSpace::vec_mut_ptr_t  Gf_ptr;
00181   if( rGf && !Gf ) {
00182     Gf_ptr = this->space_x()->create_member();
00183     Gf = Gf_ptr.get();
00184   }
00185 
00186   // Make temp D if needed
00187   mat_fcty_ptr_t::element_type::obj_ptr_t  D_ptr;
00188   if( rGf && !D ) {
00189     D_ptr = this->factory_D()->create();
00190     D = D_ptr.get();
00191   }
00192 
00193   // Remember what references are already set
00194   value_type             *f_saved  = const_cast<this_t*>(this)->get_f();
00195   VectorMutable    *c_saved  = const_cast<this_t*>(this)->get_c();
00196   VectorMutable    *Gf_saved = const_cast<this_t*>(this)->get_Gf();
00197   // Set and compute the quantities
00198   const_cast<this_t*>(this)->set_f(f);
00199   const_cast<this_t*>(this)->set_c(c);
00200   const_cast<this_t*>(this)->set_Gf(Gf);
00201   if(Gf)
00202     this->calc_Gf(x,true);
00203   if(c && recalc_c)
00204     this->calc_c(x,false);
00205   if(f)
00206     this->calc_f(x,false);
00207   // Reset the remembered references
00208   const_cast<this_t*>(this)->set_f(f_saved);
00209   const_cast<this_t*>(this)->set_c(c_saved);
00210   const_cast<this_t*>(this)->set_Gf(Gf_saved);
00211 
00212   // ////////////////////////////////////////////////////////////////////////
00213   // Compute py = -inv(C)*c and/or D = -inv(C)*N at the same time
00214   // 
00215   //        [ 1/(1-x(m+1))                      ]
00216   //        [       1/(1-x(m+2))              ]
00217   //  -inv(C) =   [               .           ]
00218   //        [                 .         ]
00219   //        [                   1/(1-x(m+m))  ]
00220   //
00221   //
00222   //        [ x(1) - 10                       ]
00223   //        [       x(2) - 10               ]
00224   //  N     =   [               .           ]
00225   //        [                 .         ]
00226   //        [                   x(m) - 10   ]
00227 
00228   MatrixSymDiagStd
00229     *D_diag = dynamic_cast<MatrixSymDiagStd*>(D);
00230 
00231   Vector::vec_ptr_t
00232     xD= x.sub_view(this->var_dep()),
00233     xI = x.sub_view(this->var_indep());
00234 
00235   int task;
00236   if     ( py  && !D )  task = 0;
00237   else if( !py &&  D )  task = 1;
00238   else if( py  &&  D )  task = 2;
00239   
00240   TEST_FOR_EXCEPT(0!=RTOp_TOp_explnlp2_calc_py_D_set_task(task,&explnlp2_calc_py_D_op.op()));
00241 
00242   const int              num_vecs = task < 2 ? 2 : 3;
00243   const Vector*          vecs[3] = { NULL, NULL, NULL };
00244   const int              num_targ_vecs = task < 2 ? 1 : 2;
00245   VectorMutable*         targ_vecs[2] = { NULL, NULL };
00246 
00247   // targ_vecs[0] will have apply_op(...) called on it.
00248   if(D) {
00249     D_diag->init_identity( *this->space_c(), 0.0 );
00250     targ_vecs[0]= &D_diag->diag();
00251   }
00252   else if(py)
00253     targ_vecs[0] = py;
00254   else
00255     TEST_FOR_EXCEPT(true); // Only local error?
00256   // targ_vecs[1] will be passed to apply_op(...)
00257   if(py && D)
00258     targ_vecs[1] = py;
00259   
00260   // vecs[...]
00261   int k = 0;
00262   vecs[k] = xD.get(); ++k;
00263   if(D)  { vecs[k] = xI.get(); ++k; }
00264   if(py) { vecs[k] = c;        ++k; }
00265 
00266   AbstractLinAlgPack::apply_op(
00267     explnlp2_calc_py_D_op, num_vecs, vecs, num_targ_vecs, num_targ_vecs?targ_vecs:NULL
00268     ,NULL
00269     );
00270 
00271   // rGf = Gf(var_indep) + D' * Gf(var_dep)
00272   if(rGf) {
00273     *rGf = *Gf->sub_view(this->var_indep());
00274     Vp_MtV( rGf, *D, BLAS_Cpp::trans, *Gf->sub_view(this->var_dep()) );
00275   }
00276 }
00277 
00278 void ExampleNLPDirect::calc_semi_newton_step(
00279   const Vector    &x
00280   ,VectorMutable  *c
00281   ,bool           recalc_c
00282   ,VectorMutable  *py
00283   ) const
00284 {
00285   // In this case just call calc_point(...).
00286   // In a more specialized application, this could be much cheaper!
00287   calc_point(x,NULL,c,recalc_c,NULL,py,NULL,NULL,NULL,NULL);
00288 }
00289 
00290 } // end namespace NLPInterfacePack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines