proxsuite 0.6.7
The Advanced Proximal Optimization Toolbox
Loading...
Searching...
No Matches
solver.hpp
Go to the documentation of this file.
1//
2// Copyright (c) 2022 INRIA
3//
6#ifndef PROXSUITE_PROXQP_SPARSE_SOLVER_HPP
7#define PROXSUITE_PROXQP_SPARSE_SOLVER_HPP
8
9#include <chrono>
10#include <cmath>
11
12#include <proxsuite/linalg/dense/core.hpp>
13#include <proxsuite/linalg/sparse/core.hpp>
14#include <proxsuite/linalg/sparse/factorize.hpp>
15#include <proxsuite/linalg/sparse/update.hpp>
16#include <proxsuite/linalg/sparse/rowmod.hpp>
17#include <proxsuite/proxqp/dense/views.hpp>
18#include <proxsuite/proxqp/settings.hpp>
19#include <proxsuite/linalg/veg/vec.hpp>
20#include "proxsuite/proxqp/results.hpp"
21#include "proxsuite/proxqp/sparse/fwd.hpp"
22#include "proxsuite/proxqp/sparse/views.hpp"
23#include "proxsuite/proxqp/sparse/model.hpp"
24#include "proxsuite/proxqp/sparse/workspace.hpp"
25#include "proxsuite/proxqp/sparse/utils.hpp"
26#include "proxsuite/proxqp/sparse/preconditioner/ruiz.hpp"
27#include "proxsuite/proxqp/sparse/preconditioner/identity.hpp"
28
29#include <iostream>
30#include <iomanip>
31#include <Eigen/IterativeLinearSolvers>
32#include <unsupported/Eigen/IterativeSolvers>
33
34namespace proxsuite {
35namespace proxqp {
36namespace sparse {
37
38template<typename T, typename I>
39void
40ldl_solve(VectorViewMut<T> sol,
41 VectorView<T> rhs,
42 isize n_tot,
43 proxsuite::linalg::sparse::MatMut<T, I> ldl,
44 Eigen::MINRES<detail::AugmentedKkt<T, I>,
45 Eigen::Upper | Eigen::Lower,
46 Eigen::IdentityPreconditioner>& iterative_solver,
47 bool do_ldlt,
48 proxsuite::linalg::veg::dynstack::DynStackMut stack,
49 T* ldl_values,
50 I* perm,
51 I* ldl_col_ptrs,
52 I const* perm_inv)
53{
54 LDLT_TEMP_VEC_UNINIT(T, work_, n_tot, stack);
55 auto rhs_e = rhs.to_eigen();
56 auto sol_e = sol.to_eigen();
57 auto zx = proxsuite::linalg::sparse::util::zero_extend;
58
59 if (do_ldlt) {
60
61 for (isize i = 0; i < n_tot; ++i) {
62 work_[i] = rhs_e[isize(zx(perm[i]))];
63 }
64
65 proxsuite::linalg::sparse::dense_lsolve<T, I>( //
66 { proxsuite::linalg::sparse::from_eigen, work_ },
67 ldl.as_const());
68
69 for (isize i = 0; i < n_tot; ++i) {
70 work_[i] /= ldl_values[isize(zx(ldl_col_ptrs[i]))];
71 }
72
73 proxsuite::linalg::sparse::dense_ltsolve<T, I>( //
74 { proxsuite::linalg::sparse::from_eigen, work_ },
75 ldl.as_const());
76
77 for (isize i = 0; i < n_tot; ++i) {
78 sol_e[i] = work_[isize(zx(perm_inv[i]))];
79 }
80 } else {
81 work_ = iterative_solver.solve(rhs_e);
82 sol_e = work_;
83 }
84}
85
86template<typename T, typename I>
87void
88ldl_iter_solve_noalias(
89 VectorViewMut<T> sol,
90 VectorView<T> rhs,
91 VectorView<T> init_guess,
92 Results<T> const& results,
93 Model<T, I> const& data,
94 isize n_tot,
95 proxsuite::linalg::sparse::MatMut<T, I> ldl,
96 Eigen::MINRES<detail::AugmentedKkt<T, I>,
97 Eigen::Upper | Eigen::Lower,
98 Eigen::IdentityPreconditioner>& iterative_solver,
99 bool do_ldlt,
100 proxsuite::linalg::veg::dynstack::DynStackMut stack,
101 T* ldl_values,
102 I* perm,
103 I* ldl_col_ptrs,
104 I const* perm_inv,
105 Settings<T> const& settings,
106 proxsuite::linalg::sparse::MatMut<T, I> kkt_active,
107 proxsuite::linalg::veg::SliceMut<bool> active_constraints)
108{
109 auto rhs_e = rhs.to_eigen();
110 auto sol_e = sol.to_eigen();
111
112 if (init_guess.dim == sol.dim) {
113 sol_e = init_guess.to_eigen();
114 } else {
115 sol_e.setZero();
116 }
117
118 LDLT_TEMP_VEC_UNINIT(T, err, n_tot, stack);
119
120 T prev_err_norm = std::numeric_limits<T>::infinity();
121
122 for (isize solve_iter = 0; solve_iter < settings.nb_iterative_refinement;
123 ++solve_iter) {
124
125 auto err_x = err.head(data.dim);
126 auto err_y = err.segment(data.dim, data.n_eq);
127 auto err_z = err.tail(data.n_in);
128
129 auto sol_x = sol_e.head(data.dim);
130 auto sol_y = sol_e.segment(data.dim, data.n_eq);
131 auto sol_z = sol_e.tail(data.n_in); // removed active set condition
132
133 err = -rhs_e;
134
135 if (solve_iter > 0) {
136 T mu_eq_neg = -results.info.mu_eq;
137 T mu_in_neg = -results.info.mu_in;
138 // switch (settings.merit_function_type) {
139 // case MeritFunctionType::GPDAL:
140 // mu_in_neg = -settings.alpha_gpdal*results.info.mu_in;
141 // break;
142 // case MeritFunctionType::PDAL:
143 // mu_in_neg = -results.info.mu_in;
144 // break;
145 // }
146 detail::noalias_symhiv_add(err, kkt_active.to_eigen(), sol_e);
147 err_x += results.info.rho * sol_x;
148 err_y += mu_eq_neg * sol_y;
149 for (isize i = 0; i < data.n_in; ++i) {
150 err_z[i] += (active_constraints[i] ? mu_in_neg : T(1)) * sol_z[i];
151 }
152 }
153
154 T err_norm = infty_norm(err);
155 if (err_norm > prev_err_norm / T(2)) {
156 break;
157 }
158 prev_err_norm = err_norm;
159
160 ldl_solve({ proxqp::from_eigen, err },
161 { proxqp::from_eigen, err },
162 n_tot,
163 ldl,
164 iterative_solver,
165 do_ldlt,
166 stack,
167 ldl_values,
168 perm,
169 ldl_col_ptrs,
170 perm_inv);
171
172 sol_e -= err;
173 }
174}
199template<typename T, typename I>
200void
201ldl_solve_in_place(
202 VectorViewMut<T> rhs,
203 VectorView<T> init_guess,
204 Results<T> const& results,
205 Model<T, I> const& data,
206 isize n_tot,
207 proxsuite::linalg::sparse::MatMut<T, I> ldl,
208 Eigen::MINRES<detail::AugmentedKkt<T, I>,
209 Eigen::Upper | Eigen::Lower,
210 Eigen::IdentityPreconditioner>& iterative_solver,
211 bool do_ldlt,
212 proxsuite::linalg::veg::dynstack::DynStackMut stack,
213 T* ldl_values,
214 I* perm,
215 I* ldl_col_ptrs,
216 I const* perm_inv,
217 Settings<T> const& settings,
218 proxsuite::linalg::sparse::MatMut<T, I> kkt_active,
219 proxsuite::linalg::veg::SliceMut<bool> active_constraints)
220{
221 LDLT_TEMP_VEC_UNINIT(T, tmp, n_tot, stack);
222 ldl_iter_solve_noalias({ proxqp::from_eigen, tmp },
223 rhs.as_const(),
224 init_guess,
225 results,
226 data,
227 n_tot,
228 ldl,
229 iterative_solver,
230 do_ldlt,
231 stack,
232 ldl_values,
233 perm,
234 ldl_col_ptrs,
235 perm_inv,
236 settings,
237 kkt_active,
238 active_constraints);
239 rhs.to_eigen() = tmp;
240}
248template<typename T, typename I>
249auto
250inner_reconstructed_matrix(proxsuite::linalg::sparse::MatMut<T, I> ldl)
251 -> DMat<T>
252{
253 auto ldl_dense = ldl.to_eigen().toDense();
254 auto l = DMat<T>(ldl_dense.template triangularView<Eigen::UnitLower>());
255 auto lt = l.transpose();
256 auto d = ldl_dense.diagonal().asDiagonal();
257 auto mat = DMat<T>(l * d * lt);
258 return mat;
259}
269template<typename T, typename I>
270auto
271reconstructed_matrix(proxsuite::linalg::sparse::MatMut<T, I> ldl,
272 I const* perm_inv,
273 isize n_tot) -> DMat<T>
274{
275 auto mat = inner_reconstructed_matrix(ldl);
276 auto mat_backup = mat;
277 for (isize i = 0; i < n_tot; ++i) {
278 for (isize j = 0; j < n_tot; ++j) {
279 mat(i, j) = mat_backup(perm_inv[i], perm_inv[j]);
280 }
281 }
282 return mat;
283}
297template<typename T, typename I>
298auto
299reconstruction_error(proxsuite::linalg::sparse::MatMut<T, I> ldl,
300 I const* perm_inv,
301 Results<T> const& results,
302 Model<T, I> const& data,
303 isize n_tot,
304 proxsuite::linalg::sparse::MatMut<T, I> kkt_active,
305 proxsuite::linalg::veg::SliceMut<bool> active_constraints)
306 -> DMat<T>
307{
308 T mu_eq_neg = -results.info.mu_eq;
309 T mu_in_neg = -results.info.mu_in;
310 auto diff = DMat<T>(
311 reconstructed_matrix(ldl, perm_inv, n_tot) -
312 DMat<T>(
313 DMat<T>(kkt_active.to_eigen()).template selfadjointView<Eigen::Upper>()));
314 diff.diagonal().head(data.dim).array() -= results.info.rho;
315 diff.diagonal().segment(data.dim, data.n_eq).array() -= mu_eq_neg;
316 for (isize i = 0; i < data.n_in; ++i) {
317 diff.diagonal()[data.dim + data.n_eq + i] -=
318 active_constraints[i] ? mu_in_neg : T(1);
319 }
320 return diff;
321}
322
323template<typename T>
331
342template<typename T, typename I, typename P>
343void
344qp_solve(Results<T>& results,
345 Model<T, I>& data,
346 const Settings<T>& settings,
347 Workspace<T, I>& work,
348 P& precond)
349{
350 if (settings.compute_timings) {
351 work.timer.stop();
352 work.timer.start();
353 }
354
355 if (work.internal
356 .dirty) // the following is used when a solve has already been executed
357 // (and without any intermediary model update)
358 {
359 proxsuite::linalg::sparse::MatMut<T, I> kkt_unscaled =
360 data.kkt_mut_unscaled();
361
362 auto kkt_top_n_rows = detail::top_rows_mut_unchecked(
363 proxsuite::linalg::veg::unsafe, kkt_unscaled, data.dim);
364
365 proxsuite::linalg::sparse::MatMut<T, I> H_unscaled =
366 detail::middle_cols_mut(kkt_top_n_rows, 0, data.dim, data.H_nnz);
367
368 proxsuite::linalg::sparse::MatMut<T, I> AT_unscaled =
369 detail::middle_cols_mut(kkt_top_n_rows, data.dim, data.n_eq, data.A_nnz);
370
371 proxsuite::linalg::sparse::MatMut<T, I> CT_unscaled =
372 detail::middle_cols_mut(
373 kkt_top_n_rows, data.dim + data.n_eq, data.n_in, data.C_nnz);
374
375 SparseMat<T, I> H_triu =
376 H_unscaled.to_eigen().template triangularView<Eigen::Upper>();
377 sparse::QpView<T, I> qp = {
378 { proxsuite::linalg::sparse::from_eigen, H_triu },
379 { proxsuite::linalg::sparse::from_eigen, data.g },
380 { proxsuite::linalg::sparse::from_eigen, AT_unscaled.to_eigen() },
381 { proxsuite::linalg::sparse::from_eigen, data.b },
382 { proxsuite::linalg::sparse::from_eigen, CT_unscaled.to_eigen() },
383 { proxsuite::linalg::sparse::from_eigen, data.l },
384 { proxsuite::linalg::sparse::from_eigen, data.u }
385 };
386
387 switch (settings.initial_guess) { // the following is used when one solve
388 // has already been executed
389 case InitialGuessStatus::EQUALITY_CONSTRAINED_INITIAL_GUESS: {
390 results.cleanup(settings);
391 break;
392 }
393 case InitialGuessStatus::COLD_START_WITH_PREVIOUS_RESULT: {
394 // keep solutions but restart workspace and results
395 results.cold_start(settings);
396 precond.scale_primal_in_place(
397 { proxsuite::proxqp::from_eigen, results.x });
398 precond.scale_dual_in_place_eq(
399 { proxsuite::proxqp::from_eigen, results.y });
400 precond.scale_dual_in_place_in(
401 { proxsuite::proxqp::from_eigen, results.z });
402 break;
403 }
404 case InitialGuessStatus::NO_INITIAL_GUESS: {
405 results.cleanup(settings);
406 break;
407 }
408 case InitialGuessStatus::WARM_START: {
409 results.cold_start(settings); // because there was already a solve,
410 // precond was already computed if set so
411 precond.scale_primal_in_place(
412 { proxsuite::proxqp::from_eigen,
413 results.x }); // it contains the value given in entry for warm start
414 precond.scale_dual_in_place_eq(
415 { proxsuite::proxqp::from_eigen, results.y });
416 precond.scale_dual_in_place_in(
417 { proxsuite::proxqp::from_eigen, results.z });
418 break;
419 }
420 case InitialGuessStatus::WARM_START_WITH_PREVIOUS_RESULT: {
421 // keep workspace and results solutions except statistics
422 results.cleanup_statistics();
423 precond.scale_primal_in_place(
424 { proxsuite::proxqp::from_eigen, results.x });
425 precond.scale_dual_in_place_eq(
426 { proxsuite::proxqp::from_eigen, results.y });
427 precond.scale_dual_in_place_in(
428 { proxsuite::proxqp::from_eigen, results.z });
429 break;
430 }
431 }
432 work.setup_impl(
433 qp,
434 data,
435 settings,
436 false,
437 precond,
438 P::scale_qp_in_place_req(
439 proxsuite::linalg::veg::Tag<T>{}, data.dim, data.n_eq, data.n_in));
440
441 } else {
442 // the following is used for a first solve after initializing or updating
443 // the Qp object
444 switch (settings.initial_guess) {
445 case InitialGuessStatus::EQUALITY_CONSTRAINED_INITIAL_GUESS: {
446 break;
447 }
448 case InitialGuessStatus::COLD_START_WITH_PREVIOUS_RESULT: {
449 precond.scale_primal_in_place(
450 { proxsuite::proxqp::from_eigen,
451 results.x }); // meaningful for when there is an upate of the model
452 // and one wants to warm start with previous result
453 precond.scale_dual_in_place_eq(
454 { proxsuite::proxqp::from_eigen, results.y });
455 precond.scale_dual_in_place_in(
456 { proxsuite::proxqp::from_eigen, results.z });
457 break;
458 }
459 case InitialGuessStatus::NO_INITIAL_GUESS: {
460 break;
461 }
462 case InitialGuessStatus::WARM_START: {
463 precond.scale_primal_in_place(
464 { proxsuite::proxqp::from_eigen, results.x });
465 precond.scale_dual_in_place_eq(
466 { proxsuite::proxqp::from_eigen, results.y });
467 precond.scale_dual_in_place_in(
468 { proxsuite::proxqp::from_eigen, results.z });
469 break;
470 }
471 case InitialGuessStatus::WARM_START_WITH_PREVIOUS_RESULT: {
472 precond.scale_primal_in_place(
473 { proxsuite::proxqp::from_eigen,
474 results.x }); // meaningful for when there is an upate of the model
475 // and one wants to warm start with previous result
476 precond.scale_dual_in_place_eq(
477 { proxsuite::proxqp::from_eigen, results.y });
478 precond.scale_dual_in_place_in(
479 { proxsuite::proxqp::from_eigen, results.z });
480 break;
481 }
482 }
483 }
484
485 if (settings.verbose) {
486 sparse::print_setup_header(settings, results, data);
487 }
488 using namespace proxsuite::linalg::veg::literals;
489 namespace util = proxsuite::linalg::sparse::util;
490 auto zx = util::zero_extend;
491
492 proxsuite::linalg::veg::dynstack::DynStackMut stack = work.stack_mut();
493
494 isize n = data.dim;
495 isize n_eq = data.n_eq;
496 isize n_in = data.n_in;
497 isize n_tot = n + n_eq + n_in;
498
499 VectorViewMut<T> x{ proxqp::from_eigen, results.x };
500 VectorViewMut<T> y{ proxqp::from_eigen, results.y };
501 VectorViewMut<T> z{ proxqp::from_eigen, results.z };
502
503 proxsuite::linalg::sparse::MatMut<T, I> kkt = data.kkt_mut();
504
505 auto kkt_top_n_rows =
506 detail::top_rows_mut_unchecked(proxsuite::linalg::veg::unsafe, kkt, n);
507
508 proxsuite::linalg::sparse::MatMut<T, I> H_scaled =
509 detail::middle_cols_mut(kkt_top_n_rows, 0, n, data.H_nnz);
510
511 proxsuite::linalg::sparse::MatMut<T, I> AT_scaled =
512 detail::middle_cols_mut(kkt_top_n_rows, n, n_eq, data.A_nnz);
513
514 proxsuite::linalg::sparse::MatMut<T, I> CT_scaled =
515 detail::middle_cols_mut(kkt_top_n_rows, n + n_eq, n_in, data.C_nnz);
516
517 auto& g_scaled_e = work.internal.g_scaled;
518 auto& b_scaled_e = work.internal.b_scaled;
519 auto& l_scaled_e = work.internal.l_scaled;
520 auto& u_scaled_e = work.internal.u_scaled;
521
522 QpViewMut<T, I> qp_scaled = {
523 H_scaled,
524 { proxsuite::linalg::sparse::from_eigen, g_scaled_e },
525 AT_scaled,
526 { proxsuite::linalg::sparse::from_eigen, b_scaled_e },
527 CT_scaled,
528 { proxsuite::linalg::sparse::from_eigen, l_scaled_e },
529 { proxsuite::linalg::sparse::from_eigen, u_scaled_e },
530 };
531
532 T const dual_feasibility_rhs_2 = infty_norm(data.g);
533
534 // auto ldl_col_ptrs = work.ldl_col_ptrs_mut();
535 auto ldl_col_ptrs = work.internal.ldl.col_ptrs.ptr_mut();
536 proxsuite::linalg::veg::Tag<I> itag;
537 proxsuite::linalg::veg::Tag<T> xtag;
538
539 bool do_ldlt = work.internal.do_ldlt;
540
541 isize ldlt_ntot = do_ldlt ? n_tot : 0;
542
543 auto _perm = stack.make_new_for_overwrite(itag, ldlt_ntot);
544
545 I* perm_inv = work.internal.ldl.perm_inv.ptr_mut();
546 I* perm = _perm.ptr_mut();
547
548 if (do_ldlt) {
549 // compute perm from perm_inv
550 for (isize i = 0; i < n_tot; ++i) {
551 perm[isize(zx(perm_inv[i]))] = I(i);
552 }
553 }
554
555 I* kkt_nnz_counts = work.internal.kkt_nnz_counts.ptr_mut();
556
557 auto& iterative_solver = *work.internal.matrix_free_solver.get();
558 isize C_active_nnz = 0;
559 switch (settings.initial_guess) {
560 case InitialGuessStatus::EQUALITY_CONSTRAINED_INITIAL_GUESS: {
561 // H and A are always active
562 for (usize j = 0; j < usize(n + n_eq); ++j) {
563 kkt_nnz_counts[isize(j)] = I(kkt.col_end(j) - kkt.col_start(j));
564 }
565 // ineq constraints initially inactive
566 for (isize j = 0; j < n_in; ++j) {
567 kkt_nnz_counts[n + n_eq + j] = 0;
568 work.active_inequalities[j] = false;
569 }
570 break;
571 }
572 case InitialGuessStatus::COLD_START_WITH_PREVIOUS_RESULT: {
573 // keep solutions + restart workspace and results except rho and mu : done
574 // in setup
575
576 // H and A are always active
577 for (usize j = 0; j < usize(n + n_eq); ++j) {
578 kkt_nnz_counts[isize(j)] = I(kkt.col_end(j) - kkt.col_start(j));
579 }
580 // keep constraints inactive from previous solution
581 for (isize j = 0; j < n_in; ++j) {
582 if (results.z(j) != 0) {
583 kkt_nnz_counts[n + n_eq + j] = I(kkt.col_end(usize(n + n_eq + j)) -
584 kkt.col_start(usize(n + n_eq + j)));
585 work.active_inequalities[j] = true;
586 C_active_nnz += kkt_nnz_counts[n + n_eq + j];
587 } else {
588 kkt_nnz_counts[n + n_eq + j] = 0;
589 work.active_inequalities[j] = false;
590 }
591 }
592 break;
593 }
594 case InitialGuessStatus::NO_INITIAL_GUESS: {
595 // already set to zero in the setup
596 // H and A are always active
597 for (usize j = 0; j < usize(n + n_eq); ++j) {
598 kkt_nnz_counts[isize(j)] = I(kkt.col_end(j) - kkt.col_start(j));
599 }
600 // ineq constraints initially inactive
601 for (isize j = 0; j < n_in; ++j) {
602 kkt_nnz_counts[n + n_eq + j] = 0;
603 work.active_inequalities[j] = false;
604 }
605 break;
606 }
607 case InitialGuessStatus::WARM_START: {
608 // keep previous solution
609
610 // H and A are always active
611 for (usize j = 0; j < usize(n + n_eq); ++j) {
612 kkt_nnz_counts[isize(j)] = I(kkt.col_end(j) - kkt.col_start(j));
613 }
614 // keep constraints inactive from previous solution
615 for (isize j = 0; j < n_in; ++j) {
616 if (results.z(j) != 0) {
617 kkt_nnz_counts[n + n_eq + j] = I(kkt.col_end(usize(n + n_eq + j)) -
618 kkt.col_start(usize(n + n_eq + j)));
619 work.active_inequalities[j] = true;
620 C_active_nnz += kkt_nnz_counts[n + n_eq + j];
621
622 } else {
623 kkt_nnz_counts[n + n_eq + j] = 0;
624 work.active_inequalities[j] = false;
625 }
626 }
627 break;
628 }
629 case InitialGuessStatus::WARM_START_WITH_PREVIOUS_RESULT: {
630 // keep workspace and results solutions except statistics
631 // H and A are always active
632 for (usize j = 0; j < usize(n + n_eq); ++j) {
633 kkt_nnz_counts[isize(j)] = I(kkt.col_end(j) - kkt.col_start(j));
634 }
635 // keep constraints inactive from previous solution
636 for (isize j = 0; j < n_in; ++j) {
637 if (results.z(j) != 0) {
638 kkt_nnz_counts[n + n_eq + j] = I(kkt.col_end(usize(n + n_eq + j)) -
639 kkt.col_start(usize(n + n_eq + j)));
640 work.active_inequalities[j] = true;
641 C_active_nnz += kkt_nnz_counts[n + n_eq + j];
642 } else {
643 kkt_nnz_counts[n + n_eq + j] = 0;
644 work.active_inequalities[j] = false;
645 }
646 }
647 break;
648 }
649 }
650
651 proxsuite::linalg::sparse::MatMut<T, I> kkt_active = {
652 proxsuite::linalg::sparse::from_raw_parts,
653 n_tot,
654 n_tot,
655 data.H_nnz + data.A_nnz + C_active_nnz,
656 kkt.col_ptrs_mut(),
657 kkt_nnz_counts,
658 kkt.row_indices_mut(),
659 kkt.values_mut(),
660 };
661
662 I* etree = work.internal.ldl.etree.ptr_mut();
663 I* ldl_nnz_counts = work.internal.ldl.nnz_counts.ptr_mut();
664 I* ldl_row_indices = work.internal.ldl.row_indices.ptr_mut();
665 T* ldl_values = work.internal.ldl.values.ptr_mut();
666 proxsuite::linalg::veg::SliceMut<bool> active_constraints =
667 work.active_inequalities.as_mut();
668
669 proxsuite::linalg::sparse::MatMut<T, I> ldl = {
670 proxsuite::linalg::sparse::from_raw_parts,
671 n_tot,
672 n_tot,
673 0,
674 ldl_col_ptrs,
675 do_ldlt ? ldl_nnz_counts : nullptr, // si do_ldlt est vrai do ldl_nnz_counts
676 ldl_row_indices,
677 ldl_values,
678 };
679
680 T bcl_eta_ext_init = pow(T(0.1), settings.alpha_bcl);
681 T bcl_eta_ext = bcl_eta_ext_init;
682 T bcl_eta_in(1);
683 T eps_in_min = std::min(settings.eps_abs, T(1e-9));
684
685 auto x_e = x.to_eigen();
686 auto y_e = y.to_eigen();
687 auto z_e = z.to_eigen();
688 sparse::refactorize<T, I>(
689 work, results, settings, kkt_active, active_constraints, data, stack, xtag);
690 switch (settings.initial_guess) {
691 case InitialGuessStatus::EQUALITY_CONSTRAINED_INITIAL_GUESS: {
692 LDLT_TEMP_VEC_UNINIT(T, rhs, n_tot, stack);
693 LDLT_TEMP_VEC_UNINIT(T, no_guess, 0, stack);
694
695 rhs.head(n) = -g_scaled_e;
696 rhs.segment(n, n_eq) = b_scaled_e;
697 rhs.segment(n + n_eq, n_in).setZero();
698
699 ldl_solve_in_place({ proxqp::from_eigen, rhs },
700 { proxqp::from_eigen, no_guess },
701 results,
702 data,
703 n_tot,
704 ldl,
705 iterative_solver,
706 do_ldlt,
707 stack,
708 ldl_values,
709 perm,
710 ldl_col_ptrs,
711 perm_inv,
712 settings,
713 kkt_active,
714 active_constraints);
715 x_e = rhs.head(n);
716 y_e = rhs.segment(n, n_eq);
717 z_e = rhs.segment(n + n_eq, n_in);
718 break;
719 }
720 case InitialGuessStatus::COLD_START_WITH_PREVIOUS_RESULT: {
721 // keep solutions but restart workspace and results
722 break;
723 }
724 case InitialGuessStatus::NO_INITIAL_GUESS: {
725 // already set to zero in the setup
726 break;
727 }
728 case InitialGuessStatus::WARM_START: {
729 // keep previous solution
730 break;
731 }
732 case InitialGuessStatus::WARM_START_WITH_PREVIOUS_RESULT: {
733 // keep workspace and results solutions except statistics
734 break;
735 }
736 }
737 T rhs_duality_gap(0);
738 T scaled_eps(settings.eps_abs);
739
740 for (isize iter = 0; iter < settings.max_iter; ++iter) {
741
742 results.info.iter_ext += 1;
743 if (iter == settings.max_iter) {
744 break;
745 }
746 T new_bcl_mu_eq = results.info.mu_eq;
747 T new_bcl_mu_in = results.info.mu_in;
748 T new_bcl_mu_eq_inv = results.info.mu_eq_inv;
749 T new_bcl_mu_in_inv = results.info.mu_in_inv;
750
751 {
752 T primal_feasibility_eq_rhs_0;
753 T primal_feasibility_in_rhs_0;
754
755 T dual_feasibility_rhs_0(0);
756 T dual_feasibility_rhs_1(0);
757 T dual_feasibility_rhs_3(0);
758
759 LDLT_TEMP_VEC_UNINIT(T, primal_residual_eq_scaled, n_eq, stack);
760 LDLT_TEMP_VEC_UNINIT(T, primal_residual_in_scaled_lo, n_in, stack);
761 LDLT_TEMP_VEC_UNINIT(T, primal_residual_in_scaled_up, n_in, stack);
762
763 LDLT_TEMP_VEC_UNINIT(T, dual_residual_scaled, n, stack);
764
765 // vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
766 auto is_primal_feasible = [&](T primal_feasibility_lhs) -> bool {
767 T rhs_pri = scaled_eps;
768 if (settings.eps_rel != 0) {
769 rhs_pri +=
770 settings.eps_rel * std::max({ primal_feasibility_eq_rhs_0,
771 primal_feasibility_in_rhs_0 });
772 }
773 return primal_feasibility_lhs <= rhs_pri;
774 };
775 auto is_dual_feasible = [&](T dual_feasibility_lhs) -> bool {
776 T rhs_dua = settings.eps_abs;
777 if (settings.eps_rel != 0) {
778 rhs_dua += settings.eps_rel * std::max({
779 dual_feasibility_rhs_0,
780 dual_feasibility_rhs_1,
781 dual_feasibility_rhs_2,
782 dual_feasibility_rhs_3,
783 });
784 }
785
786 return dual_feasibility_lhs <= rhs_dua;
787 };
788 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
789
790 VEG_BIND( // ?
791 auto,
792 (primal_feasibility_lhs, dual_feasibility_lhs),
793 detail::unscaled_primal_dual_residual(work,
794 results,
795 settings,
796 primal_residual_eq_scaled,
797 primal_residual_in_scaled_lo,
798 primal_residual_in_scaled_up,
799 dual_residual_scaled,
800 primal_feasibility_eq_rhs_0,
801 primal_feasibility_in_rhs_0,
802 dual_feasibility_rhs_0,
803 dual_feasibility_rhs_1,
804 dual_feasibility_rhs_3,
805 rhs_duality_gap,
806 precond,
807 data,
808 qp_scaled.as_const(),
809 detail::vec_mut(x_e),
810 detail::vec_mut(y_e),
811 detail::vec_mut(z_e),
812 stack));
813 /*put in debug mode
814 if (settings.verbose) {
815 std::cout << "-------- outer iteration: " << iter << " primal
816 residual "
817 << primal_feasibility_lhs
818 << " dual residual "
819 << dual_feasibility_lhs <<
820 " mu_in " << results.info.mu_in
821 << " bcl_eta_ext " <<
822 bcl_eta_ext << " bcl_eta_in "
823 << bcl_eta_in <<
824 std::endl;
825 }
826 */
827 if (settings.verbose) {
828 LDLT_TEMP_VEC_UNINIT(T, tmp, n, stack);
829 tmp.setZero();
830 detail::noalias_symhiv_add(tmp, qp_scaled.H.to_eigen(), x_e);
831 precond.unscale_dual_residual_in_place({ proxqp::from_eigen, tmp });
832
833 precond.unscale_primal_in_place({ proxqp::from_eigen, x_e });
834 precond.unscale_dual_in_place_eq({ proxqp::from_eigen, y_e });
835 precond.unscale_dual_in_place_in({ proxqp::from_eigen, z_e });
836 tmp *= 0.5;
837 tmp += data.g;
838 results.info.objValue = (tmp).dot(x_e);
839 std::cout << "\033[1;32m[outer iteration " << iter + 1 << "]\033[0m"
840 << std::endl;
841 std::cout << std::scientific << std::setw(2) << std::setprecision(2)
842 << "| primal residual=" << primal_feasibility_lhs
843 << " | dual residual=" << dual_feasibility_lhs
844 << " | duality gap=" << results.info.duality_gap
845 << " | mu_in=" << results.info.mu_in
846 << " | rho=" << results.info.rho << std::endl;
847 results.info.pri_res = primal_feasibility_lhs;
848 results.info.dua_res = dual_feasibility_lhs;
849 precond.scale_primal_in_place(VectorViewMut<T>{ from_eigen, x_e });
850 precond.scale_dual_in_place_eq(VectorViewMut<T>{ from_eigen, y_e });
851 precond.scale_dual_in_place_in(VectorViewMut<T>{ from_eigen, z_e });
852 }
853 if (is_primal_feasible(primal_feasibility_lhs) &&
854 is_dual_feasible(dual_feasibility_lhs)) {
855 if (settings.check_duality_gap) {
856 if (std::fabs(results.info.duality_gap) <=
857 settings.eps_duality_gap_abs +
858 settings.eps_duality_gap_rel * rhs_duality_gap) {
859 results.info.pri_res = primal_feasibility_lhs;
860 results.info.dua_res = dual_feasibility_lhs;
861 if (settings.primal_infeasibility_solving) {
862 results.info.status =
863 QPSolverOutput::PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE;
864 } else {
865 results.info.status = QPSolverOutput::PROXQP_SOLVED;
866 }
867 break;
868 }
869 } else {
870 results.info.pri_res = primal_feasibility_lhs;
871 results.info.dua_res = dual_feasibility_lhs;
872 if (settings.primal_infeasibility_solving) {
873 results.info.status =
874 QPSolverOutput::PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE;
875 } else {
876 results.info.status = QPSolverOutput::PROXQP_SOLVED;
877 }
878 break;
879 }
880 }
881
882 LDLT_TEMP_VEC_UNINIT(T, x_prev_e, n, stack);
883 LDLT_TEMP_VEC_UNINIT(T, y_prev_e, n_eq, stack);
884 LDLT_TEMP_VEC_UNINIT(T, z_prev_e, n_in, stack);
885 LDLT_TEMP_VEC(T, dw_prev, n_tot, stack);
886
887 x_prev_e = x_e;
888 y_prev_e = y_e;
889 z_prev_e = z_e;
890
891 // Cx + 1/mu_in * z_prev
892 primal_residual_in_scaled_up += results.info.mu_in * z_prev_e;
893 // switch (settings.merit_function_type) { NOT activated for the moment
894 // case MeritFunctionType::GPDAL:
895 // primal_residual_in_scaled_up +=
896 // (settings.alpha_gpdal - 1.) * results.info.mu_in * results.z;
897 // break;
898 // case MeritFunctionType::PDAL:
899 // break;
900 // }
901 primal_residual_in_scaled_lo = primal_residual_in_scaled_up;
902
903 // Cx - l + 1/mu_in * z_prev
904 primal_residual_in_scaled_lo -= l_scaled_e;
905
906 // Cx - u + 1/mu_in * z_prev
907 primal_residual_in_scaled_up -= u_scaled_e;
908
909 // vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
910 auto primal_dual_newton_semi_smooth = [&]() -> void {
911 for (isize iter_inner = 0; iter_inner < settings.max_iter_in;
912 ++iter_inner) {
913 LDLT_TEMP_VEC_UNINIT(T, dw, n_tot, stack);
914
915 if (iter_inner == settings.max_iter_in - 1) {
916 results.info.iter += settings.max_iter_in;
917 break;
918 }
919
920 // primal_dual_semi_smooth_newton_step
921 {
922 LDLT_TEMP_VEC_UNINIT(bool, new_active_constraints, n_in, stack);
923 auto rhs = dw;
924
925 work.active_set_low.array() =
926 primal_residual_in_scaled_lo.array() <= 0;
927 work.active_set_up.array() =
928 primal_residual_in_scaled_up.array() >= 0;
929 new_active_constraints = work.active_set_low || work.active_set_up;
930
931 // active set change
932 if (n_in > 0) {
933 bool removed = false;
934 bool added = false;
935
936 for (isize i = 0; i < n_in; ++i) {
937 bool was_active = active_constraints[i];
938 bool is_active = new_active_constraints[i];
939
940 isize idx = n + n_eq + i;
941
942 usize col_nnz =
943 zx(kkt.col_end(usize(idx))) - zx(kkt.col_start(usize(idx)));
944
945 if (is_active && !was_active) {
946 added = true;
947
948 kkt_active.nnz_per_col_mut()[idx] = I(col_nnz);
949 kkt_active._set_nnz(kkt_active.nnz() + isize(col_nnz));
950
951 if (do_ldlt) {
952 proxsuite::linalg::sparse::VecRef<T, I> new_col{
953 proxsuite::linalg::sparse::from_raw_parts,
954 n_tot,
955 isize(col_nnz),
956 kkt.row_indices() + zx(kkt.col_start(usize(idx))),
957 kkt.values() + zx(kkt.col_start(usize(idx))),
958 };
959 T mu_in_neg = -results.info.mu_in;
960 // switch (settings.merit_function_type)
961 // {
962 // case MeritFunctionType::GPDAL:
963 // mu_in_neg = -settings.alpha_gpdal * results.info.mu_in;
964 // break;
965 // case MeritFunctionType::PDAL:
966 // mu_in_neg = -results.info.mu_in;
967 // break;
968 // }
969 ldl = proxsuite::linalg::sparse::add_row(
970 ldl, etree, perm_inv, idx, new_col, mu_in_neg, stack);
971 }
972 active_constraints[i] = new_active_constraints[i];
973
974 } else if (!is_active && was_active) {
975 removed = true;
976 kkt_active.nnz_per_col_mut()[idx] = 0;
977 kkt_active._set_nnz(kkt_active.nnz() - isize(col_nnz));
978 if (do_ldlt) {
979 ldl = proxsuite::linalg::sparse::delete_row(
980 ldl, etree, perm_inv, idx, stack);
981 }
982 active_constraints[i] = new_active_constraints[i];
983 }
984 }
985
986 if (!do_ldlt) {
987 if (removed || added) {
988 refactorize(work,
989 results,
990 settings,
991 kkt_active,
992 active_constraints,
993 data,
994 stack,
995 xtag);
996 }
997 }
998 }
999 rhs.setZero();
1000 rhs.head(n) = -dual_residual_scaled;
1001 rhs.segment(n, n_eq) = -primal_residual_eq_scaled;
1002 for (isize i = 0; i < n_in; ++i) {
1003 if (work.active_set_up(i)) {
1004 rhs(n + n_eq + i) =
1005 results.info.mu_in * z_e(i) - primal_residual_in_scaled_up(i);
1006 } else if (work.active_set_low(i)) {
1007 rhs(n + n_eq + i) =
1008 results.info.mu_in * z_e(i) - primal_residual_in_scaled_lo(i);
1009 } else {
1010 rhs(n + n_eq + i) = -z_e(i);
1011 rhs.head(n) += z_e(i) * CT_scaled.to_eigen().col(i);
1012 }
1013 }
1014 // switch (settings.merit_function_type) {
1015 // case MeritFunctionType::GPDAL:
1016 // for (isize i = 0; i < n_in; ++i) {
1017 // if (work.active_set_up(i)) {
1018 // rhs(n + n_eq + i) =
1019 // -primal_residual_in_scaled_up(i) +
1020 // z_e(i) * results.info.mu_in * settings.alpha_gpdal;
1021 // } else if (work.active_set_low(i)) {
1022 // rhs(n + n_eq + i) =
1023 // -primal_residual_in_scaled_lo(i) +
1024 // z_e(i) * results.info.mu_in * settings.alpha_gpdal;
1025 // } else {
1026 // rhs(n + n_eq + i) = -z_e(i);
1027 // rhs.head(n) += z_e(i) * CT_scaled.to_eigen().col(i);
1028 // }
1029 // }
1030 // break;
1031 // case MeritFunctionType::PDAL:
1032 // for (isize i = 0; i < n_in; ++i) {
1033 // if (work.active_set_up(i)) {
1034 // rhs(n + n_eq + i) = results.info.mu_in * z_e(i) -
1035 // primal_residual_in_scaled_up(i);
1036 // } else if (work.active_set_low(i)) {
1037 // rhs(n + n_eq + i) = results.info.mu_in * z_e(i) -
1038 // primal_residual_in_scaled_lo(i);
1039 // } else {
1040 // rhs(n + n_eq + i) = -z_e(i);
1041 // rhs.head(n) += z_e(i) * CT_scaled.to_eigen().col(i);
1042 // }
1043 // }
1044 // break;
1045 // }
1046 ldl_solve_in_place(
1047 { proxqp::from_eigen, rhs },
1048 { proxqp::from_eigen,
1049 dw_prev }, // todo: MAJ dw_prev avec dw pour avoir meilleur
1050 // guess sur les solve in place
1051 results,
1052 data,
1053 n_tot,
1054 ldl,
1055 iterative_solver,
1056 do_ldlt,
1057 stack,
1058 ldl_values,
1059 perm,
1060 ldl_col_ptrs,
1061 perm_inv,
1062 settings,
1063 kkt_active,
1064 active_constraints);
1065 }
1066
1067 auto dx = dw.head(n);
1068 auto dy = dw.segment(n, n_eq);
1069 auto dz = dw.segment(n + n_eq, n_in);
1070 LDLT_TEMP_VEC(T, Hdx, n, stack);
1071 LDLT_TEMP_VEC(T, Adx, n_eq, stack);
1072 LDLT_TEMP_VEC(T, Cdx, n_in, stack);
1073
1074 LDLT_TEMP_VEC(T, ATdy, n, stack);
1075 LDLT_TEMP_VEC(T, CTdz, n, stack);
1076
1077 detail::noalias_symhiv_add(Hdx, H_scaled.to_eigen(), dx);
1078 detail::noalias_gevmmv_add(Adx, ATdy, AT_scaled.to_eigen(), dx, dy);
1079 detail::noalias_gevmmv_add(Cdx, CTdz, CT_scaled.to_eigen(), dx, dz);
1080 // switch (settings.merit_function_type) {
1081 // case MeritFunctionType::GPDAL:
1082 // Cdx.noalias() +=
1083 // (settings.alpha_gpdal - 1.) * results.info.mu_in * dz;
1084 // break;
1085 // case MeritFunctionType::PDAL:
1086 // break;
1087 // }
1088 T alpha = 1;
1089 // primal dual line search
1090 if (n_in > 0) {
1091 auto primal_dual_gradient_norm =
1092 [&](T alpha_cur) -> PrimalDualGradResult<T> {
1093 LDLT_TEMP_VEC_UNINIT(T, Cdx_active, n_in, stack);
1094 LDLT_TEMP_VEC_UNINIT(T, active_part_z, n_in, stack);
1095 {
1096 LDLT_TEMP_VEC_UNINIT(T, tmp_lo, n_in, stack);
1097 LDLT_TEMP_VEC_UNINIT(T, tmp_up, n_in, stack);
1098
1099 auto zero = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(n_in);
1100
1101 tmp_lo = primal_residual_in_scaled_lo + alpha_cur * Cdx;
1102 tmp_up = primal_residual_in_scaled_up + alpha_cur * Cdx;
1103 Cdx_active =
1104 (tmp_lo.array() < 0 || tmp_up.array() > 0).select(Cdx, zero);
1105 active_part_z = (tmp_lo.array() < 0)
1106 .select(primal_residual_in_scaled_lo, zero) +
1107 (tmp_up.array() > 0)
1108 .select(primal_residual_in_scaled_up, zero);
1109 }
1110
1111 T a = dx.dot(Hdx) + //
1112 results.info.rho * dx.squaredNorm() + //
1113 results.info.mu_eq_inv * Adx.squaredNorm() +
1114 +results.info.mu_in_inv * Cdx_active.squaredNorm() +
1115 results.info.nu * results.info.mu_eq *
1116 (results.info.mu_eq_inv * Adx - dy).squaredNorm() +
1117 results.info.nu * results.info.mu_in *
1118 (results.info.mu_in_inv * Cdx_active - dz).squaredNorm();
1119
1120 T b =
1121 x_e.dot(Hdx) + //
1122 (results.info.rho * (x_e - x_prev_e) + g_scaled_e).dot(dx) + //
1123 Adx.dot(results.info.mu_eq_inv * primal_residual_eq_scaled +
1124 y_e) + //
1125 results.info.mu_in_inv * Cdx_active.dot(active_part_z) + //
1126 results.info.nu * primal_residual_eq_scaled.dot(
1127 results.info.mu_eq_inv * Adx - dy) + //
1128 results.info.nu *
1129 (active_part_z - results.info.mu_in * z_e)
1130 .dot(results.info.mu_in_inv * Cdx_active - dz);
1131
1132 return {
1133 a,
1134 b,
1135 a * alpha_cur + b,
1136 };
1137 };
1138
1139 // auto gpdal_derivative_results =
1140 // [&](T alpha_cur) -> PrimalDualGradResult<T> {
1141 // LDLT_TEMP_VEC_UNINIT(T, Cdx_active, n_in, stack);
1142 // LDLT_TEMP_VEC_UNINIT(T, active_part_z, n_in, stack);
1143 // {
1144 // LDLT_TEMP_VEC_UNINIT(T, tmp_lo, n_in, stack);
1145 // LDLT_TEMP_VEC_UNINIT(T, tmp_up, n_in, stack);
1146
1147 // auto zero = Eigen::Matrix<T, Eigen::Dynamic, 1>::Zero(n_in);
1148
1149 // tmp_lo = primal_residual_in_scaled_lo + alpha_cur * Cdx;
1150 // tmp_up = primal_residual_in_scaled_up + alpha_cur * Cdx;
1151 // Cdx_active =
1152 // (tmp_lo.array() < 0 || tmp_up.array() > 0).select(Cdx,
1153 // zero);
1154 // active_part_z = (tmp_lo.array() < 0)
1155 // .select(primal_residual_in_scaled_lo, zero)
1156 // +
1157 // (tmp_up.array() > 0)
1158 // .select(primal_residual_in_scaled_up,
1159 // zero);
1160 // }
1161
1162 // T a = dx.dot(Hdx) + //
1163 // results.info.rho * dx.squaredNorm() + //
1164 // results.info.mu_eq_inv * Adx.squaredNorm() +
1165 // +results.info.mu_in_inv * Cdx_active.squaredNorm() /
1166 // settings.alpha_gpdal + results.info.mu_eq_inv *
1167 // (Adx - results.info.mu_eq *dy).squaredNorm() +
1168 // results.info.mu_in * (1. - settings.alpha_gpdal) *
1169 // (dz).squaredNorm();
1170
1171 // T b =
1172 // x_e.dot(Hdx) + // (results.info.rho * (x_e - x_prev_e) +
1173 // g_scaled_e).dot(dx) + // results.info.mu_eq_inv * Adx.dot(
1174 // primal_residual_eq_scaled +
1175 // y_e * results.info.mu_eq) + //
1176 // results.info.mu_in_inv * Cdx_active.dot(active_part_z) /
1177 // settings.alpha_gpdal + //
1178 // results.info.mu_eq_inv * primal_residual_eq_scaled.dot( Adx -
1179 // dy * results.info.mu_eq) + //
1180 // (z_e).dot(dz) * results.info.mu_in *
1181 // (1. - settings.alpha_gpdal);
1182
1183 // return {
1184 // a,
1185 // b,
1186 // a * alpha_cur + b,
1187 // };
1188 // };
1189
1190 LDLT_TEMP_VEC_UNINIT(T, alphas, 2 * n_in, stack);
1191 isize alphas_count = 0;
1192 const T machine_eps = std::numeric_limits<T>::epsilon();
1193
1194 for (isize i = 0; i < n_in; ++i) {
1195 T alpha_candidates[2] = {
1196 -primal_residual_in_scaled_lo(i) / (Cdx(i) + machine_eps),
1197 -primal_residual_in_scaled_up(i) / (Cdx(i) + machine_eps),
1198 };
1199
1200 for (auto alpha_candidate : alpha_candidates) {
1201 if (alpha_candidate > machine_eps) {
1202 alphas[alphas_count] = alpha_candidate;
1203 ++alphas_count;
1204 }
1205 }
1206 }
1207 std::sort(alphas.data(), alphas.data() + alphas_count);
1208 alphas_count =
1209 std::unique(alphas.data(), alphas.data() + alphas_count) -
1210 alphas.data();
1211 if (alphas_count > 0) { //&& alphas[0] <= 1
1212 auto infty = std::numeric_limits<T>::infinity();
1213
1214 T last_neg_grad = 0;
1215 T alpha_last_neg = 0;
1216 T first_pos_grad = 0;
1217 T alpha_first_pos = infty;
1218 {
1219 for (isize i = 0; i < alphas_count; ++i) {
1220 T alpha_cur = alphas[i];
1221 T gr = primal_dual_gradient_norm(alpha_cur).grad;
1222 // switch (settings.merit_function_type) {
1223 // case MeritFunctionType::GPDAL:
1224 // gr = gpdal_derivative_results(alpha_cur).grad;
1225 // break;
1226 // case MeritFunctionType::PDAL:
1227 // gr = primal_dual_gradient_norm(alpha_cur).grad;
1228 // break;
1229 // }
1230
1231 if (gr < 0) {
1232 alpha_last_neg = alpha_cur;
1233 last_neg_grad = gr;
1234 } else {
1235 first_pos_grad = gr;
1236 alpha_first_pos = alpha_cur;
1237 break;
1238 }
1239 }
1240
1241 if (alpha_last_neg == 0) {
1242 last_neg_grad =
1243 primal_dual_gradient_norm(alpha_last_neg).grad;
1244 // switch (settings.merit_function_type) {
1245 // case MeritFunctionType::GPDAL:
1246 // last_neg_grad =
1247 // gpdal_derivative_results(alpha_last_neg).grad;
1248 // break;
1249 // case MeritFunctionType::PDAL:
1250 // last_neg_grad =
1251 // primal_dual_gradient_norm(alpha_last_neg).grad;
1252 // break;
1253 // }
1254 }
1255 if (alpha_first_pos == infty) {
1256 auto res = primal_dual_gradient_norm(2 * alpha_last_neg + 1);
1257 alpha = -res.b / res.a;
1258 // switch (settings.merit_function_type) {
1259 // case MeritFunctionType::GPDAL: {
1260 // auto res =
1261 // gpdal_derivative_results(2 * alpha_last_neg + 1);
1262 // alpha = -res.b / res.a;
1263 // } break;
1264 // case MeritFunctionType::PDAL: {
1265 // auto res =
1266 // primal_dual_gradient_norm(2 * alpha_last_neg + 1);
1267 // alpha = -res.b / res.a;
1268 // } break;
1269 // }
1270 } else {
1271 alpha = alpha_last_neg -
1272 last_neg_grad * (alpha_first_pos - alpha_last_neg) /
1273 (first_pos_grad - last_neg_grad);
1274 }
1275 }
1276 } else {
1277 auto res = primal_dual_gradient_norm(T(0));
1278 alpha = -res.b / res.a;
1279 // switch (settings.merit_function_type) {
1280 // case MeritFunctionType::GPDAL: {
1281 // auto res = gpdal_derivative_results(T(0));
1282 // alpha = -res.b / res.a;
1283 // } break;
1284 // case MeritFunctionType::PDAL: {
1285 // auto res = primal_dual_gradient_norm(T(0));
1286 // alpha = -res.b / res.a;
1287 // } break;
1288 // }
1289 }
1290 }
1291 if (alpha * infty_norm(dw) < T(1e-11) && iter_inner > 0) {
1292 results.info.iter += iter_inner + 1;
1293 return;
1294 }
1295
1296 x_e += alpha * dx;
1297 y_e += alpha * dy;
1298 z_e += alpha * dz;
1299 dual_residual_scaled +=
1300 alpha * (Hdx + ATdy + CTdz + results.info.rho * dx);
1301 primal_residual_eq_scaled += alpha * (Adx - results.info.mu_eq * dy);
1302 primal_residual_in_scaled_lo += alpha * Cdx;
1303 primal_residual_in_scaled_up += alpha * Cdx;
1304
1305 T err_in = std::max({
1306 (infty_norm(helpers::negative_part(primal_residual_in_scaled_lo) +
1307 helpers::positive_part(primal_residual_in_scaled_up) -
1308 results.info.mu_in * z_e)),
1309 (infty_norm(primal_residual_eq_scaled)),
1310 (infty_norm(dual_residual_scaled)),
1311 });
1312 // switch (settings.merit_function_type) {
1313 // case MeritFunctionType::GPDAL:
1314 // err_in = std::max({
1315 // (infty_norm(
1316 // helpers::negative_part(primal_residual_in_scaled_lo) +
1317 // helpers::positive_part(primal_residual_in_scaled_up) -
1318 // settings.alpha_gpdal * results.info.mu_in * z_e)),
1319 // (infty_norm(primal_residual_eq_scaled)),
1320 // (infty_norm(dual_residual_scaled)),
1321 // });
1322 // break;
1323 // case MeritFunctionType::PDAL:
1324 // err_in = std::max({
1325 // (infty_norm(
1326 // helpers::negative_part(primal_residual_in_scaled_lo) +
1327 // helpers::positive_part(primal_residual_in_scaled_up) -
1328 // results.info.mu_in * z_e)),
1329 // (infty_norm(primal_residual_eq_scaled)),
1330 // (infty_norm(dual_residual_scaled)),
1331 // });
1332 // break;
1333 // }
1334 /* put in debug mode
1335 if (settings.verbose) {
1336 std::cout << "--inner iter " << iter_inner << " iner error "
1337 << err_in << " alpha "
1338 << alpha << " infty_norm(dw) "
1339 << infty_norm(dw) <<
1340 std::endl;
1341 }
1342 */
1343 if (settings.verbose) {
1344 std::cout << "\033[1;34m[inner iteration " << iter_inner + 1
1345 << "]\033[0m" << std::endl;
1346 std::cout << std::scientific << std::setw(2) << std::setprecision(2)
1347 << "| inner residual=" << err_in << " | alpha=" << alpha
1348 << std::endl;
1349 }
1350 if (iter_inner % settings.frequence_infeasibility_check == 0 ||
1351 settings.primal_infeasibility_solving) {
1352 // compute primal and dual infeasibility criteria
1353 bool is_primal_infeasible = proxsuite::proxqp::sparse::detail::
1354 global_primal_residual_infeasibility(
1355 VectorViewMut<T>{ from_eigen, ATdy },
1356 VectorViewMut<T>{ from_eigen, CTdz },
1357 VectorViewMut<T>{ from_eigen, dy },
1358 VectorViewMut<T>{ from_eigen, dz },
1359 qp_scaled.as_const(),
1360 settings,
1361 precond);
1362 bool is_dual_infeasible = proxsuite::proxqp::sparse::detail::
1363 global_dual_residual_infeasibility(
1364 VectorViewMut<T>{ from_eigen, Adx },
1365 VectorViewMut<T>{ from_eigen, Cdx },
1366 VectorViewMut<T>{ from_eigen, Hdx },
1367 VectorViewMut<T>{ from_eigen, dx },
1368 qp_scaled.as_const(),
1369 settings,
1370 data,
1371 precond);
1372 if (is_primal_infeasible) {
1373 results.info.status = QPSolverOutput::PROXQP_PRIMAL_INFEASIBLE;
1374 if (!settings.primal_infeasibility_solving) {
1375 results.info.iter += iter_inner + 1;
1376 dw_prev = dw;
1377 break;
1378 }
1379 } else if (is_dual_infeasible) {
1380 if (!settings.primal_infeasibility_solving) {
1381 results.info.status = QPSolverOutput::PROXQP_DUAL_INFEASIBLE;
1382 results.info.iter += iter_inner + 1;
1383 dw_prev = dw;
1384 break;
1385 }
1386 }
1387 }
1388 if (err_in <= bcl_eta_in) {
1389 results.info.iter += iter_inner + 1;
1390 return;
1391 }
1392 }
1393 };
1394 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1395
1396 primal_dual_newton_semi_smooth();
1397 if ((results.info.status == QPSolverOutput::PROXQP_PRIMAL_INFEASIBLE &&
1398 !settings.primal_infeasibility_solving) ||
1399 results.info.status == QPSolverOutput::PROXQP_DUAL_INFEASIBLE) {
1400 // certificate of infeasibility
1401 results.x = dw_prev.head(data.dim);
1402 results.y = dw_prev.segment(data.dim, data.n_eq);
1403 results.z = dw_prev.tail(data.n_in);
1404 // NOTE: values are unscaled at the end, so do the certificates of
1405 // infeasibility
1406 break;
1407 }
1408 if (scaled_eps == settings.eps_abs &&
1409 settings.primal_infeasibility_solving &&
1410 results.info.status == QPSolverOutput::PROXQP_PRIMAL_INFEASIBLE) {
1411 LDLT_TEMP_VEC(T, rhs_dim, n, stack);
1412 LDLT_TEMP_VEC(T, rhs_n_eq, n_eq, stack);
1413 LDLT_TEMP_VEC(T, rhs_n_in, n_in, stack);
1414 rhs_n_eq.setConstant(T(1));
1415 rhs_n_in.setConstant(T(1));
1416 rhs_dim.noalias() =
1417 AT_scaled.to_eigen() * rhs_n_eq + CT_scaled.to_eigen() * rhs_n_in;
1418 scaled_eps = infty_norm(rhs_dim) * settings.eps_abs;
1419 }
1420 // VEG bind : met le résultat tuple de unscaled_primal_dual_residual dans
1421 // (primal_feasibility_lhs_new, dual_feasibility_lhs_new) en guessant leur
1422 // type via auto
1423 VEG_BIND(
1424 auto,
1425 (primal_feasibility_lhs_new, dual_feasibility_lhs_new),
1426 detail::unscaled_primal_dual_residual(work,
1427 results,
1428 settings,
1429 primal_residual_eq_scaled,
1430 primal_residual_in_scaled_lo,
1431 primal_residual_in_scaled_up,
1432 dual_residual_scaled,
1433 primal_feasibility_eq_rhs_0,
1434 primal_feasibility_in_rhs_0,
1435 dual_feasibility_rhs_0,
1436 dual_feasibility_rhs_1,
1437 dual_feasibility_rhs_3,
1438 rhs_duality_gap,
1439 precond,
1440 data,
1441 qp_scaled.as_const(),
1442 detail::vec_mut(x_e),
1443 detail::vec_mut(y_e),
1444 detail::vec_mut(z_e),
1445 stack));
1446
1447 if (is_primal_feasible(primal_feasibility_lhs_new) &&
1448 is_dual_feasible(dual_feasibility_lhs_new)) {
1449 if (settings.check_duality_gap) {
1450 if (std::fabs(results.info.duality_gap) <=
1451 settings.eps_duality_gap_abs +
1452 settings.eps_duality_gap_rel * rhs_duality_gap) {
1453 results.info.pri_res = primal_feasibility_lhs_new;
1454 results.info.dua_res = dual_feasibility_lhs_new;
1455 if (settings.primal_infeasibility_solving &&
1456 results.info.status ==
1457 QPSolverOutput::PROXQP_PRIMAL_INFEASIBLE) {
1458 results.info.status =
1459 QPSolverOutput::PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE;
1460 } else {
1461 results.info.status = QPSolverOutput::PROXQP_SOLVED;
1462 }
1463 break;
1464 }
1465 } else {
1466 results.info.pri_res = primal_feasibility_lhs_new;
1467 results.info.dua_res = dual_feasibility_lhs_new;
1468 if (settings.primal_infeasibility_solving &&
1469 results.info.status == QPSolverOutput::PROXQP_PRIMAL_INFEASIBLE) {
1470 results.info.status =
1471 QPSolverOutput::PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE;
1472 } else {
1473 results.info.status = QPSolverOutput::PROXQP_SOLVED;
1474 }
1475 break;
1476 }
1477 }
1478
1479 // vvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvvv
1480 auto bcl_update = [&]() -> void {
1481 if (primal_feasibility_lhs_new <= bcl_eta_ext ||
1482 iter > settings.safe_guard) {
1483 bcl_eta_ext *= pow(results.info.mu_in, settings.beta_bcl);
1484 bcl_eta_in = std::max(bcl_eta_in * results.info.mu_in, eps_in_min);
1485
1486 } else {
1487 y_e = y_prev_e;
1488 z_e = z_prev_e;
1489 new_bcl_mu_in = std::max(
1490 results.info.mu_in * settings.mu_update_factor, settings.mu_min_in);
1491 new_bcl_mu_eq = std::max(
1492 results.info.mu_eq * settings.mu_update_factor, settings.mu_min_eq);
1493
1494 new_bcl_mu_in_inv =
1495 std::min(results.info.mu_in_inv * settings.mu_update_inv_factor,
1496 settings.mu_max_in_inv);
1497 new_bcl_mu_eq_inv =
1498 std::min(results.info.mu_eq_inv * settings.mu_update_inv_factor,
1499 settings.mu_max_eq_inv);
1500 bcl_eta_ext =
1501 bcl_eta_ext_init * pow(new_bcl_mu_in, settings.alpha_bcl);
1502 bcl_eta_in = std::max(new_bcl_mu_in, eps_in_min);
1503 }
1504 };
1505 // ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
1506 bcl_update();
1507
1508 VEG_BIND(
1509 auto,
1510 (_, dual_feasibility_lhs_new_2),
1511 detail::unscaled_primal_dual_residual(work,
1512 results,
1513 settings,
1514 primal_residual_eq_scaled,
1515 primal_residual_in_scaled_lo,
1516 primal_residual_in_scaled_up,
1517 dual_residual_scaled,
1518 primal_feasibility_eq_rhs_0,
1519 primal_feasibility_in_rhs_0,
1520 dual_feasibility_rhs_0,
1521 dual_feasibility_rhs_1,
1522 dual_feasibility_rhs_3,
1523 rhs_duality_gap,
1524 precond,
1525 data,
1526 qp_scaled.as_const(),
1527 detail::vec_mut(x_e),
1528 detail::vec_mut(y_e),
1529 detail::vec_mut(z_e),
1530 stack));
1531 proxsuite::linalg::veg::unused(_);
1532
1533 if (primal_feasibility_lhs_new >= primal_feasibility_lhs && //
1534 dual_feasibility_lhs_new_2 >= dual_feasibility_lhs && //
1535 results.info.mu_in <= T(1.E-5)) {
1536 new_bcl_mu_in = settings.cold_reset_mu_in;
1537 new_bcl_mu_eq = settings.cold_reset_mu_eq;
1538 new_bcl_mu_in_inv = settings.cold_reset_mu_in_inv;
1539 new_bcl_mu_eq_inv = settings.cold_reset_mu_eq_inv;
1540 }
1541 }
1542 if (results.info.mu_in != new_bcl_mu_in ||
1543 results.info.mu_eq != new_bcl_mu_eq) {
1544 {
1545 ++results.info.mu_updates;
1546 }
1547 /*
1548 refactorize(
1549 work,
1550 results,
1551 kkt_active,
1552 active_constraints,
1553 data,
1554 stack,
1555 xtag);
1556 */
1557 if (work.internal.do_ldlt) {
1558 isize w_values = 1; // un seul elt non nul
1559 T alpha = 0;
1560 for (isize j = 0; j < n_eq + n_in; ++j) {
1561 I row_index = I(j + n);
1562 if (j < n_eq) {
1563 alpha = results.info.mu_eq - new_bcl_mu_eq;
1564
1565 } else {
1566 if (!work.active_inequalities[j - n_eq]) {
1567 continue;
1568 }
1569 alpha = results.info.mu_in - new_bcl_mu_in;
1570 // switch (settings.merit_function_type)
1571 // {
1572 // case MeritFunctionType::GPDAL:
1573 // alpha = settings.alpha_gpdal * (results.info.mu_in -
1574 // new_bcl_mu_in); break;
1575 // case MeritFunctionType::PDAL:
1576 // alpha = results.info.mu_in - new_bcl_mu_in;
1577 // break;
1578 // }
1579 }
1580 T value = 1;
1581 proxsuite::linalg::sparse::VecRef<T, I> w{
1582 proxsuite::linalg::veg::from_raw_parts,
1583 n + n_eq + n_in,
1584 w_values,
1585 &row_index, // &: adresse de row index
1586 &value,
1587 };
1588 ldl = rank1_update(ldl, etree, perm_inv, w, alpha, stack);
1589 }
1590 } else {
1591 refactorize(work,
1592 results,
1593 settings,
1594 kkt_active,
1595 active_constraints,
1596 data,
1597 stack,
1598 xtag);
1599 }
1600 }
1601
1602 results.info.mu_eq = new_bcl_mu_eq;
1603 results.info.mu_in = new_bcl_mu_in;
1604 results.info.mu_eq_inv = new_bcl_mu_eq_inv;
1605 results.info.mu_in_inv = new_bcl_mu_in_inv;
1606 }
1607 LDLT_TEMP_VEC_UNINIT(T, tmp, n, stack);
1608 tmp.setZero();
1609 detail::noalias_symhiv_add(tmp, qp_scaled.H.to_eigen(), x_e);
1610 precond.unscale_dual_residual_in_place({ proxqp::from_eigen, tmp });
1611
1612 precond.unscale_primal_in_place({ proxqp::from_eigen, x_e });
1613 precond.unscale_dual_in_place_eq({ proxqp::from_eigen, y_e });
1614 precond.unscale_dual_in_place_in({ proxqp::from_eigen, z_e });
1615 tmp *= 0.5;
1616 tmp += data.g;
1617 results.info.objValue = (tmp).dot(x_e);
1618
1619 if (settings.compute_timings) {
1620 results.info.solve_time = work.timer.elapsed().user; // in nanoseconds
1621 results.info.run_time = results.info.solve_time + results.info.setup_time;
1622 }
1623 if (settings.verbose) {
1624 std::cout << "-------------------SOLVER STATISTICS-------------------"
1625 << std::endl;
1626 std::cout << "outer iter: " << results.info.iter_ext << std::endl;
1627 std::cout << "total iter: " << results.info.iter << std::endl;
1628 std::cout << "mu updates: " << results.info.mu_updates << std::endl;
1629 std::cout << "rho updates: " << results.info.rho_updates << std::endl;
1630 std::cout << "objective: " << results.info.objValue << std::endl;
1631 switch (results.info.status) {
1632 case QPSolverOutput::PROXQP_SOLVED: {
1633 std::cout << "status: "
1634 << "Solved" << std::endl;
1635 break;
1636 }
1637 case QPSolverOutput::PROXQP_MAX_ITER_REACHED: {
1638 std::cout << "status: "
1639 << "Maximum number of iterations reached" << std::endl;
1640 break;
1641 }
1642 case QPSolverOutput::PROXQP_PRIMAL_INFEASIBLE: {
1643 std::cout << "status: "
1644 << "Primal infeasible" << std::endl;
1645 break;
1646 }
1647 case QPSolverOutput::PROXQP_DUAL_INFEASIBLE: {
1648 std::cout << "status: "
1649 << "Dual infeasible" << std::endl;
1650 break;
1651 }
1652 case QPSolverOutput::PROXQP_SOLVED_CLOSEST_PRIMAL_FEASIBLE: {
1653 std::cout << "status: "
1654 << "Solved closest primal feasible" << std::endl;
1655 break;
1656 }
1657 case QPSolverOutput::PROXQP_NOT_RUN: {
1658 std::cout << "status: "
1659 << "Solver not run" << std::endl;
1660 break;
1661 }
1662 }
1663 if (settings.compute_timings)
1664 std::cout << "run time [μs]: " << results.info.solve_time << std::endl;
1665 std::cout << "--------------------------------------------------------"
1666 << std::endl;
1667 }
1668
1669 assert(!std::isnan(results.info.pri_res));
1670 assert(!std::isnan(results.info.dua_res));
1671 assert(!std::isnan(results.info.duality_gap));
1672
1673 work.set_dirty();
1674}
1675} // namespace sparse
1676} // namespace proxqp
1677} // namespace proxsuite
1678
1679#endif /* end of include guard PROXSUITE_PROXQP_SPARSE_SOLVER_HPP */
LDLT_TEMP_VEC
#define LDLT_TEMP_VEC(Type, Name, Rows, Stack)
Definition core.hpp:74
LDLT_TEMP_VEC_UNINIT
#define LDLT_TEMP_VEC_UNINIT(Type, Name, Rows, Stack)
Definition core.hpp:76
proxsuite::helpers::positive_part
auto positive_part(T const &expr) PROXSUITE_DEDUCE_RET((expr.array() > 0).select(expr
Returns the positive part of an expression.
proxsuite::helpers::negative_part
auto negative_part(T const &expr) PROXSUITE_DEDUCE_RET((expr.array()< 0).select(expr
Returns the negative part of an expression.
proxsuite::linalg::sparse::delete_row
auto delete_row(MatMut< T, I > ld, I *etree, I const *perm_inv, isize pos, DynStackMut stack) noexcept(false) -> MatMut< T, I >
Definition rowmod.hpp:52
proxsuite::linalg::sparse::add_row
auto add_row(MatMut< T, I > ld, I *etree, I const *perm_inv, isize pos, VecRef< T, I > new_col, proxsuite::linalg::veg::DoNotDeduce< T > diag_element, DynStackMut stack) noexcept(false) -> MatMut< T, I >
Definition rowmod.hpp:190
proxsuite::proxqp::sparse::detail::noalias_gevmmv_add
void noalias_gevmmv_add(OutL &&out_l, OutR &&out_r, A const &a, InL const &in_l, InR const &in_r)
Definition utils.hpp:261
proxsuite::proxqp::sparse::detail::global_primal_residual_infeasibility
bool global_primal_residual_infeasibility(VectorViewMut< T > ATdy, VectorViewMut< T > CTdz, VectorViewMut< T > dy, VectorViewMut< T > dz, const QpView< T, I > qp_scaled, const Settings< T > &qpsettings, const P &ruiz)
Definition utils.hpp:473
proxsuite::proxqp::sparse::detail::noalias_symhiv_add
void noalias_symhiv_add(Out &&out, A const &a, In const &in)
Definition utils.hpp:278
proxsuite::proxqp::sparse::detail::middle_cols_mut
auto middle_cols_mut(proxsuite::linalg::sparse::MatMut< T, I > mat, isize start, isize ncols, isize nnz) -> proxsuite::linalg::sparse::MatMut< T, I >
Definition utils.hpp:400
proxsuite::proxqp::sparse::detail::top_rows_mut_unchecked
auto top_rows_mut_unchecked(proxsuite::linalg::veg::Unsafe, proxsuite::linalg::sparse::MatMut< T, I > mat, isize nrows) -> proxsuite::linalg::sparse::MatMut< T, I >
Definition utils.hpp:440
proxsuite::proxqp::sparse::detail::vec_mut
auto vec_mut(V &&v) -> VecMapMut< typename proxsuite::linalg::veg::uncvref_t< V >::Scalar >
Definition utils.hpp:360
proxsuite::proxqp::sparse::detail::global_dual_residual_infeasibility
bool global_dual_residual_infeasibility(VectorViewMut< T > Adx, VectorViewMut< T > Cdx, VectorViewMut< T > Hdx, VectorViewMut< T > dx, const QpView< T, I > qp_scaled, const Settings< T > &qpsettings, const Model< T, I > &qpmodel, const P &ruiz)
Definition utils.hpp:529
proxsuite::proxqp::sparse::detail::unscaled_primal_dual_residual
auto unscaled_primal_dual_residual(Workspace< T, I > &work, Results< T > &results, const Settings< T > &settings, VecMapMut< T > primal_residual_eq_scaled, VecMapMut< T > primal_residual_in_scaled_lo, VecMapMut< T > primal_residual_in_scaled_up, VecMapMut< T > dual_residual_scaled, T &primal_feasibility_eq_rhs_0, T &primal_feasibility_in_rhs_0, T &dual_feasibility_rhs_0, T &dual_feasibility_rhs_1, T &dual_feasibility_rhs_3, T &rhs_duality_gap, const P &precond, Model< T, I > const &data, const QpView< T, I > qp_scaled, VecMapMut< T > x_e, VecMapMut< T > y_e, VecMapMut< T > z_e, proxsuite::linalg::veg::dynstack::DynStackMut stack) -> proxsuite::linalg::veg::Tuple< T, T >
Definition utils.hpp:616
proxsuite::proxqp::sparse::ldl_iter_solve_noalias
void ldl_iter_solve_noalias(VectorViewMut< T > sol, VectorView< T > rhs, VectorView< T > init_guess, Results< T > const &results, Model< T, I > const &data, isize n_tot, proxsuite::linalg::sparse::MatMut< T, I > ldl, Eigen::MINRES< detail::AugmentedKkt< T, I >, Eigen::Upper|Eigen::Lower, Eigen::IdentityPreconditioner > &iterative_solver, bool do_ldlt, proxsuite::linalg::veg::dynstack::DynStackMut stack, T *ldl_values, I *perm, I *ldl_col_ptrs, I const *perm_inv, Settings< T > const &settings, proxsuite::linalg::sparse::MatMut< T, I > kkt_active, proxsuite::linalg::veg::SliceMut< bool > active_constraints)
Definition solver.hpp:88
proxsuite::proxqp::sparse::refactorize
void refactorize(Workspace< T, I > &work, Results< T > const &results, Settings< T > const &settings, proxsuite::linalg::sparse::MatMut< T, I > kkt_active, proxsuite::linalg::veg::SliceMut< bool > active_constraints, Model< T, I > const &data, proxsuite::linalg::veg::dynstack::DynStackMut stack, proxsuite::linalg::veg::Tag< T > &xtag)
Definition workspace.hpp:33
proxsuite::proxqp::sparse::inner_reconstructed_matrix
auto inner_reconstructed_matrix(proxsuite::linalg::sparse::MatMut< T, I > ldl) -> DMat< T >
Definition solver.hpp:250
proxsuite::proxqp::sparse::ldl_solve_in_place
void ldl_solve_in_place(VectorViewMut< T > rhs, VectorView< T > init_guess, Results< T > const &results, Model< T, I > const &data, isize n_tot, proxsuite::linalg::sparse::MatMut< T, I > ldl, Eigen::MINRES< detail::AugmentedKkt< T, I >, Eigen::Upper|Eigen::Lower, Eigen::IdentityPreconditioner > &iterative_solver, bool do_ldlt, proxsuite::linalg::veg::dynstack::DynStackMut stack, T *ldl_values, I *perm, I *ldl_col_ptrs, I const *perm_inv, Settings< T > const &settings, proxsuite::linalg::sparse::MatMut< T, I > kkt_active, proxsuite::linalg::veg::SliceMut< bool > active_constraints)
Definition solver.hpp:201
proxsuite::proxqp::sparse::qp_solve
void qp_solve(Results< T > &results, Model< T, I > &data, const Settings< T > &settings, Workspace< T, I > &work, P &precond)
Definition solver.hpp:344
proxsuite::proxqp::sparse::SparseMat
Eigen::SparseMatrix< T, Eigen::ColMajor, I > SparseMat
Definition fwd.hpp:31
proxsuite::proxqp::sparse::reconstructed_matrix
auto reconstructed_matrix(proxsuite::linalg::sparse::MatMut< T, I > ldl, I const *perm_inv, isize n_tot) -> DMat< T >
Definition solver.hpp:271
proxsuite::proxqp::sparse::print_setup_header
void print_setup_header(const Settings< T > &settings, Results< T > &results, const Model< T, I > &model)
Definition utils.hpp:36
proxsuite::proxqp::sparse::ldl_solve
void ldl_solve(VectorViewMut< T > sol, VectorView< T > rhs, isize n_tot, proxsuite::linalg::sparse::MatMut< T, I > ldl, Eigen::MINRES< detail::AugmentedKkt< T, I >, Eigen::Upper|Eigen::Lower, Eigen::IdentityPreconditioner > &iterative_solver, bool do_ldlt, proxsuite::linalg::veg::dynstack::DynStackMut stack, T *ldl_values, I *perm, I *ldl_col_ptrs, I const *perm_inv)
Definition solver.hpp:40
proxsuite::proxqp::sparse::DMat
Eigen::Matrix< T, -1, -1 > DMat
Definition fwd.hpp:23
proxsuite::proxqp::sparse::reconstruction_error
auto reconstruction_error(proxsuite::linalg::sparse::MatMut< T, I > ldl, I const *perm_inv, Results< T > const &results, Model< T, I > const &data, isize n_tot, proxsuite::linalg::sparse::MatMut< T, I > kkt_active, proxsuite::linalg::veg::SliceMut< bool > active_constraints) -> DMat< T >
Definition solver.hpp:299
proxsuite::proxqp::usize
decltype(sizeof(0)) usize
Definition views.hpp:24
proxsuite::linalg::veg::_detail::_collections::VecImpl::ptr_mut
VEG_NODISCARD VEG_INLINE auto ptr_mut() VEG_NOEXCEPT -> T *
Definition vec.hpp:966
proxsuite::linalg::veg::_detail::_collections::VecImpl::as_mut
VEG_NODISCARD VEG_INLINE auto as_mut() VEG_NOEXCEPT -> SliceMut< T >
Definition vec.hpp:957
proxsuite::proxqp::Results
This class stores all the results of PROXQP solvers with sparse and dense backends.
Definition results.hpp:68
proxsuite::proxqp::Results::z
sparse::Vec< T > z
Definition results.hpp:74
proxsuite::proxqp::Results::cold_start
void cold_start(optional< Settings< T > > settings=nullopt)
Definition results.hpp:175
proxsuite::proxqp::Results::cleanup
void cleanup(optional< Settings< T > > settings=nullopt)
Definition results.hpp:149
proxsuite::proxqp::Results::y
sparse::Vec< T > y
Definition results.hpp:73
proxsuite::proxqp::Results::x
sparse::Vec< T > x
Definition results.hpp:72
proxsuite::proxqp::Settings
This class defines the settings of PROXQP solvers with sparse and dense backends.
Definition settings.hpp:89
proxsuite::proxqp::Settings::initial_guess
InitialGuessStatus initial_guess
Definition settings.hpp:123
proxsuite::proxqp::VectorViewMut::to_eigen
VEG_INLINE auto to_eigen() const -> detail::VecMapMut< T >
Definition views.hpp:677
proxsuite::proxqp::VectorViewMut::as_const
VEG_INLINE auto as_const() const noexcept -> VectorView< T >
Definition views.hpp:656
proxsuite::proxqp::VectorView::to_eigen
VEG_INLINE auto to_eigen() const -> detail::VecMap< T >
Definition views.hpp:627
proxsuite::proxqp::sparse::Model
This class stores the model of the QP problem.
Definition model.hpp:23
proxsuite::proxqp::sparse::Model::kkt_mut
auto kkt_mut() -> proxsuite::linalg::sparse::MatMut< T, I >
Definition model.hpp:96
proxsuite::proxqp::sparse::Model::kkt_mut_unscaled
auto kkt_mut_unscaled() -> proxsuite::linalg::sparse::MatMut< T, I >
Definition model.hpp:134
proxsuite::proxqp::sparse::PrimalDualGradResult::VEG_REFLECT
VEG_REFLECT(PrimalDualGradResult, a, b, grad)
proxsuite::proxqp::sparse::Workspace
This class defines the workspace of the sparse solver.
Definition workspace.hpp:120
proxsuite::proxqp::sparse::Workspace::u_scaled
Eigen::Matrix< T, Eigen::Dynamic, 1 > u_scaled
Definition workspace.hpp:136
proxsuite::proxqp::sparse::Workspace::active_inequalities
proxsuite::linalg::veg::Vec< bool > active_inequalities
Definition workspace.hpp:165
proxsuite::proxqp::sparse::Workspace::b_scaled
Eigen::Matrix< T, Eigen::Dynamic, 1 > b_scaled
Definition workspace.hpp:134
proxsuite::proxqp::sparse::Workspace::setup_impl
void setup_impl(const QpView< T, I > qp, Model< T, I > &data, const Settings< T > &settings, bool execute_or_not, P &precond, proxsuite::linalg::veg::dynstack::StackReq precond_req)
Definition workspace.hpp:333
proxsuite::proxqp::sparse::Workspace::matrix_free_solver
std::unique_ptr< Eigen::MINRES< detail::AugmentedKkt< T, I >, Eigen::Upper|Eigen::Lower, Eigen::IdentityPreconditioner > > matrix_free_solver
Definition workspace.hpp:146
proxsuite::proxqp::sparse::Workspace::l_scaled
Eigen::Matrix< T, Eigen::Dynamic, 1 > l_scaled
Definition workspace.hpp:135
proxsuite::proxqp::sparse::Workspace::g_scaled
Eigen::Matrix< T, Eigen::Dynamic, 1 > g_scaled
Definition workspace.hpp:133
proxsuite::proxqp::sparse::Workspace::internal
struct proxsuite::proxqp::sparse::Workspace::@16 internal
proxsuite::proxqp::sparse::Workspace::kkt_nnz_counts
proxsuite::linalg::veg::Vec< I > kkt_nnz_counts
Definition workspace.hpp:137
proxsuite::proxqp::sparse::Workspace::stack_mut
auto stack_mut() -> proxsuite::linalg::veg::dynstack::DynStackMut
Definition workspace.hpp:788
VEG_BIND
#define VEG_BIND(CV_Auto, Identifiers, Tuple)
Definition tuple.hpp:53