qpmad
Eigen-based C++ QP solver.
implementation.h
Go to the documentation of this file.
1 /**
2  @file
3  @author Alexander Sherikov
4 
5  @copyright 2017 Alexander Sherikov. Licensed under the Apache License,
6  Version 2.0. (see LICENSE or http://www.apache.org/licenses/LICENSE-2.0)
7 
8  @brief Solver implementation
9 */
10 
11 
12 #pragma once
13 
14 #include "common.h"
15 #include "givens.h"
16 #include "inverse.h"
17 #include "solver_parameters.h"
18 #include "constraint_status.h"
19 #include "chosen_constraint.h"
20 #include "active_set.h"
21 #include "factorization_data.h"
22 
23 #ifdef QPMAD_ENABLE_TRACING
24 # include "testing.h"
25 #endif
26 
27 
28 
29 namespace qpmad
30 {
31  template <typename t_Scalar, int t_primal_size, int t_has_bounds, int t_general_ctr_number>
32  class SolverBase
33  {
34  public:
36  {
37  OK = 0,
39  UNDEFINED = -1 // used for initialization of status variables
40  };
41 
42  template <int t_rows>
43  using Vector = Eigen::Matrix<t_Scalar, t_rows, 1>;
44  template <int t_rows, int t_cols>
45  using Matrix = Eigen::Matrix<t_Scalar, t_rows, t_cols>; /// @deprecated
46  using Scalar = t_Scalar;
47 
48 
49  protected:
50  using MatrixRef = Eigen::Ref<Eigen::Matrix<t_Scalar, Eigen::Dynamic, Eigen::Dynamic>>;
51  using VectorRef = Eigen::Ref<Eigen::Matrix<t_Scalar, Eigen::Dynamic, 1>>;
52 
53 
54  protected:
56 
59 
62 
68 
70  Eigen::Dynamic == t_general_ctr_number ?
71  Eigen::Dynamic :
72  (0 == t_has_bounds ? t_general_ctr_number :
73  (Eigen::Dynamic == t_primal_size ? Eigen::Dynamic :
74  t_general_ctr_number + t_primal_size));
75 
76 
78 
79  std::ptrdiff_t iter_counter_;
81 
84 
85  Eigen::Array<uint8_t, num_constraints_compile_time_, 1> constraints_status_;
86 
87 
88 
89  public:
91  {
92  primal_size_ = 0;
93  h_size_ = 0;
96 
97  iter_counter_ = 0;
98  machinery_initialized_ = false;
100  }
101 
102 
103  /**
104  * @brief Returns type of the Hessian produced by the latest execution
105  * of `solve()`.
106  */
108  {
109  return (hessian_type_);
110  }
111 
112 
113  /**
114  * @brief Returns number of inequality iterations during the latest
115  * execution of `solve()`.
116  */
117  std::ptrdiff_t getNumberOfInequalityIterations() const
118  {
119  return (iter_counter_);
120  }
121 
122 
123  /**
124  * @brief Returns dual variables (Lagrange multipliers) corresponding
125  * to inequality constraints. Must be called after successful
126  * `solve()`, the result is undefined if previous call to `solve()`
127  * failed.
128  *
129  * @tparam t_status_size
130  * @tparam t_dual_size
131  * @tparam t_index_size
132  *
133  * @param[out] dual dual variables
134  * @param[out] indices constraint indices corresponding to the dual
135  * variables, index 0 corresponds to the first simple bound if present
136  * or to the first general constraint otherwise
137  * @param[out] is_lower flags indicating if lower or upper bound is
138  * active
139  */
140  template <int t_status_size, int t_dual_size, int t_index_size>
142  Vector<t_dual_size> &dual,
143  Eigen::Matrix<MatrixIndex, t_index_size, 1> &indices,
144  Eigen::Matrix<bool, t_status_size, 1> &is_lower) const
145  {
147  {
148  const MatrixIndex size = active_set_.size_ - active_set_.num_equalities_;
149 
150  dual.resize(size);
151  indices.resize(size);
152  is_lower.resize(size);
153 
154  for (MatrixIndex i = active_set_.num_equalities_; i < active_set_.size_; ++i)
155  {
156  const std::size_t output_index = i - active_set_.num_equalities_;
157 
158  dual(output_index) = dual_(i);
159  indices(output_index) = active_set_.getIndex(i);
160  is_lower(output_index) =
162  }
163  }
164  else
165  {
166  dual.resize(0);
167  indices.resize(0);
168  is_lower.resize(0);
169  }
170  }
171 
172 
173  void reserve(
174  const MatrixIndex primal_size,
175  const MatrixIndex num_simple_bounds,
176  const MatrixIndex num_general_constraints)
177  {
178  reserveMachinery(primal_size, num_general_constraints);
179  constraints_status_.resize(num_simple_bounds + num_general_constraints);
180  reserveDual(primal_size);
181  }
182 
183 
184  protected:
185  template <class t_primal, class t_H, class t_h, class t_lb, class t_ub, class t_A, class t_Alb, class t_Aub>
187  t_primal &primal,
188  t_H &H,
189  const t_h &h,
190  const t_lb &lb,
191  const t_ub &ub,
192  const t_A &A,
193  const t_Alb &Alb,
194  const t_Aub &Aub,
195  const SolverParameters &param)
196  {
197  QPMAD_TRACE(std::setprecision(std::numeric_limits<double>::digits10));
198 
199  machinery_initialized_ = false;
200  iter_counter_ = 0;
201 
202 
204 
205  if (0 == h_size_)
206  {
207  // trivial unconstrained optimum
208  setZero(primal);
209 
210  if (0 == num_constraints_)
211  {
212  // trivial solution
213  return (OK);
214  }
215  }
216 
217 
218  switch (param.hessian_type_)
219  {
222  /* Falls through. */
225  // unconstrained optimum
226  if (h_size_ > 0)
227  {
228  primal = -h;
229  H.template triangularView<Eigen::Lower>().solveInPlace(primal);
230  H.transpose().template triangularView<Eigen::Upper>().solveInPlace(primal);
231  }
232  break;
233 
236  // unconstrained optimum
237  if (h_size_ > 0)
238  {
239  primal_step_direction_.noalias() = H.template triangularView<Eigen::Upper>().transpose() * -h;
240  primal.noalias() = H.template triangularView<Eigen::Upper>() * primal_step_direction_;
241  }
242  break;
243 
244  default:
245  QPMAD_UTILS_THROW("Malformed solver parameters!");
246  break;
247  }
248 
249 
250  if (0 == num_constraints_)
251  {
252  // return unconstrained optimum
253  return (OK);
254  }
255 
256 
257  // check consistency of general constraints and activate
258  // equality constraints
260  MatrixIndex num_equalities = 0;
261  for (MatrixIndex i = 0; i < num_constraints_; ++i)
262  {
264 
265  double lb_i;
266  double ub_i;
267 
268  if (true == chosen_ctr_.is_simple_)
269  {
270  lb_i = lb(i);
271  ub_i = ub(i);
272  }
273  else
274  {
278  }
279 
280 
281  if (lb_i - param.tolerance_ > ub_i)
282  {
284  QPMAD_UTILS_THROW("Inconsistent constraints (lb > ub).");
285  }
286 
287  if (std::abs(lb_i - ub_i) > param.tolerance_)
288  {
290  }
291  else
292  {
294  ++num_equalities;
295 
296 
297  if (true == chosen_ctr_.is_simple_)
298  {
299  chosen_ctr_.violation_ = lb_i - primal(i);
300  }
301  else
302  {
304  }
305 
307 
308 
309  // if 'primal_size_' constraints are already activated
310  // all other constraints are linearly dependent
311  if (active_set_.hasEmptySpace())
312  {
313  double ctr_i_dot_primal_step_direction;
314 
315  if (true == chosen_ctr_.is_simple_)
316  {
317  factorization_data_.computeEqualityPrimalStep(primal_step_direction_, i, active_set_.size_);
318 
319  ctr_i_dot_primal_step_direction = primal_step_direction_(i);
320  }
321  else
322  {
323  factorization_data_.computeEqualityPrimalStep(
326  active_set_.size_);
327 
328  ctr_i_dot_primal_step_direction =
330  }
331 
332  // if step direction is a zero vector, constraint is
333  // linearly dependent with previously added constraints
334  if (ctr_i_dot_primal_step_direction < -param.tolerance_)
335  {
336  double primal_step_length_ = chosen_ctr_.violation_ / ctr_i_dot_primal_step_direction;
337 
338  primal.noalias() += primal_step_length_ * primal_step_direction_;
339 
340  if (false
341  == factorization_data_.update(
343  {
344  QPMAD_UTILS_THROW("Failed to add an equality constraint -- is this possible?");
345  }
346  active_set_.addEquality(i);
347 
348  continue;
349  }
350  } // otherwise -- linear dependence
351 
352  // this point is reached if constraint is linearly dependent
353 
354  // check if this constraint is actually satisfied
355  if (std::abs(chosen_ctr_.violation_) > param.tolerance_)
356  {
357  // nope it is not
359  QPMAD_UTILS_THROW("Infeasible equality constraints");
360  }
361  // otherwise keep going
362  }
363  }
364 
365 
366  if (num_equalities == num_constraints_)
367  {
368  // exit early -- avoid unnecessary memory allocations
369  return (OK);
370  }
371 
372 
373  ReturnStatus return_status;
374  chooseConstraint(primal, lb, ub, A, Alb, Aub, param.tolerance_);
375 
376 
377  if (std::abs(chosen_ctr_.violation_) < param.tolerance_)
378  {
379  // all constraints are satisfied
380  return_status = OK;
381  }
382  else
383  {
384  return_status = MAXIMAL_NUMBER_OF_ITERATIONS;
385 
386 
389 
390 
391  double chosen_ctr_dot_primal_step_direction = 0.0;
392 
393  //
394  factorization_data_.computeInequalityDualStep(dual_step_direction_, chosen_ctr_, A, active_set_);
395  if (active_set_.hasEmptySpace())
396  {
397  // compute step direction in primal space
398  factorization_data_.computeInequalityPrimalStep(primal_step_direction_, active_set_);
399  chosen_ctr_dot_primal_step_direction =
401  }
402 
403 
404  // last iteration is not counted, so iter_counter_ starts with 1.
405  for (iter_counter_ = 1; (param.max_iter_ < 0) || (iter_counter_ <= param.max_iter_); ++iter_counter_)
406  {
407  QPMAD_TRACE(">>>>>>>>>" << iter_counter_ << "<<<<<<<<<");
408 #ifdef QPMAD_ENABLE_TRACING
409  testing::computeObjective(H, h, primal);
410 #endif
411  QPMAD_TRACE("||| Chosen ctr index = " << chosen_ctr_.index_);
412  QPMAD_TRACE("||| Chosen ctr dual = " << chosen_ctr_.dual_);
413  QPMAD_TRACE("||| Chosen ctr violation = " << chosen_ctr_.violation_);
414 
415 
416  // check dual feasibility
417  MatrixIndex dual_blocking_index = primal_size_;
418  double dual_step_length = std::numeric_limits<double>::infinity();
419  for (MatrixIndex i = active_set_.num_equalities_; i < active_set_.size_; ++i)
420  {
421  if (dual_step_direction_(i) < -param.tolerance_)
422  {
423  const double dual_step_length_i = -dual_(i) / dual_step_direction_(i);
424  if (dual_step_length_i < dual_step_length)
425  {
426  dual_step_length = dual_step_length_i;
427  dual_blocking_index = i;
428  }
429  }
430  }
431 
432 
433 #ifdef QPMAD_ENABLE_TRACING
435  H,
436  h,
437  primal,
438  A,
439  active_set_,
442  dual_,
444 #endif
445 
446 
447  if (active_set_.hasEmptySpace()
448  // if step direction is a zero vector, constraint is
449  // linearly dependent with previously added constraints
450  && (std::abs(chosen_ctr_dot_primal_step_direction) > param.tolerance_))
451  {
452  double step_length = -chosen_ctr_.violation_ / chosen_ctr_dot_primal_step_direction;
453 
454  QPMAD_TRACE("======================");
455  QPMAD_TRACE("||| Primal step length = " << step_length);
456  QPMAD_TRACE("||| Dual step length = " << dual_step_length);
457  QPMAD_TRACE("======================");
458 
459 
460  bool partial_step = false;
462  (step_length >= 0.0) && (dual_step_length >= 0.0),
463  "Non-negative step lengths expected.");
464  if (dual_step_length <= step_length)
465  {
466  step_length = dual_step_length;
467  partial_step = true;
468  }
469 
470 
471  primal.noalias() += step_length * primal_step_direction_;
472 
473  dual_.segment(active_set_.num_equalities_, active_set_.num_inequalities_).noalias() +=
474  step_length
475  * dual_step_direction_.segment(
476  active_set_.num_equalities_, active_set_.num_inequalities_);
477  chosen_ctr_.dual_ += step_length;
478  chosen_ctr_.violation_ += step_length * chosen_ctr_dot_primal_step_direction;
479 
480  QPMAD_TRACE("||| Chosen ctr dual = " << chosen_ctr_.dual_);
481  QPMAD_TRACE("||| Chosen ctr violation = " << chosen_ctr_.violation_);
482 
483 
484  if ((partial_step)
485  // if violation is almost zero -- assume that a full step is made
486  && (std::abs(chosen_ctr_.violation_) > param.tolerance_))
487  {
488  QPMAD_TRACE("||| PARTIAL STEP");
489  // deactivate blocking constraint
490  constraints_status_[active_set_.getIndex(dual_blocking_index)] = ConstraintStatus::INACTIVE;
491 
492  dropElementWithoutResize(dual_, dual_blocking_index, active_set_.size_);
493 
494  factorization_data_.downdate(dual_blocking_index, active_set_.size_);
495 
496  active_set_.removeInequality(dual_blocking_index);
497 
498  // compute step direction in primal & dual space
499  factorization_data_.updateStepsAfterPartialStep(
501  chosen_ctr_dot_primal_step_direction =
503  }
504  else
505  {
506  QPMAD_TRACE("||| FULL STEP");
507  // activate constraint
508  if (false
509  == factorization_data_.update(
511  {
512  QPMAD_UTILS_THROW("Failed to add an inequality constraint -- is this possible?");
513  }
514 
516  {
518  }
519  else
520  {
522  }
524  active_set_.addInequality(chosen_ctr_.index_);
525 
526  chooseConstraint(primal, lb, ub, A, Alb, Aub, param.tolerance_);
527 
528  if (std::abs(chosen_ctr_.violation_) < param.tolerance_)
529  {
530  // all constraints are satisfied
531  return_status = OK;
532  break;
533  }
534 
535  chosen_ctr_dot_primal_step_direction = 0.0;
536  factorization_data_.computeInequalityDualStep(
538  if (active_set_.hasEmptySpace())
539  {
540  // compute step direction in primal & dual space
541  factorization_data_.computeInequalityPrimalStep(primal_step_direction_, active_set_);
542  chosen_ctr_dot_primal_step_direction =
544  }
545  }
546  }
547  else
548  {
549  if (dual_blocking_index == primal_size_)
550  {
551  QPMAD_UTILS_THROW("Infeasible inequality constraints.");
552  }
553  else
554  {
555  QPMAD_TRACE("======================");
556  QPMAD_TRACE("||| Dual step length = " << dual_step_length);
557  QPMAD_TRACE("======================");
558 
559  // otherwise -- deactivate
560  dual_.segment(active_set_.num_equalities_, active_set_.num_inequalities_).noalias() +=
561  dual_step_length
562  * dual_step_direction_.segment(
563  active_set_.num_equalities_, active_set_.num_inequalities_);
564  chosen_ctr_.dual_ += dual_step_length;
565 
566  constraints_status_[active_set_.getIndex(dual_blocking_index)] = ConstraintStatus::INACTIVE;
567 
568  dropElementWithoutResize(dual_, dual_blocking_index, active_set_.size_);
569 
570  factorization_data_.downdate(dual_blocking_index, active_set_.size_);
571 
572  active_set_.removeInequality(dual_blocking_index);
573 
574  // compute step direction in primal & dual space
575  factorization_data_.updateStepsAfterPureDualStep(
577  chosen_ctr_dot_primal_step_direction =
579  }
580  }
581  }
582  }
583 
584 #ifdef QPMAD_ENABLE_TRACING
586  {
588 
591  }
592  else
593  {
594  QPMAD_TRACE("||| NO ACTIVE CONSTRAINTS");
595  }
596 #endif
597  return (return_status);
598  }
599 
600 
601  private:
602  void reserveMachinery(const MatrixIndex primal_size, const MatrixIndex num_general_constraints)
603  {
604  active_set_.initialize(primal_size);
605  primal_step_direction_.resize(primal_size);
606  general_ctr_dot_primal_.resize(num_general_constraints);
607 
608  factorization_data_.reserve(primal_size);
609  }
610 
611 
612  void reserveDual(const MatrixIndex primal_size)
613  {
614  dual_.resize(primal_size);
615  dual_step_direction_.resize(primal_size);
616  }
617 
618 
619  template <class t_MatrixType>
620  void initializeMachineryLazy(t_MatrixType &H, const bool return_inverted_cholesky_factor)
621  {
623  {
625 
626  factorization_data_.initialize(H, hessian_type_, primal_size_, return_inverted_cholesky_factor);
627  if (return_inverted_cholesky_factor)
628  {
630  }
631 
632  machinery_initialized_ = true;
633  }
634  }
635 
636 
637  template <
638  class t_Primal,
639  class t_LowerBounds,
640  class t_UpperBounds,
641  class t_Constraints,
642  class t_ConstraintsLowerBounds,
643  class t_ConstraintsUpperBounds>
645  const t_Primal &primal,
646  const t_LowerBounds &lb,
647  const t_UpperBounds &ub,
648  const t_Constraints &A,
649  const t_ConstraintsLowerBounds &Alb,
650  const t_ConstraintsUpperBounds &Aub,
651  const double tolerance)
652  {
653  chosen_ctr_.reset();
654 
655 
656  for (MatrixIndex i = 0; i < num_simple_bounds_; ++i)
657  {
659  {
660  checkConstraintViolation(i, lb(i), ub(i), primal(i));
661  }
662  }
663 
664  if ((std::abs(chosen_ctr_.violation_) < tolerance) && (num_general_constraints_ > 0))
665  {
666  general_ctr_dot_primal_.noalias() = A * primal;
668  {
670  {
672  i,
673  Alb(i - num_simple_bounds_),
674  Aub(i - num_simple_bounds_),
676  }
677  }
679  {
681  }
682  }
683 
686  }
687 
688 
690  const MatrixIndex i,
691  const double lb_i,
692  const double ub_i,
693  const double ctr_i_dot_primal)
694  {
695  double ctr_violation_i = ctr_i_dot_primal - lb_i;
696  if (ctr_violation_i < -std::abs(chosen_ctr_.violation_))
697  {
698  chosen_ctr_.violation_ = ctr_violation_i;
699  chosen_ctr_.index_ = i;
700  }
701  else
702  {
703  ctr_violation_i = ctr_i_dot_primal - ub_i;
704  if (ctr_violation_i > std::abs(chosen_ctr_.violation_))
705  {
706  chosen_ctr_.violation_ = ctr_violation_i;
707  chosen_ctr_.index_ = i;
708  }
709  }
710  }
711 
712 
713  template <class t_VectorType, class t_MatrixType>
714  double getConstraintDotPrimalStepDirection(const t_VectorType &primal_step_direction, const t_MatrixType &A)
715  const
716  {
718  {
719  return (primal_step_direction(chosen_ctr_.index_));
720  }
721  else
722  {
723  return (A.row(chosen_ctr_.general_constraint_index_) * primal_step_direction);
724  }
725  }
726 
727 
728  template <int... t_Other>
729  static void factorizeCholeskyInPlace(Eigen::Matrix<t_Scalar, t_Other...> &H)
730  {
731  const Eigen::LLT<Eigen::Ref<typename std::decay<decltype(H)>::type>, Eigen::Lower> llt(H);
733  Eigen::Success == llt.info(), "Could not perform Cholesky decomposition of the Hessian (dense).");
734  }
735 
737  {
738  const Eigen::LLT<Eigen::Ref<typename std::decay<decltype(H)>::type>, Eigen::Lower> llt(H);
740  Eigen::Success == llt.info(), "Could not perform Cholesky decomposition of the Hessian (dense).");
741  }
742 
743 #if !defined(EIGEN_MPL2_ONLY) || (EIGEN_WORLD_VERSION > 3) || (EIGEN_MAJOR_VERSION > 3)
744  template <int t_Options, typename t_StorageIndex>
745  static void factorizeCholeskyInPlace(Eigen::SparseMatrix<t_Scalar, t_Options, t_StorageIndex> &H)
746  {
747  const Eigen::SimplicialLLT<typename std::decay<decltype(H)>::type, Eigen::Lower> llt(H);
749  Eigen::Success == llt.info(), "Could not perform Cholesky decomposition of the Hessian (sparse).");
750  H = llt.matrixL();
751  }
752 #endif
753 
754 
755  template <class t_primal>
756  void setZero(t_primal &primal)
757  {
758  primal.setZero(primal_size_);
759  }
760 
761  void setZero(VectorRef primal)
762  {
763  primal.setZero();
764  }
765  };
766 } // namespace qpmad
MatrixIndex general_constraint_index_
void initializeMachineryLazy(t_MatrixType &H, const bool return_inverted_cholesky_factor)
Vector< t_primal_size > dual_step_direction_
void setZero(t_primal &primal)
void checkConstraintViolation(const MatrixIndex i, const double lb_i, const double ub_i, const double ctr_i_dot_primal)
SolverParameters::HessianType hessian_type_
ReturnStatus solveGeneric(t_primal &primal, t_H &H, const t_h &h, const t_lb &lb, const t_ub &ub, const t_A &A, const t_Alb &Alb, const t_Aub &Aub, const SolverParameters &param)
void setZero(VectorRef primal)
std::ptrdiff_t iter_counter_
MatrixIndex num_constraints_
static void factorizeCholeskyInPlace(Eigen::Matrix< t_Scalar, t_Other... > &H)
FactorizationData< t_Scalar, t_primal_size > factorization_data_
MatrixIndex num_general_constraints_
double getConstraintDotPrimalStepDirection(const t_VectorType &primal_step_direction, const t_MatrixType &A) const
std::ptrdiff_t getNumberOfInequalityIterations() const
Returns number of inequality iterations during the latest execution of solve().
MatrixIndex num_simple_bounds_
Eigen::Matrix< t_Scalar, t_rows, 1 > Vector
SolverParameters::HessianType getHessianType() const
Returns type of the Hessian produced by the latest execution of solve().
Vector< t_primal_size > primal_step_direction_
ChosenConstraint chosen_ctr_
void reserveDual(const MatrixIndex primal_size)
static void factorizeCholeskyInPlace(Eigen::SparseMatrix< t_Scalar, t_Options, t_StorageIndex > &H)
MatrixIndex primal_size_
void getInequalityDual(Vector< t_dual_size > &dual, Eigen::Matrix< MatrixIndex, t_index_size, 1 > &indices, Eigen::Matrix< bool, t_status_size, 1 > &is_lower) const
Returns dual variables (Lagrange multipliers) corresponding to inequality constraints....
Eigen::Matrix< t_Scalar, t_rows, t_cols > Matrix
Vector< t_general_ctr_number > general_ctr_dot_primal_
MatrixIndex h_size_
void chooseConstraint(const t_Primal &primal, const t_LowerBounds &lb, const t_UpperBounds &ub, const t_Constraints &A, const t_ConstraintsLowerBounds &Alb, const t_ConstraintsUpperBounds &Aub, const double tolerance)
Eigen::Ref< Eigen::Matrix< t_Scalar, Eigen::Dynamic, Eigen::Dynamic > > MatrixRef
Vector< t_primal_size > dual_
ActiveSet< t_primal_size > active_set_
void reserveMachinery(const MatrixIndex primal_size, const MatrixIndex num_general_constraints)
void reserve(const MatrixIndex primal_size, const MatrixIndex num_simple_bounds, const MatrixIndex num_general_constraints)
Eigen::Ref< Eigen::Matrix< t_Scalar, Eigen::Dynamic, 1 > > VectorRef
Eigen::Array< uint8_t, num_constraints_compile_time_, 1 > constraints_status_
static const MatrixIndex num_constraints_compile_time_
static void factorizeCholeskyInPlace(MatrixRef &H)
#define QPMAD_TRACE(info)
Definition: common.h:31
#define QPMAD_UTILS_ASSERT(condition,...)
Definition: exception.h:31
#define QPMAD_UTILS_THROW(s)
Definition: exception.h:18
#define QPMAD_UTILS_PERSISTENT_ASSERT(condition,...)
Definition: exception.h:22
double computeObjective(const Eigen::MatrixXd &H, const Eigen::VectorXd &h, const Eigen::VectorXd &primal)
Definition: testing.h:21
void printActiveSet(const t_ActiveSet &active_set, const t_ConstraintStatuses &constraints_status, const Eigen::VectorXd &dual)
Definition: testing.h:165
void checkLagrangeMultipliers(const Eigen::MatrixXd &H, const Eigen::VectorXd &h, const Eigen::VectorXd &primal, const Eigen::MatrixXd &A, const t_ActiveSet &active_set, const MatrixIndex &num_simple_bounds, const t_ConstraintStatuses &constraints_status, const Eigen::VectorXd &dual, const Eigen::VectorXd &dual_direction=Eigen::VectorXd())
Definition: testing.h:39
void dropElementWithoutResize(t_VectorType &vector, const MatrixIndex index, const MatrixIndex size)
Definition: common.h:44
EIGEN_DEFAULT_DENSE_INDEX_TYPE MatrixIndex
Definition: common.h:37