ConstrainedOptPack: C++ Tools for Constrained (and Unconstrained) Optimization  Version of the Day
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Groups Pages
ConstrainedOptPack_QPSolverRelaxedQPSchur.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 <vector>
45 
46 #include "ConstrainedOptPack_QPSolverRelaxedQPSchur.hpp"
47 #include "AbstractLinAlgPack_MatrixOp.hpp"
48 #include "AbstractLinAlgPack_LinAlgOpPack.hpp"
49 #include "AbstractLinAlgPack_VectorSpaceFactory.hpp"
50 #include "AbstractLinAlgPack_SortByDescendingAbsValue.hpp"
51 #include "AbstractLinAlgPack_VectorDenseEncap.hpp"
52 #include "AbstractLinAlgPack_VectorSpaceSerial.hpp"
53 #include "AbstractLinAlgPack_sparse_bounds.hpp"
54 #include "DenseLinAlgPack_LinAlgOpPack.hpp"
55 #include "Teuchos_dyn_cast.hpp"
56 #include "ProfileHackPack_profile_hack.hpp"
57 
58 namespace ConstrainedOptPack {
59 
61  const init_kkt_sys_ptr_t& init_kkt_sys
62  ,const constraints_ptr_t& constraints
63  ,value_type max_qp_iter_frac
64  ,value_type max_real_runtime
66  inequality_pick_policy
67  ,ELocalOutputLevel print_level
68  ,value_type bounds_tol
69  ,value_type inequality_tol
70  ,value_type equality_tol
71  ,value_type loose_feas_tol
72  ,value_type dual_infeas_tol
73  ,value_type huge_primal_step
74  ,value_type huge_dual_step
75  ,value_type bigM
76  ,value_type warning_tol
77  ,value_type error_tol
78  ,size_type iter_refine_min_iter
79  ,size_type iter_refine_max_iter
80  ,value_type iter_refine_opt_tol
81  ,value_type iter_refine_feas_tol
82  ,bool iter_refine_at_solution
83  ,value_type pivot_warning_tol
84  ,value_type pivot_singular_tol
85  ,value_type pivot_wrong_inertia_tol
86  ,bool add_equalities_initially
87  )
88  :init_kkt_sys_(init_kkt_sys)
89  ,constraints_(constraints)
90  ,max_qp_iter_frac_(max_qp_iter_frac)
91  ,max_real_runtime_(max_real_runtime)
92  ,inequality_pick_policy_(inequality_pick_policy)
93  ,print_level_(print_level)
94  ,bounds_tol_(bounds_tol)
95  ,inequality_tol_(inequality_tol)
96  ,equality_tol_(equality_tol)
97  ,loose_feas_tol_(loose_feas_tol)
98  ,dual_infeas_tol_(dual_infeas_tol)
99  ,huge_primal_step_(huge_primal_step)
100  ,huge_dual_step_(huge_dual_step)
101  ,bigM_(bigM)
102  ,warning_tol_(warning_tol)
103  ,error_tol_(error_tol)
104  ,iter_refine_min_iter_(iter_refine_min_iter)
105  ,iter_refine_max_iter_(iter_refine_max_iter)
106  ,iter_refine_opt_tol_(iter_refine_opt_tol)
107  ,iter_refine_feas_tol_(iter_refine_feas_tol)
108  ,iter_refine_at_solution_(iter_refine_at_solution)
109  ,pivot_warning_tol_(pivot_warning_tol)
110  ,pivot_singular_tol_(pivot_singular_tol)
111  ,pivot_wrong_inertia_tol_(pivot_wrong_inertia_tol)
112  ,add_equalities_initially_(add_equalities_initially)
113 {}
114 
116 {
117  this->release_memory();
118 }
119 
120 // Overridden from QPSolverRelaxed
121 
124 {
125  return qp_stats_;
126 }
127 
129 {
130  // Nothing to release!
131 }
132 
135  std::ostream* out, EOutputLevel olevel, ERunTests test_what
136  ,const Vector& g, const MatrixSymOp& G
137  ,value_type etaL
138  ,const Vector* dL, const Vector* dU
139  ,const MatrixOp* E, BLAS_Cpp::Transp trans_E, const Vector* b
140  ,const Vector* eL, const Vector* eU
141  ,const MatrixOp* F, BLAS_Cpp::Transp trans_F, const Vector* f
142  ,value_type* obj_d
143  ,value_type* eta, VectorMutable* d
144  ,VectorMutable* nu
145  ,VectorMutable* mu, VectorMutable* Ed
146  ,VectorMutable* lambda, VectorMutable* Fd
147  )
148 {
149  namespace mmp = MemMngPack;
150  using Teuchos::dyn_cast;
151  using LinAlgOpPack::V_mV;
152  typedef QPSchurPack::ConstraintsRelaxedStd constr_t;
153 
154 #ifdef PROFILE_HACK_ENABLED
155  ProfileHackPack::ProfileTiming profile_timing( "QPSolverRelaxedQPSchur::imp_solve_qp(...)" );
156 #endif
157 
158  const size_type
159  nd = g.dim(),
160  m_in = E ? BLAS_Cpp::rows(E->rows(),E->cols(),trans_E) : 0,
161  m_eq = F ? BLAS_Cpp::rows(F->rows(),F->cols(),trans_F) : 0;
162 
163  VectorDenseEncap g_de(g);
164 
165  VectorSpace::space_ptr_t
166  space_d_eta = d->space().small_vec_spc_fcty()->create_vec_spc(nd+1); // ToDo: Generalize!
167 
168  // ///////////////////////////////
169  // Setup the initial KKT system
170 
171  InitKKTSystem::i_x_free_t i_x_free;
172  InitKKTSystem::i_x_fixed_t i_x_fixed;
173  InitKKTSystem::bnd_fixed_t bnd_fixed;
174  InitKKTSystem::j_f_decomp_t j_f_decomp;
175  size_type n_R_tmp;
176  init_kkt_sys().initialize_kkt_system(
177  g,G,etaL,dL,dU,F,trans_F,f,d,nu
178  ,&n_R_tmp,&i_x_free,&i_x_fixed,&bnd_fixed,&j_f_decomp
179  ,&b_X_,&Ko_,&fo_ );
180  const size_type
181  n_R = n_R_tmp,
182  n_X = nd + 1 - n_R; // fixed variables in d and eta
183  TEUCHOS_TEST_FOR_EXCEPT( !( i_x_free.size() == 0 || i_x_free.size() >= n_R ) ); // Todo: Make an exception!
184  TEUCHOS_TEST_FOR_EXCEPT( !( i_x_fixed.size() >= n_X ) ); // Todo: Make an exception!
185  TEUCHOS_TEST_FOR_EXCEPT( !( bnd_fixed.size() >= n_X ) ); // Todo: Make and exception!
186 
187  // //////////////////////////////
188  // Initialize constraints object
189 
190  // Setup j_f_undecomp
191  const bool all_f_undecomp = F ? j_f_decomp.size() == 0 : true;
192  const size_type
193  m_undecomp = F ? f->dim() - j_f_decomp.size() : 0;
194  typedef std::vector<size_type> j_f_undecomp_t;
195  j_f_undecomp_t j_f_undecomp;
196  if( m_undecomp && !all_f_undecomp ) {
197  j_f_undecomp.resize(m_undecomp);
198  // Create a full lookup array to determine if a constraint
199  // is decomposed or not. We need this to fill the array
200  // j_f_undecomp[] (which is sorted).
201  TEUCHOS_TEST_FOR_EXCEPT(true); // ToDo: Implement this!
202  }
203 
204  // initialize constraints object
205  constraints_->initialize(
206  space_d_eta
207  ,etaL,dL,dU,E,trans_E,b,eL,eU,F,trans_F,f
208  ,m_undecomp, m_undecomp && !all_f_undecomp ? &j_f_undecomp[0] : NULL
209  ,Ed
210  ,!add_equalities_initially() // If we add equalities the the schur complement intially
211  // then we don't need to check if they are violated.
212  );
213  // ToDo: Add j_f_decomp to the above constraints class!
214 
215  // ///////////////////////////
216  // Initialize the QP object
217 
218  // g_relaxed_ = [ g; bigM ]
219  g_relaxed_.resize(nd+1);
220  g_relaxed_(1,nd) = g_de();
221  g_relaxed_(nd+1) = bigM();
222 
223  // G_relaxed_ = [ G, zeros(...); zeros(...), bigM ]
224  bigM_vec_.initialize(1); // dim == 1
225  bigM_vec_ = bigM();
226  G_relaxed_.initialize(
227  Teuchos::rcp(&dyn_cast<const MatrixSymOpNonsing>(G),false)
228  ,Teuchos::rcp(&bigM_vec_,false)
229  ,space_d_eta
230  );
231 
232  qp_.initialize(
233  g_relaxed_(),G_relaxed_,NULL
234  ,n_R, i_x_free.size() ? &i_x_free[0] : NULL
235  ,&i_x_fixed[0],&bnd_fixed[0]
236  ,b_X_(),*Ko_,fo_(),constraints_.get()
237  ,out,test_what==RUN_TESTS,warning_tol(),error_tol()
238  ,int(olevel)>=int(PRINT_ITER_VECTORS)
239  );
240 
241  // ///////////////////////////////////////////////////////
242  // Setup for a warm start (changes to initial KKT system)
243 
244  typedef std::vector<int> ij_act_change_t;
245  typedef std::vector<EBounds> bnds_t;
246  size_type num_act_change = 0; // The default is a cold start
247  const size_type max_num_act_change = 2*nd;
248  ij_act_change_t ij_act_change(max_num_act_change);
249  bnds_t bnds(max_num_act_change);
250  // Go ahead and add the equality constraints. If these are linearly
251  // dependent let's hope that QPSchur can handle this and still do a
252  // good job of things. This is a scary thing to do!
253  if( m_eq && add_equalities_initially() ) {
254  for( size_type j = 1; j <= m_eq; ++j ) {
255  ij_act_change[num_act_change] = (nd + 1) + m_in + j;
256  bnds[num_act_change] = EQUALITY;
257  ++num_act_change;
258  }
259  }
260  // We will add fixed (EQUALITY) variable bounds to the initial active set
261  // (if it is not already an intially fixed variable). If fixing a variable
262  // causes the KKT system to become singular then we are in real trouble!
263  // We should add these eairly on since they will not be freed.
264  if( dL ) {
265  const QPSchurPack::QP::x_init_t &x_init = qp_.x_init();
266  const value_type inf_bnd = this->infinite_bound();
267  VectorDenseEncap dL_de(*dL);
268  VectorDenseEncap dU_de(*dU);
269  // read iterators
271  dLU_itr( dL_de().begin(), dL_de().end()
272  ,dU_de().begin(), dU_de().end()
273  ,inf_bnd );
274  for( ; !dLU_itr.at_end(); ++dLU_itr ) {
275  if( dLU_itr.lbound() == dLU_itr.ubound() && x_init(dLU_itr.index()) == FREE ) {
276  ij_act_change[num_act_change] = dLU_itr.index();
277  bnds[num_act_change] = EQUALITY;
278  ++num_act_change;
279  }
280  }
281  }
282  // Add inequality constriants to the list from nu and mu
283  if( ( nu && nu->nz() ) || ( m_in && mu->nz() ) ) {
284  //
285  // Setup num_act_change, ij_act_change, and bnds for a warm start!
286  //
287  const size_type
288  nu_nz = nu ? nu->nz() : 0,
289  mu_nz = mu ? mu->nz() : 0;
290  // Combine all the multipliers for the bound and general inequality
291  // constraints and sort them from the largest to the smallest. Hopefully
292  // the constraints with the larger multiplier values will not be dropped
293  // from the active set.
294  SpVector gamma( nd + 1 + m_in , nu_nz + mu_nz );
295  typedef SpVector::element_type ele_t;
296  if(nu && nu_nz) {
297  VectorDenseEncap nu_de(*nu);
298  DVectorSlice::const_iterator
299  nu_itr = nu_de().begin(),
300  nu_end = nu_de().end();
301  index_type i = 1;
302  while( nu_itr != nu_end ) {
303  for( ; *nu_itr == 0.0; ++nu_itr, ++i );
304  gamma.add_element(ele_t(i,*nu_itr));
305  ++nu_itr; ++i;
306  }
307  }
308  if(mu && mu_nz) {
309  VectorDenseEncap mu_de(*mu);
310  DVectorSlice::const_iterator
311  mu_itr = mu_de().begin(),
312  mu_end = mu_de().end();
313  index_type i = 1;
314  while( mu_itr != mu_end ) {
315  for( ; *mu_itr == 0.0; ++mu_itr, ++i );
316  gamma.add_element(ele_t(i+nd,*mu_itr));
317  ++mu_itr; ++i;
318  }
319  }
320  std::sort( gamma.begin(), gamma.end()
322  // Now add the inequality constraints in decreasing order (if they are
323  // not already initially fixed variables)
324  const QPSchurPack::QP::x_init_t &x_init = qp_.x_init();
325  if(gamma.nz()) {
326  const SpVector::difference_type o = gamma.offset();
327  for( SpVector::const_iterator itr = gamma.begin(); itr != gamma.end(); ++itr ) {
328  const size_type i = itr->index() + o;
329  if( i <= nd && x_init(i) != FREE )
330  continue; // This variable is already initially fixed
331  // This is not an initially fixed variable so add it
332  ij_act_change[num_act_change] = i;
333  bnds[num_act_change]
334  = itr->value() > 0.0 ? UPPER : LOWER;
335  ++num_act_change;
336  }
337  }
338  }
339  // We need to loop through x_init() and nu() in order and look for variables
340  // that are initially fixed in x_init() but are not present in nu(). For these
341  // variables we need to free them in ij_act_change[].
342  if( dL && nu->nz() ) {
343  QPSchurPack::QP::x_init_t::const_iterator
344  x_init_itr = qp_.x_init().begin();
345  VectorDenseEncap nu_de(*nu);
346  DVectorSlice::const_iterator
347  nu_itr = nu_de().begin();
348  for( size_type i = 1; i <= nd; ++i, ++x_init_itr, ++nu_itr ) {
349  if( *x_init_itr != FREE && *x_init_itr != EQUALITY ) {
350  // This is an initially fixed upper or lower bound.
351  // Look for lagrange multiplier stating that it is
352  // still fixed.
353  if( *nu_itr != 0.0 ) {
354  // This active bound is present but lets make sure
355  // that it is still the same bound
356  if( ( *x_init_itr == LOWER && *nu_itr > 0 )
357  || ( *x_init_itr == UPPER && *nu_itr < 0 ) )
358  {
359  // The bound has changed from upper to lower or visa-versa!
360  ij_act_change[num_act_change] = i;
361  bnds[num_act_change]
362  = *nu_itr > 0.0 ? UPPER : LOWER;
363  ++num_act_change;
364  }
365  }
366  else {
367  // This initially fixed variable is not fixed in nu so lets free it!
368  ij_act_change[num_act_change] = -i;
369  bnds[num_act_change] = FREE;
370  ++num_act_change;
371  }
372  }
373  }
374  }
375  // Consider the relaxation variable!
376  if( *eta > etaL) {
377  ij_act_change[num_act_change] = -int(nd+1);
378  bnds[num_act_change] = FREE;
379  ++num_act_change;
380  }
381 
382  // Set the output level
383  QPSchur::EOutputLevel qpschur_olevel;
384  switch( print_level() ) {
385  case USE_INPUT_ARG: {
386  // Use the input print level
387  switch( olevel ) {
388  case PRINT_NONE:
389  qpschur_olevel = QPSchur::NO_OUTPUT;
390  break;
391  case PRINT_BASIC_INFO:
392  qpschur_olevel = QPSchur::OUTPUT_BASIC_INFO;
393  break;
394  case PRINT_ITER_SUMMARY:
395  qpschur_olevel = QPSchur::OUTPUT_ITER_SUMMARY;
396  break;
397  case PRINT_ITER_STEPS:
398  qpschur_olevel = QPSchur::OUTPUT_ITER_STEPS;
399  break;
400  case PRINT_ITER_ACT_SET:
401  case PRINT_ITER_VECTORS:
402  qpschur_olevel = QPSchur::OUTPUT_ACT_SET;
403  break;
404  case PRINT_EVERY_THING:
405  qpschur_olevel = QPSchur::OUTPUT_ITER_QUANTITIES;
406  break;
407  default:
409  }
410  break;
411  }
412  case NO_OUTPUT:
413  qpschur_olevel = QPSchur::NO_OUTPUT;
414  break;
415  case OUTPUT_BASIC_INFO:
416  qpschur_olevel = QPSchur::OUTPUT_BASIC_INFO;
417  break;
418  case OUTPUT_ITER_SUMMARY:
419  qpschur_olevel = QPSchur::OUTPUT_ITER_SUMMARY;
420  break;
421  case OUTPUT_ITER_STEPS:
422  qpschur_olevel = QPSchur::OUTPUT_ITER_STEPS;
423  break;
424  case OUTPUT_ACT_SET:
425  qpschur_olevel = QPSchur::OUTPUT_ACT_SET;
426  break;
427  case OUTPUT_ITER_QUANTITIES:
428  qpschur_olevel = QPSchur::OUTPUT_ITER_QUANTITIES;
429  break;
430  default:
432  }
433 
434  //
435  // Set options for ConstraintsRelaxedStd.
436  //
437  if( bounds_tol() > 0.0 )
438  constraints_->bounds_tol(bounds_tol());
439  if( inequality_tol() > 0.0 )
440  constraints_->inequality_tol(inequality_tol());
441  if( equality_tol() > 0.0 )
442  constraints_->equality_tol(equality_tol());
443  constraints_->inequality_pick_policy(inequality_pick_policy());
444 
445  //
446  // Set options for QPSchur.
447  //
448  qp_solver_.max_iter(static_cast<index_type>(max_qp_iter_frac()*nd) );
449  qp_solver_.max_real_runtime( max_real_runtime() );
450  qp_solver_.feas_tol( constraints_->bounds_tol() ); // Let's assume the bound tolerance is the tightest
451  if(loose_feas_tol() > 0.0)
452  qp_solver_.loose_feas_tol( loose_feas_tol() );
453  else
454  qp_solver_.loose_feas_tol( 10.0 * qp_solver_.feas_tol() );
455  if(dual_infeas_tol() > 0.0)
456  qp_solver_.dual_infeas_tol( dual_infeas_tol() );
457  if(huge_primal_step() > 0.0)
458  qp_solver_.huge_primal_step( huge_primal_step() );
459  if(huge_dual_step() > 0.0)
460  qp_solver_.huge_dual_step( huge_dual_step() );
461  qp_solver_.set_schur_comp( QPSchur::schur_comp_ptr_t( &schur_comp_, false ) );
462  qp_solver_.warning_tol( warning_tol() );
463  qp_solver_.error_tol( error_tol() );
464  qp_solver_.iter_refine_min_iter( iter_refine_min_iter() );
465  qp_solver_.iter_refine_max_iter( iter_refine_max_iter() );
466  qp_solver_.iter_refine_opt_tol( iter_refine_opt_tol() );
467  qp_solver_.iter_refine_feas_tol( iter_refine_feas_tol() );
468  qp_solver_.iter_refine_at_solution( iter_refine_at_solution() );
469  qp_solver_.pivot_tols(
470  MatrixSymAddDelUpdateable::PivotTolerances(
471  pivot_warning_tol(), pivot_singular_tol(), pivot_wrong_inertia_tol()
472  ));
473 
474  //
475  // Solve the QP with QPSchur
476  //
477  DVector _x(nd+1); // solution vector [ d; eta ]
478  SpVector _mu; // lagrange multipliers for variable bounds [ nu; kappa ]
479  SpVector _lambda_breve; // solution for extra general constraints [ mu; lambda ]
480  size_type qp_iter = 0, num_adds = 0, num_drops = 0;
482  solve_returned
483  = qp_solver_.solve_qp(
484  qp_
485  ,num_act_change, num_act_change ? &ij_act_change[0] : NULL
486  ,num_act_change ? &bnds[0] : NULL
487  ,out, qpschur_olevel
488  ,test_what==RUN_TESTS ? QPSchur::RUN_TESTS : QPSchur::NO_TESTS
489  ,&_x(), &_mu, NULL, &_lambda_breve
490  ,&qp_iter, &num_adds, &num_drops
491  );
492 
493  // Set the solution
494 
495  // d
496  (VectorDenseMutableEncap(*d))() = _x(1,nd);
497  // nu
498  if( nu ) {
499  *nu = 0.0;
500  const SpVector::difference_type o = _mu.offset();
501  if( _mu.nz() ) {
502  for(SpVector::const_iterator _mu_itr = _mu.begin(); _mu_itr != _mu.end(); ++_mu_itr) {
503  typedef SpVector::element_type ele_t;
504  if( _mu_itr->index() + o <= nd ) // don't add multiplier for eta <= etaL
505  nu->set_ele( _mu_itr->index() + o, _mu_itr->value() );
506  }
507  }
508  }
509  // mu, lambda
510  if( m_in || m_eq ) {
511  *eta = _x(nd+1); // must be non-null
512  *mu = 0.0;
513  const SpVector::difference_type o = _lambda_breve.offset();
514  if(_lambda_breve.nz()) {
515  for(SpVector::const_iterator itr = _lambda_breve.begin();
516  itr != _lambda_breve.end();
517  ++itr)
518  {
519  typedef SpVector::element_type ele_t;
520  if( itr->index() + o <= m_in ) {
521  mu->set_ele( itr->index() + o, itr->value() );
522  }
523  else {
524  lambda->set_ele( itr->index() + o - m_in, itr->value() );
525  }
526  }
527  }
528  }
529  // obj_d (This could be updated within QPSchur in the future)
530  if(obj_d) {
531  // obj_d = g'*d + 1/2 * d' * G * g
532  *obj_d = AbstractLinAlgPack::dot(g,*d)
534  }
535  // Ed
536  if(Ed && E) {
537  switch(constraints_->inequality_pick_policy()) {
538  case constr_t::ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY:
539  if(solve_returned == QPSchur::OPTIMAL_SOLUTION)
540  break; // Ed already computed (see ConstraintsRelaxedStd::pick_violated())
541  case constr_t::ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY:
542  break; // Ed already computed (see ConstraintsRelaxedStd::pick_violated())
543  default:
544  // We need to compute Ed
545  LinAlgOpPack::V_MtV( Ed, *E, trans_E, *d );
546  }
547  }
548  // Fd (This could be updated within ConstraintsRelaxedStd in the future)
549  if(Fd) {
550  LinAlgOpPack::V_MtV( Fd, *F, trans_F, *d );
551  }
552  // Set the QP statistics
553  QPSolverStats::ESolutionType solution_type;
554  QPSolverStats::EConvexity convexity = QPSolverStats::CONVEX;
555  switch( solve_returned ) {
556  case QPSchur::OPTIMAL_SOLUTION:
557  solution_type = QPSolverStats::OPTIMAL_SOLUTION;
558  break;
559  case QPSchur::MAX_ITER_EXCEEDED:
560  solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
561  break;
562  case QPSchur::MAX_RUNTIME_EXEEDED_FAIL:
563  solution_type = QPSolverStats::SUBOPTIMAL_POINT;
564  break;
565  case QPSchur::MAX_RUNTIME_EXEEDED_DUAL_FEAS:
566  solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
567  break;
568  case QPSchur::MAX_ALLOWED_STORAGE_EXCEEDED:
569  solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
570  break;
571  case QPSchur::INFEASIBLE_CONSTRAINTS:
572  case QPSchur::NONCONVEX_QP:
573  convexity = QPSolverStats::NONCONVEX;
574  case QPSchur::DUAL_INFEASIBILITY:
575  case QPSchur::SUBOPTIMAL_POINT:
576  solution_type = QPSolverStats::SUBOPTIMAL_POINT;
577  break;
578  default:
580  }
581  qp_stats_.set_stats(
582  solution_type,convexity,qp_iter,num_adds,num_drops
583  , num_act_change > 0 || n_X > 1, *eta > 0.0 );
584 
585  return qp_stats_.solution_type();
586 }
587 
588 } // end namespace ConstrainedOptPack
void set_stats(ESolutionType solution_type, EConvexity convexity, int num_qp_iter, int num_adds, int num_drops, bool warm_start, bool infeasible_qp)
Initialize the statistics.
ERunTests
Enumeration for if to run internal tests or not.
QPSolverStats::ESolutionType imp_solve_qp(std::ostream *out, EOutputLevel olevel, ERunTests test_what, const Vector &g, const MatrixSymOp &G, value_type etaL, const Vector *dL, const Vector *dU, const MatrixOp *E, BLAS_Cpp::Transp trans_E, const Vector *b, const Vector *eL, const Vector *eU, const MatrixOp *F, BLAS_Cpp::Transp trans_F, const Vector *f, value_type *obj_d, value_type *eta, VectorMutable *d, VectorMutable *nu, VectorMutable *mu, VectorMutable *Ed, VectorMutable *lambda, VectorMutable *Fd)
void initialize(const DVectorSlice &g, const MatrixSymOp &G, const MatrixOp *A, size_type n_R, const size_type i_x_free[], const size_type i_x_fixed[], const EBounds bnd_fixed[], const DVectorSlice &b_X, const MatrixSymOpNonsing &Ko, const DVectorSlice &fo, Constraints *constraints, std::ostream *out=NULL, bool test_setup=false, value_type warning_tol=1e-10, value_type error_tol=1e-5, bool print_all_warnings=false)
Initialize.
size_type rows(size_type rows, size_type cols, BLAS_Cpp::Transp _trans)
EConvexity
Enumeration for the type of projected QP on output.
T_To & dyn_cast(T_From &from)
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
void initialize(const G_ptr_t &G_ptr, const vec_mut_ptr_t &M_diag_ptr, const space_ptr_t &space=Teuchos::null)
Initialize the Hessian and the relaxation terms.
void V_mV(VectorMutable *v_lhs, const V &V_rhs)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
size_t size_type
void pivot_tols(MSADU::PivotTolerances pivot_tols)
Set the tolerances to use when updating the schur complement.
Class for storing statistics about a run of a (active set?) QP solver.
value_type dot(const Vector &v_rhs1, const Vector &v_rhs2)
Constraints subclass that is used to represent generic varaible bounds, and general inequality and eq...
virtual ESolveReturn solve_qp(QP &qp, size_type num_act_change, const int ij_act_change[], const EBounds bnds[], std::ostream *out, EOutputLevel output_level, ERunTests test_what, DVectorSlice *x, SpVector *mu, DVectorSlice *lambda, SpVector *lambda_breve, size_type *iter, size_type *num_adds, size_type *num_drops)
Solve a QP.
EOutputLevel
Enumeration for the amount of output to create from solve_qp().
Transp
ESolutionType
Enumeration for the type of point returned from solve_qp(...).
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
QPSolverRelaxedQPSchur(const init_kkt_sys_ptr_t &init_kkt_sys=Teuchos::null, const constraints_ptr_t &constraints=Teuchos::rcp(new QPSchurPack::ConstraintsRelaxedStd), value_type max_qp_iter_frac=10.0, value_type max_real_runtime=1e+20, QPSchurPack::ConstraintsRelaxedStd::EInequalityPickPolicy inequality_pick_policy=QPSchurPack::ConstraintsRelaxedStd::ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY, ELocalOutputLevel print_level=USE_INPUT_ARG, value_type bounds_tol=-1.0, value_type inequality_tol=-1.0, value_type equality_tol=-1.0, value_type loose_feas_tol=-1.0, value_type dual_infeas_tol=-1.0, value_type huge_primal_step=-1.0, value_type huge_dual_step=-1.0, value_type bigM=1e+10, value_type warning_tol=1e-10, value_type error_tol=1e-5, size_type iter_refine_min_iter=1, size_type iter_refine_max_iter=3, value_type iter_refine_opt_tol=1e-12, value_type iter_refine_feas_tol=1e-12, bool iter_refine_at_solution=true, value_type pivot_warning_tol=1e-8, value_type pivot_singular_tol=1e-11, value_type pivot_wrong_inertia_tol=1e-11, bool add_equalities_initially=true)