58 namespace ConstrainedOptPack {
61 const init_kkt_sys_ptr_t& init_kkt_sys
62 ,
const constraints_ptr_t& constraints
66 inequality_pick_policy
82 ,
bool iter_refine_at_solution
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
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()
227 Teuchos::rcp(&dyn_cast<const MatrixSymOpNonsing>(G),
false)
234 ,n_R, i_x_free.size() ? &i_x_free[0] : NULL
235 ,&i_x_fixed[0],&bnd_fixed[0]
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;
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();
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
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 ) {
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() ) {
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() );
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() );
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()
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
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:
541 case constr_t::ADD_MOST_VIOLATED_BOUNDS_AND_INEQUALITY:
555 switch( solve_returned ) {
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
Iterate through a set of sparse bounds.
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.
AbstractLinAlgPack::size_type size_type
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)
RTOp_index_type index_type
VectorMutableDense bigM_vec_
SparseVector< SparseElement< index_type, value_type >, std::allocator< SparseElement< index_type, value_type > > > SpVector
QPSchurPack::QPInitFixedFreeStd qp_
~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)
Return rows of a possible transposed matrix.
EConvexity
Enumeration for the type of projected QP on output.
InitKKTSystem::Ko_ptr_t Ko_
int resize(OrdinalType length_in)
Helper class that takes care of timing.
std::vector< EBounds > bnd_fixed_t
MatrixSymHessianRelaxNonSing G_relaxed_
EOutputLevel
Output level.
value_type transVtMtV(const Vector &v_rhs1, const MatrixOp &M_rhs2, BLAS_Cpp::Transp trans_rhs2, const Vector &v_rhs3)
result = v_rhs1' * op(M_rhs2) * 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
MatrixSymAddDelBunchKaufman schur_comp_
void V_mV(VectorMutable *v_lhs, const V &V_rhs)
v_lhs = - V_rhs.
TEUCHOS_DEPRECATED RCP< T > rcp(T *p, Dealloc_T dealloc, bool owns_mem)
T_To & dyn_cast(T_From &from)
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)
result = v_rhs1' * 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().
void V_MtV(VectorMutable *v_lhs, const MatrixOp &M_rhs1, BLAS_Cpp::Transp trans_rhs1, const V &V_rhs2)
v_lhs = op(M_rhs1) * V_rhs2.
AbstractLinAlgPack::value_type value_type
std::vector< size_type > i_x_fixed_t
Utility class for a ranged check vector.
Function object class for sorting a sparse vectors in descending order by abs(v(i)).
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)