MOOCHO  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
NLPInterfacePack_ExampleNLPDirect.cpp
1 // @HEADER
2 // ***********************************************************************
3 //
4 // Moocho: Multi-functional Object-Oriented arCHitecture for Optimization
5 // Copyright (2003) Sandia Corporation
6 //
7 // Under terms of Contract DE-AC04-94AL85000, there is a non-exclusive
8 // license for use of this work by or on behalf of the U.S. Government.
9 //
10 // Redistribution and use in source and binary forms, with or without
11 // modification, are permitted provided that the following conditions are
12 // met:
13 //
14 // 1. Redistributions of source code must retain the above copyright
15 // notice, this list of conditions and the following disclaimer.
16 //
17 // 2. Redistributions in binary form must reproduce the above copyright
18 // notice, this list of conditions and the following disclaimer in the
19 // documentation and/or other materials provided with the distribution.
20 //
21 // 3. Neither the name of the Corporation nor the names of the
22 // contributors may be used to endorse or promote products derived from
23 // this software without specific prior written permission.
24 //
25 // THIS SOFTWARE IS PROVIDED BY SANDIA CORPORATION "AS IS" AND ANY
26 // EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE
27 // IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
28 // PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL SANDIA CORPORATION OR THE
29 // CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
30 // EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
31 // PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
32 // PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
33 // LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
34 // NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
35 // SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
36 //
37 // Questions? Contact Roscoe A. Bartlett (rabartl@sandia.gov)
38 //
39 // ***********************************************************************
40 // @HEADER
41 
42 #include <assert.h>
43 
44 #include <stdexcept>
45 
46 #include "NLPInterfacePack_ExampleNLPDirect.hpp"
47 #include "ExampleNLPDirectRTOps.h"
48 #include "AbstractLinAlgPack_BasisSystemComposite.hpp"
49 #include "AbstractLinAlgPack_VectorMutable.hpp"
50 #include "AbstractLinAlgPack_MatrixSymDiagStd.hpp"
51 #include "AbstractLinAlgPack_VectorStdOps.hpp"
52 #include "AbstractLinAlgPack_VectorAuxiliaryOps.hpp"
53 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
54 #include "RTOpPack_RTOpC.hpp"
55 #include "Teuchos_dyn_cast.hpp"
56 #include "Teuchos_Assert.hpp"
57 #include "Teuchos_AbstractFactoryStd.hpp"
58 
59 namespace {
60 
61 static RTOpPack::RTOpC explnlp2_calc_py_D_op;
62 
63 class init_rtop_server_t {
64 public:
65  init_rtop_server_t() {
67  }
68 };
69 init_rtop_server_t init_rtop_server;
70 
71 } // end namespace
72 
73 namespace NLPInterfacePack {
74 
76  const VectorSpace::space_ptr_t& vec_space
77  ,value_type xo
78  ,bool has_bounds
79  ,bool dep_bounded
80  )
81  :ExampleNLPObjGrad(vec_space,xo,has_bounds,dep_bounded),initialized_(false)
82 {
83  namespace rcp = MemMngPack;
84 
85  // Create the factory object for D
90  ,Teuchos::rcp(
92  );
93 }
94 
95 // Overridden public members from NLP
96 
97 void ExampleNLPDirect::initialize(bool test_setup)
98 {
99 
100  if( initialized_ ) {
101  NLPDirect::initialize(test_setup);
102  ExampleNLPObjGrad::initialize(test_setup);
103  return;
104  }
105 
106  NLPDirect::initialize(test_setup);
107  ExampleNLPObjGrad::initialize(test_setup);
108 
109  initialized_ = true;
110 }
111 
113 {
114  return initialized_;
115 }
116 
117 // Overridden public members from NLPDirect
118 
120 {
122 }
123 
125 {
127 }
128 
131 {
132  return factory_D_;
133 }
134 
136  const Vector &x
137  ,value_type *f
138  ,VectorMutable *c
139  ,bool recalc_c
140  ,VectorMutable *Gf
141  ,VectorMutable *py
142  ,VectorMutable *rGf
143  ,MatrixOp *GcU
144  ,MatrixOp *D
145  ,MatrixOp *Uz
146  ) const
147 {
148  using Teuchos::dyn_cast;
149  using LinAlgOpPack::Vp_MtV;
150 
151  assert_is_initialized();
152 
153  const size_type
154  n = this->n();
155 
156  // Validate the input
157 
158 #ifdef TEUCHOS_DEBUG
160  x.dim() != n, std::invalid_argument
161  ,"ExampleNLPDirect::calc_point(...), Error x.dim() = " << x.dim()
162  << " != n = " << n );
164  c && !this->space_c()->is_compatible(c->space()), std::invalid_argument
165  ,"ExampleNLPDirect::calc_point(...), Error c is not compatible" );
166  TEUCHOS_TEST_FOR_EXCEPTION(
167  Gf && !this->space_x()->is_compatible(Gf->space()), std::invalid_argument
168  ,"ExampleNLPDirect::calc_point(...), Error, Gf is not compatible" );
169  TEUCHOS_TEST_FOR_EXCEPTION(
170  py && !this->space_x()->sub_space(this->var_dep())->is_compatible(py->space()), std::invalid_argument
171  ,"ExampleNLPDirect::calc_point(...), Error, py is not compatible" );
172  TEUCHOS_TEST_FOR_EXCEPTION(
173  rGf && !this->space_x()->sub_space(this->var_dep())->is_compatible(rGf->space()), std::invalid_argument
174  ,"ExampleNLPDirect::calc_point(...), Error, py is not compatible" );
175  TEUCHOS_TEST_FOR_EXCEPTION(
176  GcU, std::invalid_argument
177  ,"ExampleNLPDirect::calc_point(...), Error, there are no undecomposed equalities" );
178  TEUCHOS_TEST_FOR_EXCEPTION(
179  D && !dynamic_cast<MatrixSymDiagStd*>(D), std::invalid_argument
180  ,"ExampleNLPDirect::calc_point(...), Error, D is not compatible" );
181  TEUCHOS_TEST_FOR_EXCEPTION(
182  py!=NULL && c==NULL, std::invalid_argument
183  ,"ExampleNLPDirect::calc_point(...) : "
184  "Error, if py!=NULL then c!=NULL must also be true" );
185 #endif
186 
187  // ///////////////////////////////////
188  // Compute f(x), c(x) and Gf(x)
189 
190  typedef ExampleNLPDirect this_t;
191 
192  // Make temp Gf if needed
193  VectorSpace::vec_mut_ptr_t Gf_ptr;
194  if( rGf && !Gf ) {
195  Gf_ptr = this->space_x()->create_member();
196  Gf = Gf_ptr.get();
197  }
198 
199  // Make temp D if needed
200  mat_fcty_ptr_t::element_type::obj_ptr_t D_ptr;
201  if( rGf && !D ) {
202  D_ptr = this->factory_D()->create();
203  D = D_ptr.get();
204  }
205 
206  // Remember what references are already set
207  value_type *f_saved = const_cast<this_t*>(this)->get_f();
208  VectorMutable *c_saved = const_cast<this_t*>(this)->get_c();
209  VectorMutable *Gf_saved = const_cast<this_t*>(this)->get_Gf();
210  // Set and compute the quantities
211  const_cast<this_t*>(this)->set_f(f);
212  const_cast<this_t*>(this)->set_c(c);
213  const_cast<this_t*>(this)->set_Gf(Gf);
214  if(Gf)
215  this->calc_Gf(x,true);
216  if(c && recalc_c)
217  this->calc_c(x,false);
218  if(f)
219  this->calc_f(x,false);
220  // Reset the remembered references
221  const_cast<this_t*>(this)->set_f(f_saved);
222  const_cast<this_t*>(this)->set_c(c_saved);
223  const_cast<this_t*>(this)->set_Gf(Gf_saved);
224 
225  // ////////////////////////////////////////////////////////////////////////
226  // Compute py = -inv(C)*c and/or D = -inv(C)*N at the same time
227  //
228  // [ 1/(1-x(m+1)) ]
229  // [ 1/(1-x(m+2)) ]
230  // -inv(C) = [ . ]
231  // [ . ]
232  // [ 1/(1-x(m+m)) ]
233  //
234  //
235  // [ x(1) - 10 ]
236  // [ x(2) - 10 ]
237  // N = [ . ]
238  // [ . ]
239  // [ x(m) - 10 ]
240 
241  MatrixSymDiagStd
242  *D_diag = dynamic_cast<MatrixSymDiagStd*>(D);
243 
244  Vector::vec_ptr_t
245  xD= x.sub_view(this->var_dep()),
246  xI = x.sub_view(this->var_indep());
247 
248  int task;
249  if ( py && !D ) task = 0;
250  else if( !py && D ) task = 1;
251  else if( py && D ) task = 2;
252 
253  TEUCHOS_TEST_FOR_EXCEPT(0!=RTOp_TOp_explnlp2_calc_py_D_set_task(task,&explnlp2_calc_py_D_op.op()));
254 
255  const int num_vecs = task < 2 ? 2 : 3;
256  const Vector* vecs[3] = { NULL, NULL, NULL };
257  const int num_targ_vecs = task < 2 ? 1 : 2;
258  VectorMutable* targ_vecs[2] = { NULL, NULL };
259 
260  // targ_vecs[0] will have apply_op(...) called on it.
261  if(D) {
262  D_diag->init_identity( *this->space_c(), 0.0 );
263  targ_vecs[0]= &D_diag->diag();
264  }
265  else if(py)
266  targ_vecs[0] = py;
267  else
268  TEUCHOS_TEST_FOR_EXCEPT(true); // Only local error?
269  // targ_vecs[1] will be passed to apply_op(...)
270  if(py && D)
271  targ_vecs[1] = py;
272 
273  // vecs[...]
274  int k = 0;
275  vecs[k] = xD.get(); ++k;
276  if(D) { vecs[k] = xI.get(); ++k; }
277  if(py) { vecs[k] = c; ++k; }
278 
279  AbstractLinAlgPack::apply_op(
280  explnlp2_calc_py_D_op, num_vecs, vecs, num_targ_vecs, num_targ_vecs?targ_vecs:NULL
281  ,NULL
282  );
283 
284  // rGf = Gf(var_indep) + D' * Gf(var_dep)
285  if(rGf) {
286  *rGf = *Gf->sub_view(this->var_indep());
287  Vp_MtV( rGf, *D, BLAS_Cpp::trans, *Gf->sub_view(this->var_dep()) );
288  }
289 }
290 
292  const Vector &x
293  ,VectorMutable *c
294  ,bool recalc_c
295  ,VectorMutable *py
296  ) const
297 {
298  // In this case just call calc_point(...).
299  // In a more specialized application, this could be much cheaper!
300  calc_point(x,NULL,c,recalc_c,NULL,py,NULL,NULL,NULL,NULL);
301 }
302 
303 } // end namespace NLPInterfacePack
void initialize(bool test_setup)
virtual void calc_c(const Vector &x, bool newx=true) const
void calc_semi_newton_step(const Vector &x, VectorMutable *c, bool recalc_c, VectorMutable *py) const
void set_factories(const mat_sym_fcty_ptr_t &factory_transDtD, const mat_sym_nonsing_fcty_ptr_t &factory_S)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
int RTOp_TOp_explnlp2_calc_py_D_set_task(int task, struct RTOp_RTOp *op)
T_To & dyn_cast(T_From &from)
T * get() const
virtual VectorMutable * get_Gf()
ExampleNLPDirect(const VectorSpace::space_ptr_t &vec_space, value_type xo, bool has_bounds, bool dep_bounded)
Constructor.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
int RTOp_TOp_explnlp2_calc_py_D_construct(int task, struct RTOp_RTOp *op)
virtual void calc_f(const Vector &x, bool newx=true) const
size_t size_type
virtual void set_c(VectorMutable *c)
Simple example NLP subclass to illustrate how to implement the NLPObjGrad interface for a specialized...
virtual value_type * get_f()
RTOp_RTOp & op()
Simple example NLP subclass to illustrate how to implement the NLPDirect interface for a specialized ...
virtual VectorMutable * get_c()
virtual void set_Gf(VectorMutable *Gf)
virtual void set_f(value_type *f)
virtual void calc_Gf(const Vector &x, bool newx=true) const
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void calc_point(const Vector &x, value_type *f, VectorMutable *c, bool recalc_c, VectorMutable *Gf, VectorMutable *py, VectorMutable *rGf, MatrixOp *GcU, MatrixOp *D, MatrixOp *Uz) const
virtual VectorMutable & c()