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 // Redistribution and use in source and binary forms, with or without
00011 // modification, are permitted provided that the following conditions are
00012 // met:
00013 //
00014 // 1. Redistributions of source code must retain the above copyright
00015 // notice, this list of conditions and the following disclaimer.
00016 //
00017 // 2. Redistributions in binary form must reproduce the above copyright
00018 // notice, this list of conditions and the following disclaimer in the
00019 // documentation and/or other materials provided with the distribution.
00020 //
00021 // 3. Neither the name of the Corporation nor the names of the
00022 // contributors may be used to endorse or promote products derived from
00023 // this software without specific prior written permission.
00024 //
00025 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
00026 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
00027 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
00028 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
00029 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
00030 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
00031 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
00032 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
00033 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
00034 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
00035 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
00036 //
00037 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov) 
00038 // 
00039 // ***********************************************************************
00040 // @HEADER
00041 
00042 #include <assert.h>
00043 
00044 #include <stdexcept>
00045 
00046 #include "NLPInterfacePack_ExampleNLPDirect.hpp"
00047 #include "ExampleNLPDirectRTOps.h"
00048 #include "AbstractLinAlgPack_BasisSystemComposite.hpp"
00049 #include "AbstractLinAlgPack_VectorMutable.hpp"
00050 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
00051 #include "AbstractLinAlgPack_VectorStdOps.hpp"
00052 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
00053 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
00054 #include "RTOpPack_RTOpC.hpp"
00055 #include "Teuchos_dyn_cast.hpp"
00056 #include "Teuchos_Assert.hpp"
00057 #include "Teuchos_AbstractFactoryStd.hpp"
00058 
00059 namespace {
00060 
00061 static RTOpPack::RTOpC explnlp2_calc_py_D_op;
00062 
00063 class init_rtop_server_t {
00064 public:
00065   init_rtop_server_t() {
00066     TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_explnlp2_calc_py_D_construct(0,&explnlp2_calc_py_D_op.op()));
00067   }
00068 }; 
00069 init_rtop_server_t  init_rtop_server;
00070 
00071 } // end namespace
00072 
00073 namespace NLPInterfacePack {
00074 
00075 ExampleNLPDirect::ExampleNLPDirect(
00076   const VectorSpace::space_ptr_t&  vec_space
00077   ,value_type                      xo
00078   ,bool                            has_bounds
00079   ,bool                            dep_bounded
00080   )
00081   :ExampleNLPObjGrad(vec_space,xo,has_bounds,dep_bounded),initialized_(false)
00082 {
00083   namespace rcp = MemMngPack;
00084 
00085   // Create the factory object for D
00086   factory_D_ = Teuchos::rcp(new Teuchos::AbstractFactoryStd<MatrixOp,MatrixSymDiagStd>());
00087   NLPDirect::set_factories(
00088     Teuchos::rcp(
00089       new Teuchos::AbstractFactoryStd<MatrixSymOp,MatrixSymDiagStd>())               // D'*D
00090     ,Teuchos::rcp(
00091       new Teuchos::AbstractFactoryStd<MatrixSymOpNonsing,MatrixSymDiagStd>())    // S
00092     );
00093 }
00094 
00095 // Overridden public members from NLP
00096 
00097 void ExampleNLPDirect::initialize(bool test_setup)
00098 {
00099 
00100   if( initialized_ ) {
00101     NLPDirect::initialize(test_setup);
00102     ExampleNLPObjGrad::initialize(test_setup);
00103     return;
00104   }
00105 
00106   NLPDirect::initialize(test_setup);
00107   ExampleNLPObjGrad::initialize(test_setup);
00108 
00109   initialized_ = true;
00110 }
00111 
00112 bool ExampleNLPDirect::is_initialized() const
00113 {
00114   return initialized_;
00115 }
00116 
00117 // Overridden public members from NLPDirect
00118 
00119 Range1D ExampleNLPDirect::var_dep() const
00120 {
00121   return ExampleNLPObjGrad::var_dep();
00122 }
00123 
00124 Range1D ExampleNLPDirect::var_indep() const
00125 {
00126   return ExampleNLPObjGrad::var_indep();
00127 }
00128 
00129 const NLPDirect::mat_fcty_ptr_t
00130 ExampleNLPDirect::factory_D() const
00131 {
00132   return factory_D_;
00133 }
00134 
00135 void ExampleNLPDirect::calc_point(
00136   const Vector     &x
00137   ,value_type      *f
00138   ,VectorMutable   *c
00139   ,bool            recalc_c
00140   ,VectorMutable   *Gf
00141   ,VectorMutable   *py
00142   ,VectorMutable   *rGf
00143   ,MatrixOp        *GcU
00144   ,MatrixOp        *D
00145   ,MatrixOp        *Uz
00146   ) const
00147 {
00148   using Teuchos::dyn_cast;
00149   using LinAlgOpPack::Vp_MtV;
00150 
00151   assert_is_initialized();
00152 
00153   const size_type
00154     n = this->n();
00155 
00156   // Validate the input
00157 
00158 #ifdef TEUCHOS_DEBUG
00159   TEUCHOS_TEST_FOR_EXCEPTION(
00160     x.dim() != n, std::invalid_argument
00161     ,"ExampleNLPDirect::calc_point(...), Error x.dim() = " << x.dim()
00162     << " != n = " << n );
00163   TEUCHOS_TEST_FOR_EXCEPTION(
00164     c && !this->space_c()->is_compatible(c->space()), std::invalid_argument
00165     ,"ExampleNLPDirect::calc_point(...), Error c is not compatible" );
00166   TEUCHOS_TEST_FOR_EXCEPTION(
00167     Gf && !this->space_x()->is_compatible(Gf->space()), std::invalid_argument
00168     ,"ExampleNLPDirect::calc_point(...), Error, Gf is not compatible" );
00169   TEUCHOS_TEST_FOR_EXCEPTION(
00170     py && !this->space_x()->sub_space(this->var_dep())->is_compatible(py->space()), std::invalid_argument
00171     ,"ExampleNLPDirect::calc_point(...), Error, py is not compatible" );
00172   TEUCHOS_TEST_FOR_EXCEPTION(
00173     rGf && !this->space_x()->sub_space(this->var_dep())->is_compatible(rGf->space()), std::invalid_argument
00174     ,"ExampleNLPDirect::calc_point(...), Error, py is not compatible" );
00175   TEUCHOS_TEST_FOR_EXCEPTION(
00176     GcU, std::invalid_argument
00177     ,"ExampleNLPDirect::calc_point(...), Error, there are no undecomposed equalities" );
00178   TEUCHOS_TEST_FOR_EXCEPTION(
00179     D && !dynamic_cast<MatrixSymDiagStd*>(D), std::invalid_argument
00180     ,"ExampleNLPDirect::calc_point(...), Error, D is not compatible" );
00181   TEUCHOS_TEST_FOR_EXCEPTION(
00182     py!=NULL && c==NULL, std::invalid_argument
00183     ,"ExampleNLPDirect::calc_point(...) : "
00184     "Error, if py!=NULL then c!=NULL must also be true" );
00185 #endif
00186 
00187   // ///////////////////////////////////
00188   // Compute f(x), c(x) and Gf(x)
00189 
00190   typedef ExampleNLPDirect  this_t;
00191 
00192   // Make temp Gf if needed
00193   VectorSpace::vec_mut_ptr_t  Gf_ptr;
00194   if( rGf && !Gf ) {
00195     Gf_ptr = this->space_x()->create_member();
00196     Gf = Gf_ptr.get();
00197   }
00198 
00199   // Make temp D if needed
00200   mat_fcty_ptr_t::element_type::obj_ptr_t  D_ptr;
00201   if( rGf && !D ) {
00202     D_ptr = this->factory_D()->create();
00203     D = D_ptr.get();
00204   }
00205 
00206   // Remember what references are already set
00207   value_type             *f_saved  = const_cast<this_t*>(this)->get_f();
00208   VectorMutable    *c_saved  = const_cast<this_t*>(this)->get_c();
00209   VectorMutable    *Gf_saved = const_cast<this_t*>(this)->get_Gf();
00210   // Set and compute the quantities
00211   const_cast<this_t*>(this)->set_f(f);
00212   const_cast<this_t*>(this)->set_c(c);
00213   const_cast<this_t*>(this)->set_Gf(Gf);
00214   if(Gf)
00215     this->calc_Gf(x,true);
00216   if(c && recalc_c)
00217     this->calc_c(x,false);
00218   if(f)
00219     this->calc_f(x,false);
00220   // Reset the remembered references
00221   const_cast<this_t*>(this)->set_f(f_saved);
00222   const_cast<this_t*>(this)->set_c(c_saved);
00223   const_cast<this_t*>(this)->set_Gf(Gf_saved);
00224 
00225   // ////////////////////////////////////////////////////////////////////////
00226   // Compute py = -inv(C)*c and/or D = -inv(C)*N at the same time
00227   // 
00228   //        [ 1/(1-x(m+1))                      ]
00229   //        [       1/(1-x(m+2))              ]
00230   //  -inv(C) =   [               .           ]
00231   //        [                 .         ]
00232   //        [                   1/(1-x(m+m))  ]
00233   //
00234   //
00235   //        [ x(1) - 10                       ]
00236   //        [       x(2) - 10               ]
00237   //  N     =   [               .           ]
00238   //        [                 .         ]
00239   //        [                   x(m) - 10   ]
00240 
00241   MatrixSymDiagStd
00242     *D_diag = dynamic_cast<MatrixSymDiagStd*>(D);
00243 
00244   Vector::vec_ptr_t
00245     xD= x.sub_view(this->var_dep()),
00246     xI = x.sub_view(this->var_indep());
00247 
00248   int task;
00249   if     ( py  && !D )  task = 0;
00250   else if( !py &&  D )  task = 1;
00251   else if( py  &&  D )  task = 2;
00252   
00253   TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_explnlp2_calc_py_D_set_task(task,&explnlp2_calc_py_D_op.op()));
00254 
00255   const int              num_vecs = task < 2 ? 2 : 3;
00256   const Vector*          vecs[3] = { NULL, NULL, NULL };
00257   const int              num_targ_vecs = task < 2 ? 1 : 2;
00258   VectorMutable*         targ_vecs[2] = { NULL, NULL };
00259 
00260   // targ_vecs[0] will have apply_op(...) called on it.
00261   if(D) {
00262     D_diag->init_identity( *this->space_c(), 0.0 );
00263     targ_vecs[0]= &D_diag->diag();
00264   }
00265   else if(py)
00266     targ_vecs[0] = py;
00267   else
00268     TEUCHOS_TEST_FOR_EXCEPT(true); // Only local error?
00269   // targ_vecs[1] will be passed to apply_op(...)
00270   if(py && D)
00271     targ_vecs[1] = py;
00272   
00273   // vecs[...]
00274   int k = 0;
00275   vecs[k] = xD.get(); ++k;
00276   if(D)  { vecs[k] = xI.get(); ++k; }
00277   if(py) { vecs[k] = c;        ++k; }
00278 
00279   AbstractLinAlgPack::apply_op(
00280     explnlp2_calc_py_D_op, num_vecs, vecs, num_targ_vecs, num_targ_vecs?targ_vecs:NULL
00281     ,NULL
00282     );
00283 
00284   // rGf = Gf(var_indep) + D' * Gf(var_dep)
00285   if(rGf) {
00286     *rGf = *Gf->sub_view(this->var_indep());
00287     Vp_MtV( rGf, *D, BLAS_Cpp::trans, *Gf->sub_view(this->var_dep()) );
00288   }
00289 }
00290 
00291 void ExampleNLPDirect::calc_semi_newton_step(
00292   const Vector    &x
00293   ,VectorMutable  *c
00294   ,bool           recalc_c
00295   ,VectorMutable  *py
00296   ) const
00297 {
00298   // In this case just call calc_point(...).
00299   // In a more specialized application, this could be much cheaper!
00300   calc_point(x,NULL,c,recalc_c,NULL,py,NULL,NULL,NULL,NULL);
00301 }
00302 
00303 } // end namespace NLPInterfacePack
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Defines