NLPInterfacePack_ExampleNLPDirect.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 <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     m = n/2;
00143 
00144   // Validate the input
00145 
00146 #ifdef TEUCHOS_DEBUG
00147   TEST_FOR_EXCEPTION(
00148     x.dim() != n, std::invalid_argument
00149     ,"ExampleNLPDirect::calc_point(...), Error x.dim() = " << x.dim()
00150     << " != n = " << n );
00151   TEST_FOR_EXCEPTION(
00152     c && !this->space_c()->is_compatible(c->space()), std::invalid_argument
00153     ,"ExampleNLPDirect::calc_point(...), Error c is not compatible" );
00154   TEST_FOR_EXCEPTION(
00155     Gf && !this->space_x()->is_compatible(Gf->space()), std::invalid_argument
00156     ,"ExampleNLPDirect::calc_point(...), Error, Gf is not compatible" );
00157   TEST_FOR_EXCEPTION(
00158     py && !this->space_x()->sub_space(this->var_dep())->is_compatible(py->space()), std::invalid_argument
00159     ,"ExampleNLPDirect::calc_point(...), Error, py is not compatible" );
00160   TEST_FOR_EXCEPTION(
00161     rGf && !this->space_x()->sub_space(this->var_dep())->is_compatible(rGf->space()), std::invalid_argument
00162     ,"ExampleNLPDirect::calc_point(...), Error, py is not compatible" );
00163   TEST_FOR_EXCEPTION(
00164     GcU, std::invalid_argument
00165     ,"ExampleNLPDirect::calc_point(...), Error, there are no undecomposed equalities" );
00166   TEST_FOR_EXCEPTION(
00167     D && !dynamic_cast<MatrixSymDiagStd*>(D), std::invalid_argument
00168     ,"ExampleNLPDirect::calc_point(...), Error, D is not compatible" );
00169   TEST_FOR_EXCEPTION(
00170     py!=NULL && c==NULL, std::invalid_argument
00171     ,"ExampleNLPDirect::calc_point(...) : "
00172     "Error, if py!=NULL then c!=NULL must also be true" );
00173 #endif
00174 
00175   // ///////////////////////////////////
00176   // Compute f(x), c(x) and Gf(x)
00177 
00178   typedef ExampleNLPDirect  this_t;
00179 
00180   // Make temp Gf if needed
00181   VectorSpace::vec_mut_ptr_t  Gf_ptr;
00182   if( rGf && !Gf ) {
00183     Gf_ptr = this->space_x()->create_member();
00184     Gf = Gf_ptr.get();
00185   }
00186 
00187   // Make temp D if needed
00188   mat_fcty_ptr_t::element_type::obj_ptr_t  D_ptr;
00189   if( rGf && !D ) {
00190     D_ptr = this->factory_D()->create();
00191     D = D_ptr.get();
00192   }
00193 
00194   // Remember what references are already set
00195   value_type             *f_saved  = const_cast<this_t*>(this)->get_f();
00196   VectorMutable    *c_saved  = const_cast<this_t*>(this)->get_c();
00197   VectorMutable    *Gf_saved = const_cast<this_t*>(this)->get_Gf();
00198   // Set and compute the quantities
00199   const_cast<this_t*>(this)->set_f(f);
00200   const_cast<this_t*>(this)->set_c(c);
00201   const_cast<this_t*>(this)->set_Gf(Gf);
00202   if(Gf)
00203     this->calc_Gf(x,true);
00204   if(c && recalc_c)
00205     this->calc_c(x,false);
00206   if(f)
00207     this->calc_f(x,false);
00208   // Reset the remembered references
00209   const_cast<this_t*>(this)->set_f(f_saved);
00210   const_cast<this_t*>(this)->set_c(c_saved);
00211   const_cast<this_t*>(this)->set_Gf(Gf_saved);
00212 
00213   // ////////////////////////////////////////////////////////////////////////
00214   // Compute py = -inv(C)*c and/or D = -inv(C)*N at the same time
00215   // 
00216   //        [ 1/(1-x(m+1))                      ]
00217   //        [       1/(1-x(m+2))              ]
00218   //  -inv(C) =   [               .           ]
00219   //        [                 .         ]
00220   //        [                   1/(1-x(m+m))  ]
00221   //
00222   //
00223   //        [ x(1) - 10                       ]
00224   //        [       x(2) - 10               ]
00225   //  N     =   [               .           ]
00226   //        [                 .         ]
00227   //        [                   x(m) - 10   ]
00228 
00229   MatrixSymDiagStd
00230     *D_diag = dynamic_cast<MatrixSymDiagStd*>(D);
00231 
00232   Vector::vec_ptr_t
00233     xD= x.sub_view(this->var_dep()),
00234     xI = x.sub_view(this->var_indep());
00235 
00236   int task;
00237   if     ( py  && !D )  task = 0;
00238   else if( !py &&  D )  task = 1;
00239   else if( py  &&  D )  task = 2;
00240   
00241   TEST_FOR_EXCEPT(0!=RTOp_TOp_explnlp2_calc_py_D_set_task(task,&explnlp2_calc_py_D_op.op()));
00242 
00243   const int              num_vecs = task < 2 ? 2 : 3;
00244   const Vector*          vecs[3] = { NULL, NULL, NULL };
00245   const int              num_targ_vecs = task < 2 ? 1 : 2;
00246   VectorMutable*         targ_vecs[2] = { NULL, NULL };
00247 
00248   // targ_vecs[0] will have apply_op(...) called on it.
00249   if(D) {
00250     D_diag->init_identity( *this->space_c(), 0.0 );
00251     targ_vecs[0]= &D_diag->diag();
00252   }
00253   else if(py)
00254     targ_vecs[0] = py;
00255   else
00256     TEST_FOR_EXCEPT(true); // Only local error?
00257   // targ_vecs[1] will be passed to apply_op(...)
00258   if(py && D)
00259     targ_vecs[1] = py;
00260   
00261   // vecs[...]
00262   int k = 0;
00263   vecs[k] = xD.get(); ++k;
00264   if(D)  { vecs[k] = xI.get(); ++k; }
00265   if(py) { vecs[k] = c;        ++k; }
00266 
00267   AbstractLinAlgPack::apply_op(
00268     explnlp2_calc_py_D_op, num_vecs, vecs, num_targ_vecs, num_targ_vecs?targ_vecs:NULL
00269     ,NULL
00270     );
00271 
00272   // rGf = Gf(var_indep) + D' * Gf(var_dep)
00273   if(rGf) {
00274     *rGf = *Gf->sub_view(this->var_indep());
00275     Vp_MtV( rGf, *D, BLAS_Cpp::trans, *Gf->sub_view(this->var_dep()) );
00276   }
00277 }
00278 
00279 void ExampleNLPDirect::calc_semi_newton_step(
00280   const Vector    &x
00281   ,VectorMutable  *c
00282   ,bool           recalc_c
00283   ,VectorMutable  *py
00284   ) const
00285 {
00286   // In this case just call calc_point(...).
00287   // In a more specialized application, this could be much cheaper!
00288   calc_point(x,NULL,c,recalc_c,NULL,py,NULL,NULL,NULL,NULL);
00289 }
00290 
00291 } // end namespace NLPInterfacePack

Generated on Wed May 12 21:38:01 2010 for MOOCHO by  doxygen 1.4.7