proxsuite 0.6.7
The Advanced Proximal Optimization Toolbox
Loading...
Searching...
No Matches
workspace.hpp
Go to the documentation of this file.
1//
2// Copyright (c) 2022 INRIA
3//
7#ifndef PROXSUITE_PROXQP_DENSE_WORKSPACE_HPP
8#define PROXSUITE_PROXQP_DENSE_WORKSPACE_HPP
9
10#include <Eigen/Core>
11#include <proxsuite/linalg/dense/ldlt.hpp>
12#include <proxsuite/proxqp/timings.hpp>
13#include <proxsuite/linalg/veg/vec.hpp>
14#include <proxsuite/proxqp/settings.hpp>
15namespace proxsuite {
16namespace proxqp {
17namespace dense {
21
24template<typename T>
25struct Workspace
26{
27
29 proxsuite::linalg::dense::Ldlt<T> ldl{};
30 proxsuite::linalg::veg::Vec<unsigned char> ldl_stack;
31 Timer<T> timer;
32
34 Mat<T> H_scaled;
35 Vec<T> g_scaled;
36 Mat<T> A_scaled;
37 Mat<T> C_scaled;
38 Vec<T> b_scaled;
39 Vec<T> u_scaled;
40 Vec<T> l_scaled;
41
42 Vec<T> u_box_scaled;
43 Vec<T> l_box_scaled;
44 Vec<T> i_scaled;
45
47
48 Vec<T> x_prev;
49 Vec<T> y_prev;
50 Vec<T> z_prev;
51
53 Mat<T> kkt;
54
56 VecISize current_bijection_map;
57 VecISize new_bijection_map;
58
59 VecBool active_set_up;
60 VecBool active_set_low;
61 VecBool active_inequalities;
62
64
65 Vec<T> Hdx;
66 Vec<T> Cdx;
67 Vec<T> Adx;
68
69 Vec<T> active_part_z;
70 proxsuite::linalg::veg::Vec<T> alphas;
71
73 Vec<T> dw_aug;
74 Vec<T> rhs;
75 Vec<T> err;
76
78
79 T dual_feasibility_rhs_2;
80 T correction_guess_rhs_g;
81 T correction_guess_rhs_b;
82 T alpha;
83
84 Vec<T> dual_residual_scaled;
85 Vec<T> primal_residual_in_scaled_up;
86
87 Vec<T> primal_residual_in_scaled_up_plus_alphaCdx;
88 Vec<T> primal_residual_in_scaled_low_plus_alphaCdx;
89 Vec<T> CTz;
90
91 bool constraints_changed;
92 bool dirty;
93 bool refactorize;
94 bool proximal_parameter_update;
95 bool is_initialized;
96
97 sparse::isize n_c; // final number of active inequalities
104 Workspace(isize dim = 0,
105 isize n_eq = 0,
106 isize n_in = 0,
107 bool box_constraints = false,
108 DenseBackend dense_backend = DenseBackend::PrimalDualLDLT)
109 : ldl{}
110 , H_scaled(dim, dim)
111 , g_scaled(dim)
112 , A_scaled(n_eq, dim)
113 , C_scaled(n_in, dim)
114 , b_scaled(n_eq)
115 , u_scaled(n_in)
116 , l_scaled(n_in)
117 , x_prev(dim)
118 , y_prev(n_eq)
119 , Hdx(dim)
120 , Adx(n_eq)
121 , dual_residual_scaled(dim)
122 , CTz(dim)
123 , constraints_changed(false)
124 , dirty(false)
125 , refactorize(false)
126 , proximal_parameter_update(false)
127 , is_initialized(false)
128 {
129
130 if (box_constraints) {
131 u_box_scaled.resize(dim);
132 u_box_scaled.setZero();
133 l_box_scaled.resize(dim);
134 l_box_scaled.setZero();
135 i_scaled.resize(dim);
136 i_scaled.setOnes();
137
138 z_prev.resize(dim + n_in);
139 // TODO appropriate heuristic for automatic choice
140 switch (dense_backend) {
141 case DenseBackend::PrimalDualLDLT:
142 kkt.resize(dim + n_eq, dim + n_eq);
143 ldl.reserve_uninit(dim + n_eq + n_in + dim);
144 ldl_stack.resize_for_overwrite(
145 proxsuite::linalg::veg::dynstack::StackReq(
146 // optimize here
147 proxsuite::linalg::dense::Ldlt<T>::factorize_req(dim + n_eq +
148 n_in + dim) |
149
150 (proxsuite::linalg::dense::temp_vec_req(
151 proxsuite::linalg::veg::Tag<T>{}, n_eq + n_in + dim) &
152 proxsuite::linalg::veg::dynstack::StackReq{
153 isize{ sizeof(isize) } * (n_eq + n_in + dim),
154 alignof(isize) } &
155 proxsuite::linalg::dense::Ldlt<T>::diagonal_update_req(
156 dim + n_eq + n_in + dim, n_eq + n_in + dim)) |
157
158 (proxsuite::linalg::dense::temp_mat_req(
159 proxsuite::linalg::veg::Tag<T>{},
160 dim + n_eq + n_in + dim,
161 n_in + dim) &
162 proxsuite::linalg::dense::Ldlt<T>::insert_block_at_req(
163 dim + n_eq + n_in + dim, n_in + dim)) |
164
165 proxsuite::linalg::dense::Ldlt<T>::solve_in_place_req(dim + n_eq +
166 n_in + dim))
167 // TODO optimize here
168 .alloc_req());
169 break;
170 case DenseBackend::PrimalLDLT:
171 kkt.resize(dim, dim);
172 ldl.reserve_uninit(dim);
173 ldl_stack.resize_for_overwrite(
174 proxsuite::linalg::veg::dynstack::StackReq(
175
176 proxsuite::linalg::dense::Ldlt<T>::factorize_req(dim) |
177 // check simplification possible
178 (proxsuite::linalg::dense::temp_vec_req(
179 proxsuite::linalg::veg::Tag<T>{}, n_eq + n_in + dim) &
180 proxsuite::linalg::veg::dynstack::StackReq{
181 isize{ sizeof(isize) } * (n_eq + n_in + dim),
182 alignof(isize) } &
183 proxsuite::linalg::dense::Ldlt<T>::diagonal_update_req(
184 dim + n_eq + n_in + dim, n_eq + n_in + dim)) |
185
186 (proxsuite::linalg::dense::temp_mat_req(
187 proxsuite::linalg::veg::Tag<T>{},
188 dim + n_eq + n_in + dim,
189 n_in + dim) &
190 proxsuite::linalg::dense::Ldlt<T>::insert_block_at_req(
191 dim + n_eq + n_in + dim, n_in + dim)) |
192 // end check
193 proxsuite::linalg::dense::Ldlt<T>::solve_in_place_req(dim))
194
195 .alloc_req());
196 break;
197 case DenseBackend::Automatic:
198 break;
199 }
200 current_bijection_map.resize(n_in + dim);
201 new_bijection_map.resize(n_in + dim);
202 for (isize i = 0; i < n_in + dim; i++) {
203 current_bijection_map(i) = i;
204 new_bijection_map(i) = i;
205 }
206 active_set_up.resize(n_in + dim);
207 active_set_low.resize(n_in + dim);
208 active_inequalities.resize(n_in + dim);
209 active_part_z.resize(n_in + dim);
210 dw_aug.resize(dim + n_eq + n_in + dim);
211 rhs.resize(dim + n_eq + n_in + dim);
212 err.resize(dim + n_eq + n_in + dim);
213 primal_residual_in_scaled_up.resize(dim + n_in);
214 primal_residual_in_scaled_up_plus_alphaCdx.resize(dim + n_in);
215 primal_residual_in_scaled_low_plus_alphaCdx.resize(dim + n_in);
216 Cdx.resize(n_in + dim);
217 alphas.reserve(2 * n_in + 2 * dim);
218 } else {
219 z_prev.resize(n_in);
220
221 switch (dense_backend) {
222 case DenseBackend::PrimalDualLDLT:
223 kkt.resize(dim + n_eq, dim + n_eq);
224 ldl.reserve_uninit(dim + n_eq + n_in);
225 ldl_stack.resize_for_overwrite(
226 proxsuite::linalg::veg::dynstack::StackReq(
227 // todo optimize here
228 proxsuite::linalg::dense::Ldlt<T>::factorize_req(dim + n_eq +
229 n_in) |
230
231 (proxsuite::linalg::dense::temp_vec_req(
232 proxsuite::linalg::veg::Tag<T>{}, n_eq + n_in) &
233 proxsuite::linalg::veg::dynstack::StackReq{
234 isize{ sizeof(isize) } * (n_eq + n_in), alignof(isize) } &
235 proxsuite::linalg::dense::Ldlt<T>::diagonal_update_req(
236 dim + n_eq + n_in, n_eq + n_in)) |
237
238 (proxsuite::linalg::dense::temp_mat_req(
239 proxsuite::linalg::veg::Tag<T>{}, dim + n_eq + n_in, n_in) &
240 proxsuite::linalg::dense::Ldlt<T>::insert_block_at_req(
241 dim + n_eq + n_in, n_in)) |
242
243 proxsuite::linalg::dense::Ldlt<T>::solve_in_place_req(dim + n_eq +
244 n_in))
245 // end todo optimize here
246 .alloc_req());
247 break;
248 case DenseBackend::PrimalLDLT:
249 kkt.resize(dim, dim);
250 ldl.reserve_uninit(dim);
251 ldl_stack.resize_for_overwrite(
252 proxsuite::linalg::veg::dynstack::StackReq(
253
254 proxsuite::linalg::dense::Ldlt<T>::factorize_req(dim) |
255 // check if it can be more simplified
256 (proxsuite::linalg::dense::temp_vec_req(
257 proxsuite::linalg::veg::Tag<T>{}, n_eq + n_in) &
258 proxsuite::linalg::veg::dynstack::StackReq{
259 isize{ sizeof(isize) } * (n_eq + n_in), alignof(isize) } &
260 proxsuite::linalg::dense::Ldlt<T>::diagonal_update_req(
261 dim + n_eq + n_in, n_eq + n_in)) |
262 (proxsuite::linalg::dense::temp_mat_req(
263 proxsuite::linalg::veg::Tag<T>{}, dim + n_eq + n_in, n_in) &
264 proxsuite::linalg::dense::Ldlt<T>::insert_block_at_req(
265 dim + n_eq + n_in, n_in)) |
266 // end check
267 proxsuite::linalg::dense::Ldlt<T>::solve_in_place_req(dim))
268
269 .alloc_req());
270 break;
271 case DenseBackend::Automatic:
272 break;
273 }
274
275 current_bijection_map.resize(n_in);
276 new_bijection_map.resize(n_in);
277 for (isize i = 0; i < n_in; i++) {
278 current_bijection_map(i) = i;
279 new_bijection_map(i) = i;
280 }
281 active_set_up.resize(n_in);
282 active_set_low.resize(n_in);
283 active_inequalities.resize(n_in);
284 active_part_z.resize(n_in);
285 dw_aug.resize(dim + n_eq + n_in);
286 rhs.resize(dim + n_eq + n_in);
287 err.resize(dim + n_eq + n_in);
288 primal_residual_in_scaled_up.resize(n_in);
289 primal_residual_in_scaled_up_plus_alphaCdx.resize(n_in);
290 primal_residual_in_scaled_low_plus_alphaCdx.resize(n_in);
291 Cdx.resize(n_in);
292 alphas.reserve(2 * n_in);
293 }
294
295 H_scaled.setZero();
296 g_scaled.setZero();
297 A_scaled.setZero();
298 C_scaled.setZero();
299 b_scaled.setZero();
300 u_scaled.setZero();
301 l_scaled.setZero();
302 x_prev.setZero();
303 y_prev.setZero();
304 z_prev.setZero();
305 kkt.setZero();
306 Hdx.setZero();
307 Cdx.setZero();
308 Adx.setZero();
309 active_part_z.setZero();
310 dw_aug.setZero();
311 rhs.setZero();
312 err.setZero();
313
314 dual_feasibility_rhs_2 = 0;
315 correction_guess_rhs_g = 0;
316 correction_guess_rhs_b = 0;
317 alpha = 1.;
318
319 dual_residual_scaled.setZero();
320 primal_residual_in_scaled_up.setZero();
321
322 primal_residual_in_scaled_up_plus_alphaCdx.setZero();
323 primal_residual_in_scaled_low_plus_alphaCdx.setZero();
324 CTz.setZero();
325 n_c = 0;
326 }
330 void cleanup(const bool box_constraints)
331 {
332 isize n_in = C_scaled.rows();
333 isize dim = H_scaled.rows();
334 H_scaled.setZero();
335 g_scaled.setZero();
336 A_scaled.setZero();
337 C_scaled.setZero();
338 b_scaled.setZero();
339 u_scaled.setZero();
340 l_scaled.setZero();
341 Hdx.setZero();
342 Cdx.setZero();
343 Adx.setZero();
344 active_part_z.setZero();
345 dw_aug.setZero();
346 rhs.setZero();
347 err.setZero();
348
349 alpha = 1.;
350
351 dual_residual_scaled.setZero();
352 primal_residual_in_scaled_up.setZero();
353
354 primal_residual_in_scaled_up_plus_alphaCdx.setZero();
355 primal_residual_in_scaled_low_plus_alphaCdx.setZero();
356 CTz.setZero();
357
358 x_prev.setZero();
359 y_prev.setZero();
360 z_prev.setZero();
361 isize n_constraints(n_in);
362 if (box_constraints) {
363 n_constraints += dim;
364 }
365 for (isize i = 0; i < n_constraints; i++) {
366 current_bijection_map(i) = i;
367 new_bijection_map(i) = i;
368 active_inequalities(i) = false;
369 }
370
371 constraints_changed = false;
372 dirty = false;
373 refactorize = false;
374 proximal_parameter_update = false;
375 is_initialized = false;
376 n_c = 0;
377 }
378};
379} // namespace dense
380} // namespace proxqp
381} // namespace proxsuite
382
383#endif /* end of include guard PROXSUITE_PROXQP_DENSE_WORKSPACE_HPP */
proxsuite::linalg::dense::temp_vec_req
auto temp_vec_req(proxsuite::linalg::veg::Tag< T >, isize rows) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
Definition core.hpp:862
proxsuite::linalg::dense::temp_mat_req
auto temp_mat_req(proxsuite::linalg::veg::Tag< T >, isize rows, isize cols) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
Definition core.hpp:850
proxsuite::proxqp::dense::VecISize
Eigen::Matrix< isize, DYN, 1 > VecISize
Definition fwd.hpp:48
proxsuite::proxqp::dense::VecBool
Eigen::Matrix< bool, DYN, 1 > VecBool
Definition fwd.hpp:51
proxsuite::proxqp::dense::Vec
Eigen::Matrix< T, DYN, 1 > Vec
Definition fwd.hpp:24
proxsuite::proxqp::dense::Mat
Eigen::Matrix< T, DYN, DYN, l > Mat
Definition fwd.hpp:31
proxsuite::linalg::dense::Ldlt::diagonal_update_req
static auto diagonal_update_req(isize n, isize r) noexcept -> proxsuite::linalg::veg::dynstack::StackReq
Definition ldlt.hpp:484
proxsuite::proxqp::Timer
This class mimics the way "boost/timer/timer.hpp" operates while using the modern std::chrono library...
Definition timings.hpp:36
proxsuite::proxqp::dense::Workspace
This class defines the workspace of the dense solver.
Definition workspace.hpp:26
proxsuite::proxqp::dense::Workspace::ldl_stack
proxsuite::linalg::veg::Vec< unsigned char > ldl_stack
Definition workspace.hpp:30
proxsuite::proxqp::dense::Workspace::cleanup
void cleanup(const bool box_constraints)
Definition workspace.hpp:330
proxsuite::proxqp::dense::Workspace::Workspace
Workspace(isize dim=0, isize n_eq=0, isize n_in=0, bool box_constraints=false, DenseBackend dense_backend=DenseBackend::PrimalDualLDLT)
Definition workspace.hpp:104
proxsuite::proxqp::dense::Workspace::alphas
proxsuite::linalg::veg::Vec< T > alphas
Definition workspace.hpp:70
proxsuite::proxqp::dense::Workspace::ldl
proxsuite::linalg::dense::Ldlt< T > ldl
Definition workspace.hpp:29