MOOCHO (Single Doxygen Collection)  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Groups Pages
NLPInterfacePack_ExampleNLPBanded.cpp
Go to the documentation of this file.
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 <math.h>
43 
47 #include "Teuchos_Assert.hpp"
48 
49 namespace NLPInterfacePack {
50 
51 // Constructors / initializers
52 
54  size_type nD
55  ,size_type nI
56  ,size_type bw
57  ,size_type mU
58  ,size_type mI
59  ,value_type xo
60  ,value_type xDl
61  ,value_type xDu
62  ,value_type xIl
63  ,value_type xIu
64  ,value_type hl
65  ,value_type hu
66  ,bool nlp_selects_basis
67  ,value_type diag_scal
68  ,value_type diag_vary
69  ,bool sym_basis
70  ,value_type f_offset
71  ,value_type co
72  ,bool ignore_constraints
73  )
74  :is_initialized_(false)
75  ,nlp_selects_basis_(nlp_selects_basis)
76  ,basis_selection_was_given_(false)
77  ,has_var_bounds_(false)
78  ,f_offset_(f_offset)
79  ,nD_(nD)
80  ,nI_(nI)
81  ,bw_(bw)
82  ,mU_(mU)
83  ,mI_(mI)
84  ,ignore_constraints_(ignore_constraints)
85  ,diag_scal_(diag_scal)
86  ,diag_vary_(diag_vary)
87  ,fu_( sym_basis ? 3 : 6 )
88 {
89 #ifdef TEUCHOS_DEBUG
90  const char msg_err_head[] = "ExampleNLPBanded::ExampleNLPBanded(...) : Error";
92  nD <= 0, std::invalid_argument
93  ,msg_err_head<<"!" );
95  nI <= 0 || nD < nI, std::invalid_argument
96  ,msg_err_head<<"!" );
98  bw < 1 || nD < bw, std::invalid_argument
99  ,msg_err_head<<"!" );
101  mU < 0, std::invalid_argument
102  ,msg_err_head<<"!" );
104  mI < 0, std::invalid_argument
105  ,msg_err_head<<"!" );
107  mU != 0, std::invalid_argument
108  ,msg_err_head<<", can't handle undecomposed equalities yet!" );
109 #endif
110  Gc_orig_nz_ = nD_ * ( 1 + 2*(bw_-1) + nI_ ); // Overestimate, ToDo: compute the exact value!
111  Gh_orig_nz_ = 2*mI_;
112  //
114  xl_orig_.resize(xinit_orig_.dim());
115  xu_orig_.resize(xinit_orig_.dim());
118  co_orig_.resize(nD_ + mU_);
119  //
120  xinit_orig_ = xo;
121  xl_orig_(1,nD_) = xDl;
122  xu_orig_(1,nD_) = xDu;
123  xl_orig_(nD_+1,nD_+nI_) = xIl;
124  xu_orig_(nD_+1,nD_+nI_) = xIu;
125  if( mI_ ) {
126  hl_orig_ = hl;
127  hu_orig_ = hu;
128  }
129  co_orig_ = co;
130  //
132  if( xDl > -inf || xDu < +inf || xIl > -inf || xIu < +inf )
133  has_var_bounds_ = true;
134  else
135  has_var_bounds_ = false;
136 }
137 
138 // Overridden public members from NLP
139 
140 void ExampleNLPBanded::initialize(bool test_setup)
141 {
142  if(is_initialized_) {
144  return;
145  }
146  // Nothing to initialize?
148  is_initialized_ = true;
149 }
150 
152 {
153  return is_initialized_;
154 }
155 
157 {
158  return +1e+20; // Functions defined everywhere!
159 }
160 
161 // Overridden protected methods from NLPSerialPreprocess
162 
164 {
165  return !is_initialized_;
166 }
167 
169 {
170  return nD_ + nI_;
171 }
172 
174 {
175  if(ignore_constraints_) return 0;
176  return nD_ + mU_;
177 }
178 
180 {
181  if(ignore_constraints_) return 0;
182  return mI_;
183 }
184 
186 {
187  return xinit_orig_();
188 }
189 
191 {
192  return has_var_bounds_;
193 }
194 
196 {
197  return xl_orig_();
198 }
199 
201 {
202  return xu_orig_();
203 }
204 
206 {
207  return hl_orig_();
208 }
209 
211 {
212  return hu_orig_();
213 }
214 
216  const DVectorSlice &x_full
217  ,bool newx
218  ,const ZeroOrderInfoSerial &zero_order_info
219  ) const
220 {
221  inform_new_point(newx);
222  const DVectorSlice x_orig = x_full(1,imp_n_orig());
223  *zero_order_info.f = f_offset_ + ( 1.0 / 2.0 ) * DenseLinAlgPack::dot( x_orig, x_orig );
224 }
225 
227  const DVectorSlice &x_full
228  ,bool newx
229  ,const ZeroOrderInfoSerial &zero_order_info
230  ) const
231 {
232  inform_new_point(newx);
233  if(c_orig_updated_)
234  return; // c(x) is already computed in *zero_order_info.c
235  TEUCHOS_TEST_FOR_EXCEPT( !( zero_order_info.c ) );
236  DVector
237  &c = *zero_order_info.c;
238  const size_type
239  num_I_per_D = nD_ / nI_, // Integer division (rounds down)
240  I_remainder = nD_ % nI_;
241  size_type j = 0;
242  const value_type
243  ds_alpha = nD_ > 1 ? diag_scal_ * (diag_vary_ - 1.0)/(nD_ - 1.0) : 0.0,
244  ds_beta = diag_scal_;
245  for( size_type q_i = 1; q_i <= nI_; ++q_i ) {
246  const size_type num_I_per_D_local = num_I_per_D + ( q_i <= I_remainder ? 1 : 0 );
247  for( size_type q_k = 0; q_k < num_I_per_D_local; ++q_k ) {
248  ++j;
249  const size_type
250  klu = ( j - bw_ >= 0 ? bw_-1 : j-1 ),
251  kuu = ( j + bw_ - 1 <= nD_ ? bw_-1 : nD_-j );
252  const value_type
253  ds_j = ds_alpha * (j-1) + ds_beta;
254  value_type
255  &c_j = (c(j) = ds_j * x_full(j));
256  {for( size_type k = 1; k <= klu; ++k ) {
257  c_j -= (3.0 / k) * x_full(j-k);
258  }}
259  {for( size_type k = 1; k <= kuu; ++k ) {
260  c_j -= (fu_ / k) * x_full(j+k);
261  }}
262  const value_type
263  term = x_full(nD_ + q_i) + 1;
264  c_j *= (term * term);
265  c_j += co_orig_(j);
266  }
267  }
268  c_orig_updated_ = true;
269 }
270 
272  const DVectorSlice &x_full
273  ,bool newx
274  ,const ZeroOrderInfoSerial &zero_order_info
275  ) const
276 {
277  inform_new_point(newx);
278  DVector
279  &h = *zero_order_info.h;
280  const size_type
281  num_I_per_D = nD_ / nI_, // Integer division (rounds down)
282  I_remainder = nD_ % nI_;
283  size_type jI = 0;
284  for( size_type q_i = 1; q_i <= nI_; ++q_i ) {
285  const value_type x_q = x_full(nD_ + q_i);
286  const size_type num_I_per_D_local = num_I_per_D + ( q_i <= I_remainder ? 1 : 0 );
287  for( size_type q_k = 0; q_k < num_I_per_D_local; ++q_k ) {
288  ++jI;
289  if( jI > mI_ ) return;
290  h(jI) = x_full(jI) - x_q;
291  }
292  }
293  }
294 
296  const DVectorSlice &x_full
297  ,bool newx
298  ,const ObjGradInfoSerial &obj_grad_info
299  ) const
300 {
301  inform_new_point(newx);
302  const Range1D var_orig(1,imp_n_orig());
303  (*obj_grad_info.Gf)(var_orig) = x_full(var_orig);
304 }
305 
307  IVector *var_perm_full
308  ,IVector *equ_perm_full
309  ,size_type *rank_full
310  ,size_type *rank
311  )
312 {
314  return false; // Already gave this basis selection.
315  // Select the first nD variables as basis variables which gives
316  // a nice banded matrix for the basis matrix C.
317  // Also, if the general inequality constraints are begin
318  // converted to equalities with slacks, make the slack variables
319  // basic variables also (after the nD variables).
320 #ifdef TEUCHOS_DEBUG
321  TEUCHOS_TEST_FOR_EXCEPT( !( var_perm_full ) );
322  TEUCHOS_TEST_FOR_EXCEPT( !( equ_perm_full ) );
323  TEUCHOS_TEST_FOR_EXCEPT( !( rank ) );
324 #endif
325  const size_type n_orig = nD_ + nI_;
326  const size_type n_full = n_orig + mI_;
327  const size_type m_full = nD_ + mI_;
328  var_perm_full->resize(n_full);
329  equ_perm_full->resize(m_full);
330  if( mI_ ) {
331  index_type k, i_perm = 1;;
332  // basic variables
333  for( k = 1; k <= nD_; ++k, ++i_perm ) // Put xD variables first
334  (*var_perm_full)(i_perm) = k;
335  for( k = 1; k <= mI_; ++k, ++i_perm ) // Put slacks after xD
336  (*var_perm_full)(i_perm) = n_orig + k;
337  // non-basic variables
338  for( k = 1; k <= nI_; ++k, ++i_perm ) // Followed by nI
339  (*var_perm_full)(i_perm) = nD_ + k;
340  }
341  else {
342  DenseLinAlgPack::identity_perm( var_perm_full );
343  }
344  DenseLinAlgPack::identity_perm( equ_perm_full );
345  // Count the number of fixed basic variables to reduce
346  // the rank of the basis.
347  index_type num_fixed = 0;
348  for( index_type k = 1; k <= nD_; ++k ) {
349  if( xl_orig_(k) == xu_orig_(k) )
350  ++num_fixed;
351  }
352  if( mI_ ) {
353  for( index_type k = 1; k <= mI_; ++k ) {
354  if( hl_orig_(k) == hu_orig_(k) )
355  ++num_fixed;
356  }
357  }
358  // Set the rank
359  *rank_full = m_full;
360  *rank = m_full - num_fixed;
362  return true;
363 }
364 
366  const DVectorSlice &x_orig
367  ,const DVectorSlice *lambda_orig
368  ,const DVectorSlice *lambdaI_orig
369  ,const DVectorSlice *nu_orig
370  ,bool is_optimal
371  )
372 {
373  // ToDo: Do something with the final soltuion?
374 }
375 
377 {
378  return nlp_selects_basis_;
379 }
380 
381 // Overridden protected methods from NLPSerialPreprocessExplJac
382 
384 {
385  return Gc_orig_nz_;
386 }
387 
389 {
390  return Gh_orig_nz_;
391 }
392 
394  const DVectorSlice& x_full, bool newx
395  , const FirstOrderExplInfo& first_order_expl_info
396  ) const
397 {
398  inform_new_point(newx);
399  // Compute c(x) if not already (will compute in temp if needed)
400  this->imp_calc_c_orig( x_full, newx, zero_order_orig_info() );
401  TEUCHOS_TEST_FOR_EXCEPT( !( first_order_expl_info.c ) );
402  DVector
403  &c = *first_order_expl_info.c;
404  // Get references/pointers to data for Gc to be computed/updated.
405  index_type
406  &Gc_nz = *first_order_expl_info.Gc_nz;
407  value_type
408  *Gc_val = &(*first_order_expl_info.Gc_val)[0];
409  index_type
410  *Gc_ivect = ( first_order_expl_info.Gc_ivect
411  ? &(*first_order_expl_info.Gc_ivect)[0] : NULL ),
412  *Gc_jvect = ( first_order_expl_info.Gc_jvect
413  ? &(*first_order_expl_info.Gc_jvect)[0] : NULL );
414  TEUCHOS_TEST_FOR_EXCEPT( !( (Gc_ivect != NULL) == (Gc_jvect != NULL) ) );
415  // Set nonzeros for Gc (in sorted compressed column format w.r.t., i.e. grouped by constraints)
416  const size_type
417  num_I_per_D = nD_ / nI_, // Integer division (rounds down)
418  I_remainder = nD_ % nI_;
419  Gc_nz = 0;
420  size_type j = 0;
421  const value_type
422  ds_alpha = nD_ > 1 ? diag_scal_ * (diag_vary_ - 1.0)/(nD_ - 1.0) : 0.0,
423  ds_beta = diag_scal_;
424  for( size_type q_i = 1; q_i <= nI_; ++q_i ) {
425  const value_type
426  x_q = x_full(nD_ + q_i),
427  x_q_term = (x_q + 1) * (x_q + 1);
428  const size_type num_I_per_D_local = num_I_per_D + ( q_i <= I_remainder ? 1 : 0 );
429  for( size_type q_k = 0; q_k < num_I_per_D_local; ++q_k ) {
430  ++j;
431  const size_type
432  klu = ( j - bw_ >= 0 ? bw_-1 : j-1 ),
433  kuu = ( j + bw_ - 1 <= nD_ ? bw_-1 : nD_-j );
434  const value_type
435  ds_j = ds_alpha * (j-1) + ds_beta;
436  //
437  {for( index_type k = klu; k >= 1; --k ) {
438  ++Gc_nz;
439  *Gc_val++ = -3.0 / k * x_q_term;
440  if(Gc_ivect) {
441  *Gc_ivect++ = j - k;
442  *Gc_jvect++ = j;
443  }
444  }}
445  //
446  ++Gc_nz;
447  *Gc_val++ = ds_j * x_q_term;
448  if(Gc_ivect) {
449  *Gc_ivect++ = j;
450  *Gc_jvect++ = j;
451  }
452  //
453  {for( index_type k = 1; k <= kuu; ++k ) {
454  ++Gc_nz;
455  *Gc_val++ = -fu_ / k * x_q_term;
456  if(Gc_ivect) {
457  *Gc_ivect++ = j + k;
458  *Gc_jvect++ = j;
459  }
460  }}
461  //
462  ++Gc_nz;
463  *Gc_val++ = 2.0 * (c(j) - co_orig_(j)) / (x_q + 1);
464  if(Gc_ivect) {
465  *Gc_ivect++ = nD_ + q_i;
466  *Gc_jvect++ = j;
467  }
468  }
469  }
470 }
471 
473  const DVectorSlice& x_full, bool newx
474  , const FirstOrderExplInfo& first_order_expl_info
475  ) const
476 {
477  inform_new_point(newx);
478  // Get references/pointers to data for Gh to be computed/updated.
479  index_type
480  &Gh_nz = *first_order_expl_info.Gh_nz;
481  value_type
482  *Gh_val = &(*first_order_expl_info.Gh_val)[0];
483  index_type
484  *Gh_ivect = ( first_order_expl_info.Gh_ivect
485  ? &(*first_order_expl_info.Gh_ivect)[0] : NULL ),
486  *Gh_jvect = ( first_order_expl_info.Gh_jvect
487  ? &(*first_order_expl_info.Gh_jvect)[0] : NULL );
488  TEUCHOS_TEST_FOR_EXCEPT( !( (Gh_ivect != NULL) == (Gh_jvect != NULL) ) );
489  // Set nonzeros for Gh (in sorted compressed column format w.r.t., i.e. grouped by constraints)
490  const size_type
491  num_I_per_D = nD_ / nI_, // Integer division (rounds down)
492  I_remainder = nD_ % nI_;
493  Gh_nz = 0;
494  size_type jI = 0;
495  for( size_type q_i = 1; q_i <= nI_; ++q_i ) {
496  const size_type nD_q_i = nD_ + q_i;
497  const size_type num_I_per_D_local = num_I_per_D + ( q_i <= I_remainder ? 1 : 0 );
498  for( size_type q_k = 0; q_k < num_I_per_D_local; ++q_k ) {
499  ++jI;
500  if( jI > mI_ ) return;
501  // w.r.t. x(jI)
502  ++Gh_nz;
503  *Gh_val++ = 1.0;
504  if(Gh_ivect) {
505  *Gh_ivect++ = jI;
506  *Gh_jvect++ = jI;
507  }
508  // w.r.t. x(nD+q(jI))
509  ++Gh_nz;
510  *Gh_val++ = -1.0;
511  if(Gh_ivect) {
512  *Gh_ivect++ = nD_q_i;
513  *Gh_jvect++ = jI;
514  }
515  }
516  }
517 }
518 
519 // private
520 
522 {
523  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implemenet!
524 }
525 
527 {
528  if(newx) {
529  c_orig_updated_ = false;
530  }
531 }
532 
533 } // end namespace NLPInterfacePack
AbstractLinAlgPack::size_type size_type
value_type * f
Pointer to objective function f (may be NULL if not set)
DVector * h
Pointer to constraints residual h (may be NULL if not set)
#define TEUCHOS_TEST_FOR_EXCEPTION(throw_exception_test, Exception, msg)
int resize(OrdinalType length_in)
DVectorSlice x_full() const
Give reference to current x_full.
bool imp_get_next_basis(IVector *var_perm_full, IVector *equ_perm_full, size_type *rank_full, size_type *rank)
Struct for objective and constriants (pointer) as serial vectors.
Struct for zero and explicit first order quantities that subclass must fill in.
void imp_calc_h_orig(const DVectorSlice &x_full, bool newx, const ZeroOrderInfoSerial &zero_order_info) const
DVector * Gf
Gradient of objective function Gf (may be NULL if not set)
void imp_calc_Gc_orig(const DVectorSlice &x_full, bool newx, const FirstOrderExplInfo &first_order_expl_info) const
DVector * c
Pointer to constraints residual c (may be NULL if not set)
void identity_perm(IVector *perm)
DenseLinAlgPack::VectorSliceTmpl< value_type > DVectorSlice
const char inf
void imp_calc_c_orig(const DVectorSlice &x_full, bool newx, const ZeroOrderInfoSerial &zero_order_info) const
AbstractLinAlgPack::value_type value_type
Struct for serial gradient (objective), objective and constriants (pointers)
void imp_calc_Gf_orig(const DVectorSlice &x_full, bool newx, const ObjGradInfoSerial &obj_grad_info) const
void imp_calc_f_orig(const DVectorSlice &x_full, bool newx, const ZeroOrderInfoSerial &zero_order_info) const
void imp_report_orig_final_solution(const DVectorSlice &x_orig, const DVectorSlice *lambda_orig, const DVectorSlice *lambdaI_orig, const DVectorSlice *nu_orig, bool is_optimal)
static value_type infinite_bound()
Value for an infinite bound.
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
void imp_calc_Gh_orig(const DVectorSlice &x_full, bool newx, const FirstOrderExplInfo &first_order_expl_info) const
value_type dot(const DVectorSlice &vs_rhs1, const DVectorSlice &vs_rhs2)
result = vs_rhs1' * vs_rhs2 (BLAS xDOT)
virtual VectorMutable & c()
Returns non-const *this->get_c().