1+ #include < iostream>
2+
13#include " eiquadprog/eiquadprog-fast.hpp"
4+ #define TRACE_SOLVER
25
36namespace eiquadprog {
47namespace solvers {
58
6-
79EiquadprogFast::EiquadprogFast () {
810 m_maxIter = DEFAULT_MAX_ITER;
9- q = 0 ; // size of the active set A (containing the indices of the active constraints)
11+ q = 0 ; // size of the active set A (containing the indices
12+ // of the active constraints)
1013 is_inverse_provided_ = false ;
1114 m_nVars = 0 ;
1215 m_nEqCon = 0 ;
@@ -182,7 +185,6 @@ void EiquadprogFast::delete_constraint(MatrixXd& R, MatrixXd& J, VectorXi& A, Ve
182185 }
183186}
184187
185-
186188EiquadprogFast_status EiquadprogFast::solve_quadprog (const MatrixXd& Hess, const VectorXd& g0, const MatrixXd& CE,
187189 const VectorXd& ce0, const MatrixXd& CI, const VectorXd& ci0,
188190 VectorXd& x) {
@@ -232,7 +234,8 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
232234 R.setZero (nVars, nVars);
233235 R_norm = 1.0 ;
234236
235- /* compute the inverse of the factorized matrix Hess^-1, this is the initial value for H */
237+ /* compute the inverse of the factorized matrix Hess^-1,
238+ this is the initial value for H */
236239 // m_J = L^-T
237240 if (!is_inverse_provided_) {
238241 START_PROFILER_EIQUADPROG_FAST (EIQUADPROG_FAST_CHOLESKY_INVERSE);
@@ -253,7 +256,8 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
253256 /* c1 * c2 is an estimate for cond(Hess) */
254257
255258 /*
256- * Find the unconstrained minimizer of the quadratic form 0.5 * x Hess x + g0 x
259+ * Find the unconstrained minimizer of the quadratic
260+ * form 0.5 * x Hess x + g0 x
257261 * this is a feasible point in the dual space
258262 * x = Hess^-1 * g0
259263 */
@@ -296,10 +300,12 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
296300 print_vector (" d" , d, nVars);
297301#endif
298302
299- /* compute full step length t2: i.e., the minimum step in primal space s.t. the contraint
300- becomes feasible */
303+ /* compute full step length t2: i.e.,
304+ the minimum step in primal space s.t. the contraint
305+ becomes feasible */
301306 t2 = 0.0 ;
302- if (std::abs (z.dot (z)) > std::numeric_limits<double >::epsilon ()) // i.e. z != 0
307+ if (std::abs (z.dot (z)) > std::numeric_limits<double >::epsilon ())
308+ // i.e. z != 0
303309 t2 = (-np.dot (x) - ce0 (i)) / z.dot (np);
304310
305311 x += t2 * z;
@@ -326,7 +332,8 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
326332 for (i = 0 ; i < nIneqCon; i++) iai (i) = static_cast <VectorXi::Scalar>(i);
327333
328334#ifdef USE_WARM_START
329- // DEBUG_STREAM("Gonna warm start using previous active set:\n"<<A.transpose()<<"\n")
335+ // DEBUG_STREAM("Gonna warm start using previous active
336+ // set:\n"<<A.transpose()<<"\n")
330337 for (i = nEqCon; i < q; i++) {
331338 iai (i - nEqCon) = -1 ;
332339 ip = A (i);
@@ -335,8 +342,9 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
335342 update_z (z, m_J, d, iq);
336343 update_r (R, r, d, iq);
337344
338- /* compute full step length t2: i.e., the minimum step in primal space s.t. the contraint
339- becomes feasible */
345+ /* compute full step length t2: i.e.,
346+ the minimum step in primal space s.t. the contraint
347+ becomes feasible */
340348 t2 = 0.0 ;
341349 if (std::abs (z.dot (z)) > std::numeric_limits<double >::epsilon ()) // i.e. z != 0
342350 t2 = (-np.dot (x) - ci0 (ip)) / z.dot (np);
@@ -408,7 +416,8 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
408416 if (std::abs (psi) <= static_cast <double >(nIneqCon) * std::numeric_limits<double >::epsilon () * c1 * c2 * 100.0 ) {
409417 /* numerically there are not infeasibilities anymore */
410418 q = iq;
411- // DEBUG_STREAM("Optimal active set:\n"<<A.head(iq).transpose()<<"\n\n")
419+ // DEBUG_STREAM("Optimal active
420+ // set:\n"<<A.head(iq).transpose()<<"\n\n")
412421 return EIQUADPROG_FAST_OPTIMAL;
413422 }
414423
@@ -419,7 +428,8 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
419428
420429l2: /* Step 2: check for feasibility and determine a new S-pair */
421430 START_PROFILER_EIQUADPROG_FAST (EIQUADPROG_FAST_STEP_2);
422- // find constraint with highest violation (what about normalizing constraints?)
431+ // find constraint with highest violation
432+ // (what about normalizing constraints?)
423433 for (i = 0 ; i < nIneqCon; i++) {
424434 if (s (i) < ss && iai (i) != -1 && iaexcl (i)) {
425435 ss = s (i);
@@ -449,7 +459,8 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
449459
450460l2a: /* Step 2a: determine step direction */
451461 START_PROFILER_EIQUADPROG_FAST (EIQUADPROG_FAST_STEP_2A);
452- /* compute z = H np: the step direction in the primal space (through m_J, see the paper) */
462+ /* compute z = H np: the step direction in the primal space
463+ (through m_J, see the paper) */
453464 compute_d (d, m_J, np);
454465 // update_z(z, m_J, d, iq);
455466 if (iq >= nVars) {
@@ -458,7 +469,8 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
458469 } else {
459470 update_z (z, m_J, d, iq);
460471 }
461- /* compute N* np (if q > 0): the negative of the step direction in the dual space */
472+ /* compute N* np (if q > 0): the negative of the
473+ step direction in the dual space */
462474 update_r (R, r, d, iq);
463475#ifdef TRACE_SOLVER
464476 std::cerr << " Step direction z" << std::endl;
@@ -473,7 +485,8 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
473485 /* Step 2b: compute step length */
474486 START_PROFILER_EIQUADPROG_FAST (EIQUADPROG_FAST_STEP_2B);
475487 l = 0 ;
476- /* Compute t1: partial step length (maximum step in dual space without violating dual feasibility */
488+ /* Compute t1: partial step length (maximum step in dual
489+ space without violating dual feasibility */
477490 t1 = inf; /* +inf */
478491 /* find the index l s.t. it reaches the minimum of u+(x) / r */
479492 // l: index of constraint to drop (maybe)
@@ -484,8 +497,10 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
484497 l = A (k);
485498 }
486499 }
487- /* Compute t2: full step length (minimum step in primal space such that the constraint ip becomes feasible */
488- if (std::abs (z.dot (z)) > std::numeric_limits<double >::epsilon ()) // i.e. z != 0
500+ /* Compute t2: full step length (minimum step in primal
501+ space such that the constraint ip becomes feasible */
502+ if (std::abs (z.dot (z)) > std::numeric_limits<double >::epsilon ())
503+ // i.e. z != 0
489504 t2 = -s (ip) / z.dot (np);
490505 else
491506 t2 = inf; /* +inf */
@@ -591,5 +606,3 @@ EiquadprogFast_status EiquadprogFast::solve_quadprog(const MatrixXd& Hess, const
591606
592607} /* namespace solvers */
593608} /* namespace eiquadprog */
594-
595-
0 commit comments