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"
58 namespace ConstrainedOptPack {
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
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
76 ,value_type warning_tol
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
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)
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)
136 ,
const Vector& g,
const MatrixSymOp& G
138 ,
const Vector* dL,
const Vector* dU
140 ,
const Vector* eL,
const Vector* eU
143 ,value_type* eta, VectorMutable* d
145 ,VectorMutable* mu, VectorMutable* Ed
146 ,VectorMutable* lambda, VectorMutable* Fd
149 namespace mmp = MemMngPack;
154 #ifdef PROFILE_HACK_ENABLED
163 VectorDenseEncap g_de(g);
165 VectorSpace::space_ptr_t
166 space_d_eta = d->space().small_vec_spc_fcty()->create_vec_spc(nd+1);
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
191 const bool all_f_undecomp = F ? j_f_decomp.size() == 0 :
true;
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);
205 constraints_->initialize(
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
210 ,!add_equalities_initially()
219 g_relaxed_.resize(nd+1);
220 g_relaxed_(1,nd) = g_de();
221 g_relaxed_(nd+1) = bigM();
224 bigM_vec_.initialize(1);
227 Teuchos::rcp(&dyn_cast<const MatrixSymOpNonsing>(G),
false)
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)
244 typedef std::vector<int> ij_act_change_t;
245 typedef std::vector<EBounds> bnds_t;
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);
253 if( m_eq && add_equalities_initially() ) {
255 ij_act_change[num_act_change] = (nd + 1) + m_in + j;
256 bnds[num_act_change] = EQUALITY;
266 const value_type inf_bnd = this->infinite_bound();
267 VectorDenseEncap dL_de(*dL);
268 VectorDenseEncap dU_de(*dU);
271 dLU_itr( dL_de().begin(), dL_de().end()
272 ,dU_de().begin(), dU_de().end()
274 for( ; !dLU_itr.
at_end(); ++dLU_itr ) {
276 ij_act_change[num_act_change] = dLU_itr.
index();
277 bnds[num_act_change] = EQUALITY;
283 if( ( nu && nu->nz() ) || ( m_in && mu->nz() ) ) {
288 nu_nz = nu ? nu->nz() : 0,
289 mu_nz = mu ? mu->nz() : 0;
294 SpVector gamma( nd + 1 + m_in , nu_nz + mu_nz );
295 typedef SpVector::element_type ele_t;
297 VectorDenseEncap nu_de(*nu);
298 DVectorSlice::const_iterator
299 nu_itr = nu_de().begin(),
300 nu_end = nu_de().end();
302 while( nu_itr != nu_end ) {
303 for( ; *nu_itr == 0.0; ++nu_itr, ++i );
304 gamma.add_element(ele_t(i,*nu_itr));
309 VectorDenseEncap mu_de(*mu);
310 DVectorSlice::const_iterator
311 mu_itr = mu_de().begin(),
312 mu_end = mu_de().end();
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));
320 std::sort( gamma.begin(), gamma.end()
326 const SpVector::difference_type o = gamma.offset();
327 for( SpVector::const_iterator itr = gamma.begin(); itr != gamma.end(); ++itr ) {
329 if( i <= nd && x_init(i) != FREE )
332 ij_act_change[num_act_change] = i;
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 ) {
353 if( *nu_itr != 0.0 ) {
356 if( ( *x_init_itr ==
LOWER && *nu_itr > 0 )
357 || ( *x_init_itr ==
UPPER && *nu_itr < 0 ) )
360 ij_act_change[num_act_change] = i;
368 ij_act_change[num_act_change] = -i;
369 bnds[num_act_change] = FREE;
377 ij_act_change[num_act_change] = -int(nd+1);
378 bnds[num_act_change] = FREE;
384 switch( print_level() ) {
385 case USE_INPUT_ARG: {
389 qpschur_olevel = QPSchur::NO_OUTPUT;
391 case PRINT_BASIC_INFO:
392 qpschur_olevel = QPSchur::OUTPUT_BASIC_INFO;
394 case PRINT_ITER_SUMMARY:
395 qpschur_olevel = QPSchur::OUTPUT_ITER_SUMMARY;
397 case PRINT_ITER_STEPS:
398 qpschur_olevel = QPSchur::OUTPUT_ITER_STEPS;
400 case PRINT_ITER_ACT_SET:
401 case PRINT_ITER_VECTORS:
402 qpschur_olevel = QPSchur::OUTPUT_ACT_SET;
404 case PRINT_EVERY_THING:
405 qpschur_olevel = QPSchur::OUTPUT_ITER_QUANTITIES;
413 qpschur_olevel = QPSchur::NO_OUTPUT;
415 case OUTPUT_BASIC_INFO:
416 qpschur_olevel = QPSchur::OUTPUT_BASIC_INFO;
418 case OUTPUT_ITER_SUMMARY:
419 qpschur_olevel = QPSchur::OUTPUT_ITER_SUMMARY;
421 case OUTPUT_ITER_STEPS:
422 qpschur_olevel = QPSchur::OUTPUT_ITER_STEPS;
425 qpschur_olevel = QPSchur::OUTPUT_ACT_SET;
427 case OUTPUT_ITER_QUANTITIES:
428 qpschur_olevel = QPSchur::OUTPUT_ITER_QUANTITIES;
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());
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() );
451 if(loose_feas_tol() > 0.0)
452 qp_solver_.loose_feas_tol( loose_feas_tol() );
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() );
470 MatrixSymAddDelUpdateable::PivotTolerances(
471 pivot_warning_tol(), pivot_singular_tol(), pivot_wrong_inertia_tol()
479 SpVector _lambda_breve;
480 size_type qp_iter = 0, num_adds = 0, num_drops = 0;
485 ,num_act_change, num_act_change ? &ij_act_change[0] : NULL
486 ,num_act_change ? &bnds[0] : NULL
488 ,test_what==RUN_TESTS ? QPSchur::RUN_TESTS : QPSchur::NO_TESTS
489 ,&_x(), &_mu, NULL, &_lambda_breve
490 ,&qp_iter, &num_adds, &num_drops
496 (VectorDenseMutableEncap(*d))() = _x(1,nd);
500 const SpVector::difference_type o = _mu.offset();
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 )
505 nu->set_ele( _mu_itr->index() + o, _mu_itr->value() );
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();
519 typedef SpVector::element_type ele_t;
520 if( itr->index() + o <= m_in ) {
521 mu->set_ele( itr->index() + o, itr->value() );
524 lambda->set_ele( itr->index() + o - m_in, itr->value() );
537 switch(constraints_->inequality_pick_policy()) {
538 case constr_t::ADD_BOUNDS_THEN_MOST_VIOLATED_INEQUALITY:
539 if(solve_returned == QPSchur::OPTIMAL_SOLUTION)
541 case constr_t::ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY:
545 LinAlgOpPack::V_MtV( Ed, *E, trans_E, *d );
550 LinAlgOpPack::V_MtV( Fd, *F, trans_F, *d );
555 switch( solve_returned ) {
556 case QPSchur::OPTIMAL_SOLUTION:
557 solution_type = QPSolverStats::OPTIMAL_SOLUTION;
559 case QPSchur::MAX_ITER_EXCEEDED:
560 solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
562 case QPSchur::MAX_RUNTIME_EXEEDED_FAIL:
563 solution_type = QPSolverStats::SUBOPTIMAL_POINT;
565 case QPSchur::MAX_RUNTIME_EXEEDED_DUAL_FEAS:
566 solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
568 case QPSchur::MAX_ALLOWED_STORAGE_EXCEEDED:
569 solution_type = QPSolverStats::DUAL_FEASIBLE_POINT;
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;
582 solution_type,convexity,qp_iter,num_adds,num_drops
583 , num_act_change > 0 || n_X > 1, *eta > 0.0 );
value_type ubound() const
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.
const x_init_t & x_init() const
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)
~QPSolverRelaxedQPSchur()
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)
std::vector< EBounds > bnd_fixed_t
EOutputLevel
Output level.
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.
value_type lbound() const
void V_mV(VectorMutable *v_lhs, const V &V_rhs)
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
QPSolverStats get_qp_stats() const
void pivot_tols(MSADU::PivotTolerances pivot_tols)
Set the tolerances to use when updating the schur complement.
ESolveReturn
solve_qp return values
Class for storing statistics about a run of a (active set?) QP solver.
ELocalOutputLevel
Output level.
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().
std::vector< size_type > i_x_fixed_t
Utility class for a ranged check vector.
std::vector< size_type > j_f_decomp_t
std::vector< size_type > i_x_free_t
ESolutionType
Enumeration for the type of point returned from solve_qp(...).
#define TEUCHOS_TEST_FOR_EXCEPT(throw_exception_test)
ESolutionType solution_type() const
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)