Line data Source code
1 : /*
2 : * The MIT License (MIT)
3 : *
4 : * Copyright (c) 2017 Pantelis Sopasakis (https://alphaville.github.io),
5 : * Krina Menounou (https://www.linkedin.com/in/krinamenounou),
6 : * Panagiotis Patrinos (http://homes.esat.kuleuven.be/~ppatrino)
7 : * Copyright (c) 2012 Brendan O'Donoghue (bodonoghue85@gmail.com)
8 : *
9 : * Permission is hereby granted, free of charge, to any person obtaining a copy
10 : * of this software and associated documentation files (the "Software"), to deal
11 : * in the Software without restriction, including without limitation the rights
12 : * to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
13 : * copies of the Software, and to permit persons to whom the Software is
14 : * furnished to do so, subject to the following conditions:
15 : *
16 : * The above copyright notice and this permission notice shall be included in all
17 : * copies or substantial portions of the Software.
18 : *
19 : * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
20 : * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
21 : * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
22 : * AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
23 : * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
24 : * OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
25 : * SOFTWARE.
26 : *
27 : */
28 : #include "scs.h"
29 : #include "normalize.h"
30 : #include "directions.h"
31 : #include "linsys/amatrix.h"
32 : #include <time.h>
33 :
34 : /* if verbose print summary output every this num iterations */
35 : #define SCS_PRINT_INTERVAL 100
36 : /* check for convergence every this num iterations */
37 : #define SCS_CONVERGED_INTERVAL 1
38 :
39 : /* tolerance at which we declare problem indeterminate */
40 : #define SCS_INDETERMINATE_TOL 1e-9
41 :
42 : /**
43 : * \brief Structure to hold residual information (unnormalized)
44 : */
45 : struct scs_residuals {
46 : /**
47 : * \brief The last iteration when the residuals were updated.
48 : */
49 : scs_int last_iter;
50 : /**
51 : * \brief Dual residual
52 : *
53 : * \f[
54 : * \text{resdual} = \frac{E(A'y + \tau c)}{\tau(1+\|c\|)\cdot \text{scale}_c\cdot \text{scale}}
55 : * \f]
56 : */
57 : scs_float res_dual;
58 : /**
59 : * \brief Primal residual
60 : *
61 : * \f[
62 : * \text{respri} = \frac{\|D(Ax+s-\tau b)\|}{\tau(1+\|b\|)(\text{scale}_b\cdot \text{scale})}
63 : * \f]
64 : */
65 : scs_float res_pri;
66 : /**
67 : * \brief Infeasibility residual
68 : *
69 : * \f[
70 : * \text{infres} = -\frac{\|Db\| \|EA'y\|}{b'y \cdot \text{scale}}
71 : * \f]
72 : */
73 : scs_float res_infeas;
74 : /**
75 : * \brief Unboundedness
76 : *
77 : * \f[
78 : * \text{unbdd} = -\frac{\|Ec\| \|D(Ax+s)}{c'x\cdot \text{scale}}
79 : * \f]
80 : */
81 : scs_float res_unbdd;
82 : /**
83 : * \brief Relative duality gap defined as
84 : *
85 : * \f[
86 : * \text{relgap} = \frac{c'x + b'y}{1+|c'x|+|b'y|}
87 : * \f]
88 : */
89 : scs_float rel_gap;
90 : scs_float cTx_by_tau; /* not divided by tau */
91 : scs_float bTy_by_tau; /* not divided by tau */
92 : /**
93 : * Variable \f$\tau\f$ (\f$\bar{\tau}\f$ in SuperSCS)
94 : */
95 : scs_float tau; /* for superSCS it's tau_b */
96 : /**
97 : * Variable \f$\kappa\f$ (\f$\bar{\kappa}\f$ in SuperSCS)
98 : */
99 : scs_float kap; /* for superSCS it's kap_b */
100 : };
101 : /* printing header */
102 : static const char *SCS_HEADER[] = {
103 : " Iter ", " pri res ", " dua res ", " rel gap ",
104 : " pri obj ", " dua obj ", " kap/tau ", " FPR ", " time (s)"
105 : };
106 : static const scs_int scs_hspace = 9;
107 : static const scs_int scs_header_length = 9;
108 : static const scs_int scs_header_line_length = 87;
109 :
110 : static scs_int scs_isnan(scs_float x) {
111 42 : return (isnan(x)); /* `isnan` works both for `float` and `double` types */
112 : }
113 :
114 321 : static ScsDirectionCache * scs_init_direction_cache(scs_int memory, scs_int l, scs_int print_mode, ScsDirectionType dir_type) {
115 321 : ScsDirectionCache * RESTRICT cache = scs_calloc(1, sizeof (*cache));
116 321 : scs_int length_S = 0, length_U = 0, length_S_minus_U = 0, length_t = 0, length_ws = 0;
117 :
118 :
119 321 : if (cache == SCS_NULL) {
120 : /* LCOV_EXCL_START */
121 : scs_special_print(print_mode, stderr, "ERROR: allocating YSCache failure\n");
122 : return SCS_NULL;
123 : /* LCOV_EXCL_STOP */
124 : }
125 :
126 321 : cache->ls_wspace_length = 0;
127 321 : cache->current_mem = 0;
128 :
129 321 : switch (dir_type) {
130 : case restarted_broyden:
131 : /* we allocate one extra memory position because it's needed */
132 318 : length_S = (1 + memory) * l;
133 318 : length_U = (1 + memory) * l;
134 318 : break;
135 : case anderson_acceleration:
136 3 : length_S = memory * l;
137 3 : length_U = memory * l;
138 3 : length_S_minus_U = memory * l;
139 : /* -----------------------------------------------------------------
140 : * Although t is a vector of dimension 'memory', we allocate space
141 : * of length 'l' because we first need to store the residual 'R'
142 : * therein (see directions.c: scs_compute_dir_anderson).
143 : * Note: we give some additional space because the BLAS documentation
144 : * says that we should give **at least** this much memory.
145 : * ----------------------------------------------------------------- */
146 3 : length_t = l;
147 : #ifdef SVD_ACTIVATED
148 3 : cache->ls_wspace_length = 1000 + scs_svd_workspace_size(l, memory);
149 : /* -----------------------------------------------------------------
150 : * In Anderson's acceleration, we solve a least squares problem
151 : * using lapack's SVD-based 'sgelss' (see svdls). This requires a
152 : * workspace whose size is given by 'cache->ls_wspace_length' (above).
153 : * To that we add another 'memory' memory positions to store the
154 : * singular values of the LHS matrix and another 'l*memory'
155 : * memory positions to make a copy of the Y-cache; this is
156 : * because 'sgelss' modified the LHS (i.e., the Y-cache).
157 : * ----------------------------------------------------------------- */
158 3 : length_ws = (cache->ls_wspace_length) + memory * (1 + l);
159 : #else
160 : cache->ls_wspace_length = 1000 + scs_qr_workspace_size(l, memory);
161 : /* -----------------------------------------------------------------
162 : * In Anderson's acceleration, we solve a least squares problem
163 : * using lapack's SVD-based 'sgelss' (see svdls). This requires a
164 : * workspace whose size is given by 'cache->ls_wspace_length' (above).
165 : * To that we add another 'l*memory'
166 : * memory positions to make a copy of the Y-cache; this is
167 : * because 'sgelss' modified the LHS (i.e., the Y-cache).
168 : * ----------------------------------------------------------------- */
169 : length_ws = (cache->ls_wspace_length) + memory * l;
170 : #endif
171 3 : break;
172 : default:
173 : break;
174 : }
175 :
176 321 : cache->S = scs_malloc(length_S * sizeof (scs_float));
177 321 : cache->U = scs_malloc(length_U * sizeof (scs_float));
178 321 : cache->S_minus_Y = scs_malloc(length_S_minus_U * sizeof (scs_float));
179 321 : cache->t = scs_malloc(length_t * sizeof (scs_float));
180 321 : cache->ls_wspace = scs_malloc(length_ws * sizeof (scs_float));
181 :
182 : /* the cache must know its memory length */
183 321 : cache->mem = memory;
184 :
185 : /* initial active memory is 0 */
186 321 : scs_reset_direction_cache(cache);
187 321 : return cache;
188 : }
189 :
190 326 : static void scs_free_direction_cache(ScsDirectionCache * RESTRICT cache) {
191 326 : if (cache == SCS_NULL)
192 326 : return;
193 321 : scs_free(cache->S);
194 321 : scs_free(cache->U);
195 321 : scs_free(cache->S_minus_Y);
196 321 : scs_free(cache->t);
197 321 : scs_free(cache->ls_wspace);
198 321 : scs_free(cache);
199 : }
200 :
201 329 : static void scs_free_work(ScsWork * RESTRICT work) {
202 658 : if (work == SCS_NULL) return;
203 329 : scs_free(work->u);
204 329 : scs_free(work->v);
205 329 : scs_free(work->u_t);
206 329 : scs_free(work->u_prev);
207 329 : scs_free(work->h);
208 329 : scs_free(work->g);
209 329 : scs_free(work->b);
210 329 : scs_free(work->c);
211 329 : scs_free(work->pr);
212 329 : scs_free(work->dr);
213 329 : if (work->scal != SCS_NULL) {
214 326 : scs_free(work->scal->D);
215 326 : scs_free(work->scal->E);
216 326 : scs_free(work->scal);
217 : }
218 329 : scs_free(work->u_b);
219 :
220 329 : if (work->stgs->do_super_scs == 1) {
221 326 : scs_free(work->R);
222 326 : scs_free(work->R_prev);
223 326 : scs_free(work->dir);
224 326 : scs_free(work->dut);
225 326 : scs_free(work->Sk);
226 326 : scs_free(work->Yk);
227 326 : scs_free(work->wu);
228 326 : scs_free(work->wu_t);
229 326 : scs_free(work->wu_b);
230 326 : scs_free(work->Rwu);
231 326 : scs_free_direction_cache(work->direction_cache);
232 326 : scs_free(work->s_b);
233 326 : scs_free(work->H);
234 : }
235 329 : scs_free(work);
236 : }
237 :
238 : /* LCOV_EXCL_START */
239 : static void scs_print_init_header(
240 : const ScsData * RESTRICT data,
241 : const ScsCone * RESTRICT cone) {
242 : scs_int i;
243 : ScsSettings * RESTRICT stgs = data->stgs;
244 : char *RESTRICT coneStr = scs_get_cone_header(cone);
245 : char *RESTRICT linSysMethod = scs_get_linsys_method(data->A, data->stgs);
246 : FILE * RESTRICT stream = stgs->output_stream;
247 : scs_int print_mode = stgs->do_override_streams;
248 : for (i = 0; i < scs_header_line_length; ++i) {
249 : scs_special_print(print_mode, stream, "-");
250 : }
251 : scs_special_print(print_mode, stream, "\n\tSCS v%s - Superlinear Splitting Conic Solver (SuperSCS)\n\t"
252 : "Web: https://kul-forbes.github.io/scs\n\t"
253 : "(c) P. Sopasakis, K. Menounou, P. Patrinos, KU Leuven, 2017-8\n\t"
254 : "(c) Brendan O'Donoghue, Stanford University, 2012-2016\n",
255 : scs_version());
256 : for (i = 0; i < scs_header_line_length; ++i) {
257 : scs_special_print(print_mode, stream, "-");
258 : }
259 : scs_special_print(print_mode, stream, "\n");
260 : if (linSysMethod) {
261 : scs_special_print(print_mode, stream, "Lin-sys: %s\n", linSysMethod);
262 : scs_free(linSysMethod);
263 : }
264 : if (stgs->normalize) {
265 : scs_special_print(print_mode, stream, "eps = %.2e, alpha = %.2f, max_iters = %i, normalize = %i, "
266 : "scale = %2.2f\n",
267 : stgs->eps, stgs->alpha, (int) stgs->max_iters,
268 : (int) stgs->normalize, stgs->scale);
269 : } else {
270 : scs_special_print(print_mode, stream, "eps = %.2e, alpha = %.2f, max_iters = %i, normalize = %i\n",
271 : stgs->eps, stgs->alpha, (int) stgs->max_iters,
272 : (int) stgs->normalize);
273 : }
274 : scs_special_print(print_mode, stream, "do_super_scs = %i, direction = %i, "
275 : "memory = %i\n", (int) stgs->do_super_scs, (int) stgs->direction,
276 : (int) stgs->memory);
277 : scs_special_print(print_mode, stream, "Variables n = %i, constraints m = %i\n", (int) data->n, (int) data->m);
278 : scs_special_print(print_mode, stream, "%s", coneStr);
279 : scs_free(coneStr);
280 : #ifdef MATLAB_MEX_FILE
281 : mexEvalString("drawnow;");
282 : #endif
283 : }
284 :
285 : /* LCOV_EXCL_STOP */
286 :
287 12 : static void scs_populate_on_failure(
288 : scs_int m,
289 : scs_int n,
290 : ScsSolution * RESTRICT sol,
291 : ScsInfo * RESTRICT info,
292 : scs_int statusVal,
293 : const char *RESTRICT msg) {
294 12 : if (info) {
295 12 : info->relGap = NAN;
296 12 : info->resPri = NAN;
297 12 : info->resDual = NAN;
298 12 : info->pobj = NAN;
299 12 : info->dobj = NAN;
300 12 : info->iter = -1;
301 12 : info->statusVal = statusVal;
302 12 : info->solveTime = NAN;
303 12 : strncpy(info->status, msg, SCS_INFO_STATUS_MSG_LENGTH);
304 : }
305 12 : if (sol) {
306 12 : if (n > 0) {
307 12 : if (sol->x == SCS_NULL)
308 1 : sol->x = scs_malloc(sizeof (scs_float) * n);
309 12 : scs_scale_array(sol->x, NAN, n);
310 : }
311 12 : if (m > 0) {
312 12 : if (sol->y == SCS_NULL)
313 1 : sol->y = scs_malloc(sizeof (scs_float) * m);
314 12 : scs_scale_array(sol->y, NAN, m);
315 12 : if (sol->s == SCS_NULL)
316 1 : sol->s = scs_malloc(sizeof (scs_float) * m);
317 12 : scs_scale_array(sol->s, NAN, m);
318 : }
319 : }
320 12 : }
321 :
322 12 : static scs_int scs_failure(
323 : ScsWork * RESTRICT work,
324 : scs_int m,
325 : scs_int n,
326 : ScsSolution * RESTRICT sol,
327 : ScsInfo * RESTRICT info,
328 : scs_int stint,
329 : const char *RESTRICT msg,
330 : const char *RESTRICT ststr,
331 : scs_int print_mode) {
332 12 : scs_int status = stint;
333 12 : scs_populate_on_failure(m, n, sol, info, status, ststr);
334 12 : scs_special_print(print_mode, stderr, "Failure:%s\n", msg);
335 12 : endInterruptListener();
336 12 : return status;
337 : }
338 :
339 2 : static void scs_warm_start_vars(
340 : ScsWork * RESTRICT work,
341 : const ScsSolution * RESTRICT sol) {
342 : scs_int i;
343 2 : scs_int n = work->n;
344 2 : scs_int m = work->m;
345 2 : scs_float * RESTRICT Ax = SCS_NULL;
346 2 : scs_float * RESTRICT ATy = SCS_NULL;
347 :
348 2 : if (work->stgs->do_super_scs == 0) {
349 1 : memset(work->v, 0, n * sizeof (scs_float));
350 1 : memcpy(work->u, sol->x, n * sizeof (scs_float));
351 1 : memcpy(&(work->u[n]), sol->y, m * sizeof (scs_float));
352 1 : memcpy(&(work->v[n]), sol->s, m * sizeof (scs_float));
353 1 : work->v[n + m] = 0.0;
354 1 : work->u[n + m] = 1.0;
355 : } else {
356 1 : memcpy(work->u_t, sol->x, n * sizeof (scs_float));
357 1 : memcpy(&(work->u_t[n]), sol->y, m * sizeof (scs_float));
358 1 : work->u_t[n + m] = 1.0;
359 : }
360 : #ifndef NOVALIDATE
361 2 : if (work->stgs->do_super_scs == 0) {
362 17 : for (i = 0; i < n + m + 1; ++i) {
363 34 : if (scs_isnan(work->u[i])) {
364 0 : work->u[i] = 0;
365 : }
366 34 : if (scs_isnan(work->v[i])) {
367 0 : work->v[i] = 0;
368 : }
369 : }
370 : } else {
371 8 : for (i = 0; i < n + m + 1; ++i) {
372 16 : if (scs_isnan(work->u_t[i])) {
373 0 : work->u_t[i] = 0;
374 : }
375 : }
376 : }
377 : #endif
378 2 : if (work->stgs->normalize) {
379 1 : scs_normalize_warm_start(work);
380 : }
381 2 : if (work->stgs->do_super_scs) {
382 1 : Ax = scs_calloc(m, sizeof (scs_float));
383 1 : ATy = scs_calloc(n, sizeof (scs_float));
384 :
385 1 : scs_accum_by_a(work->A, work->p, work->u_t, Ax); /* Ax_t = A*x_t */
386 1 : scs_accum_by_a_trans(work->A, work->p, &(work->u_t[n]), ATy); /* ATy_t = AT*y_t */
387 4 : for (i = 0; i < n; ++i) {
388 : /* rho_x*x_t + ATy_t + c*tau_t */
389 3 : work->u[i] = work->u_t[i] + ATy[i] + work->c[i] * work->u_t[n + m];
390 : }
391 4 : for (i = 0; i < m; ++i) {
392 : /* -Ax_t + y_t + b*tau_t */
393 4 : work->u[i + n] = -Ax[i] + work->u_t[i + n] + work->b[i] * work->u_t[n + m];
394 : }
395 : /* -cTx_t - BTy_t + tau_t */
396 2 : work->u[n + m] = -scs_inner_product(work->c, work->u_t, work->n)
397 1 : - scs_inner_product(work->b, &(work->u_t[n]), work->m)
398 1 : + work->u_t[n + m];
399 : }
400 2 : if (Ax != SCS_NULL) {
401 1 : scs_free(Ax);
402 : }
403 2 : if (ATy != SCS_NULL) {
404 1 : scs_free(ATy);
405 : }
406 2 : }
407 :
408 201 : static scs_float scs_calc_primal_resid(
409 : ScsWork * RESTRICT work,
410 : const scs_float * RESTRICT x,
411 : const scs_float * RESTRICT s,
412 : const scs_float tau,
413 : scs_float * RESTRICT nmAxs) {
414 : scs_int i;
415 201 : scs_float pres = 0, scale, *RESTRICT pr = work->pr;
416 201 : *nmAxs = 0;
417 201 : memset(pr, 0, work->m * sizeof (scs_float));
418 201 : scs_accum_by_a(work->A, work->p, x, pr);
419 201 : scs_add_scaled_array(pr, s, work->m, 1.0); /* pr = Ax + s */
420 1941208 : for (i = 0; i < work->m; ++i) {
421 1941007 : scale =
422 1941007 : work->stgs->normalize ? work->scal->D[i] / (work->sc_b * work->stgs->scale) : 1;
423 1941007 : scale = scale * scale;
424 1941007 : *nmAxs += (pr[i] * pr[i]) * scale;
425 1941007 : pres += (pr[i] - work->b[i] * tau) * (pr[i] - work->b[i] * tau) * scale;
426 : }
427 201 : *nmAxs = SQRTF(*nmAxs);
428 201 : return SQRTF(pres); /* norm(Ax + s - b * tau), for superSCS: norm(Ax_b + s_b - b * tau_b) */
429 : }
430 :
431 201 : static scs_float scs_calc_dual_resid(
432 : ScsWork * RESTRICT work,
433 : const scs_float * RESTRICT y,
434 : const scs_float tau,
435 : scs_float * RESTRICT nmATy) {
436 : scs_int i;
437 201 : scs_float dres = 0, scale, *dr = work->dr;
438 201 : *nmATy = 0;
439 201 : memset(dr, 0, work->n * sizeof (scs_float));
440 201 : scs_accum_by_a_trans(work->A, work->p, y, dr); /* dr = A'y */
441 970668 : for (i = 0; i < work->n; ++i) {
442 970467 : scale =
443 970467 : work->stgs->normalize ? (work->scal->E[i] / (work->sc_c * work->stgs->scale)) : 1;
444 970467 : scale = scale * scale;
445 970467 : *nmATy += (dr[i] * dr[i]) * scale;
446 970467 : dres += (dr[i] + work->c[i] * tau) * (dr[i] + work->c[i] * tau) * scale;
447 : }
448 201 : *nmATy = SQRTF(*nmATy);
449 201 : return SQRTF(dres); /* norm(A'y + c * tau), for superSCS: norm(A'y_b + c * tau_b)*/
450 : }
451 :
452 : /* calculates un-normalized quantities */
453 203 : static void scs_calc_residuals(
454 : ScsWork * RESTRICT work,
455 : struct scs_residuals * RESTRICT res,
456 : scs_int iter) {
457 : scs_float * RESTRICT x;
458 : scs_float * RESTRICT y;
459 : scs_float * RESTRICT s;
460 : scs_float nmpr_tau;
461 : scs_float nmdr_tau;
462 : scs_float nmAxs_tau;
463 : scs_float nmATy_tau;
464 : scs_float cTx, bTy;
465 203 : scs_int n = work->n, m = work->m;
466 :
467 : /* checks if the residuals are unchanged by checking iteration */
468 203 : if (res->last_iter == iter) {
469 2 : return;
470 : }
471 201 : res->last_iter = iter;
472 :
473 201 : s = &(work->v[work->n]);
474 201 : x = work->u;
475 201 : y = &(work->u[work->n]);
476 :
477 201 : res->tau = ABS(work->u[n + m]);
478 402 : res->kap = ABS(work->v[n + m]) /
479 201 : (work->stgs->normalize ? (work->stgs->scale * work->sc_c * work->sc_b) : 1);
480 :
481 201 : nmpr_tau = scs_calc_primal_resid(work, x, s, res->tau, &nmAxs_tau);
482 201 : nmdr_tau = scs_calc_dual_resid(work, y, res->tau, &nmATy_tau);
483 :
484 201 : res->bTy_by_tau =
485 402 : scs_inner_product(y, work->b, m) /
486 201 : (work->stgs->normalize ? (work->stgs->scale * work->sc_c * work->sc_b) : 1);
487 201 : res->cTx_by_tau =
488 402 : scs_inner_product(x, work->c, n) /
489 201 : (work->stgs->normalize ? (work->stgs->scale * work->sc_c * work->sc_b) : 1);
490 :
491 201 : res->res_infeas =
492 201 : res->bTy_by_tau < 0 ? work->nm_b * nmATy_tau / -res->bTy_by_tau : NAN;
493 201 : res->res_unbdd =
494 201 : res->cTx_by_tau < 0 ? work->nm_c * nmAxs_tau / -res->cTx_by_tau : NAN;
495 :
496 201 : bTy = res->bTy_by_tau / res->tau;
497 201 : cTx = res->cTx_by_tau / res->tau;
498 :
499 201 : res->res_pri = nmpr_tau / (1 + work->nm_b) / res->tau;
500 201 : res->res_dual = nmdr_tau / (1 + work->nm_c) / res->tau;
501 201 : res->rel_gap = ABS(cTx + bTy) / (1 + ABS(cTx) + ABS(bTy));
502 : }
503 :
504 30542 : static void scs_calc_residuals_superscs(
505 : ScsWork * RESTRICT work,
506 : struct scs_residuals * RESTRICT residuals,
507 : scs_int iter) {
508 : scs_float * RESTRICT xb;
509 : scs_float * RESTRICT yb;
510 : scs_float * RESTRICT sb;
511 : scs_float cTx;
512 : scs_float bTy;
513 30542 : scs_float * RESTRICT pr = work->pr;
514 30542 : scs_float * RESTRICT dr = work->dr;
515 30542 : scs_int n = work->n;
516 30542 : scs_int m = work->m;
517 : scs_int i;
518 : scs_float norm_D_A_x_plus_s; /* norm of D*(Ax+s), intermediate variable */
519 : scs_float norm_E_Atran_yb; /* norm of E*A'*y, intermediate variable */
520 : scs_float tmp__c_times_x; /* c'x */
521 : scs_float tmp__b_times_yb; /* b'y */
522 30542 : const scs_float temp1 = work->sc_b * work->stgs->scale; /* auxiliary variable #1 */
523 30542 : const scs_float temp2 = work->sc_c * temp1; /* auxiliary variable #2 */
524 30542 : const scs_float temp3 = work->sc_c * work->stgs->scale; /* auxiliary variable #3 */
525 :
526 :
527 : /* checks if the residuals are unchanged by checking iteration */
528 30542 : if (residuals->last_iter == iter) {
529 30542 : return;
530 : }
531 29901 : residuals->last_iter = iter;
532 :
533 29901 : sb = work->s_b;
534 29901 : xb = work->u_b;
535 29901 : yb = &(work->u_b[n]);
536 :
537 29901 : residuals->kap = work->kap_b;
538 29901 : residuals->tau = work->u_b[n + m]; /* it's actually tau_b */
539 29901 : memset(pr, 0, work->m * sizeof (scs_float)); /* pr = 0 */
540 29901 : memset(dr, 0, work->n * sizeof (scs_float)); /* dr = 0 */
541 :
542 29901 : scs_accum_by_a(work->A, work->p, xb, pr); /* pr = A xb */
543 29901 : scs_add_scaled_array(pr, sb, work->m, 1.0); /* pr = A xb + sb */
544 : /* --- compute ||D(Ax + s)|| --- */
545 29901 : norm_D_A_x_plus_s = 0;
546 29901 : if (work->stgs->normalize) {
547 557968 : for (i = 0; i < m; ++i) {
548 557968 : scs_float tmp = work->scal->D[i] * pr[i];
549 557968 : norm_D_A_x_plus_s += tmp * tmp;
550 : }
551 : } else {
552 100 : norm_D_A_x_plus_s = scs_norm_squared(pr, m);
553 : }
554 29901 : norm_D_A_x_plus_s = SQRTF(norm_D_A_x_plus_s);
555 29901 : scs_add_scaled_array(pr, work->b, m, -residuals->tau); /* pr = A xb + sb - b taub */
556 :
557 29901 : scs_accum_by_a_trans(work->A, work->p, yb, dr); /* dr = A' yb */
558 :
559 :
560 : /* --- compute ||E A' yb|| --- */
561 29901 : norm_E_Atran_yb = 0.0;
562 29901 : if (work->stgs->normalize) {
563 290822 : for (i = 0; i < n; ++i) {
564 290822 : scs_float tmp = work->scal->E[i] * dr[i];
565 290822 : norm_E_Atran_yb += tmp * tmp;
566 : }
567 : } else {
568 100 : norm_E_Atran_yb = scs_norm_squared(dr, n);
569 : }
570 29901 : norm_E_Atran_yb = SQRTF(norm_E_Atran_yb);
571 29901 : scs_add_scaled_array(dr, work->c, work->n, residuals->tau); /* dr = A' yb + c taub */
572 :
573 : /*
574 : * bTy_by_tau = b'yb / (scale*sc_c*sc_b)
575 : * cTx_by_tau = c'xb / (scale*sc_c*sc_b)
576 : */
577 29901 : tmp__b_times_yb = scs_inner_product(yb, work->b, m);
578 29901 : residuals->bTy_by_tau = tmp__b_times_yb / (work->stgs->normalize ? (temp2) : 1);
579 29901 : tmp__c_times_x = scs_inner_product(xb, work->c, n);
580 29901 : residuals->cTx_by_tau = tmp__c_times_x / (work->stgs->normalize ? (temp2) : 1);
581 :
582 : /*
583 : * bTy = b'yb / (scale*sc_c*sc_b) / taub
584 : * cTx = c'xb / (scale*sc_c*sc_b) / taub
585 : */
586 29901 : bTy = residuals->bTy_by_tau / residuals->tau;
587 29901 : cTx = residuals->cTx_by_tau / residuals->tau;
588 :
589 : /* PRIMAL RESIDUAL */
590 29901 : if (work->stgs->normalize) {
591 29801 : residuals->res_pri = 0;
592 587769 : for (i = 0; i < m; ++i) {
593 557968 : scs_float tmp = work->scal->D[i] * pr[i];
594 557968 : residuals->res_pri += (tmp * tmp);
595 : }
596 29801 : residuals->res_pri = SQRTF(residuals->res_pri) / residuals->tau;
597 29801 : residuals->res_pri /= ((1 + work->nm_b) * temp1);
598 : } else {
599 100 : residuals->res_pri = scs_norm(pr, m) / residuals->tau;
600 100 : residuals->res_pri /= (1 + work->nm_b);
601 : }
602 :
603 : /* DUAL RESIDUAL */
604 29901 : if (work->stgs->normalize) {
605 29801 : residuals->res_dual = 0;
606 320623 : for (i = 0; i < n; ++i) {
607 290822 : scs_float tmp = work->scal->E[i] * dr[i];
608 290822 : residuals->res_dual += (tmp * tmp);
609 : }
610 29801 : residuals->res_dual = SQRTF(residuals->res_dual) / residuals->tau;
611 29801 : residuals->res_dual /= ((1 + work->nm_c) * temp3);
612 : } else {
613 100 : residuals->res_dual = scs_norm(dr, n) / residuals->tau;
614 100 : residuals->res_dual /= (1 + work->nm_c);
615 : }
616 :
617 : /* UNBOUNDEDNESS */
618 29901 : if (tmp__c_times_x < 0) {
619 28371 : scs_float norm_Ec = 0;
620 28371 : if (work->stgs->normalize) {
621 98759 : for (i = 0; i < n; ++i) {
622 98759 : scs_float tmp = work->scal->E[i] * work->c[i];
623 98759 : norm_Ec += (tmp * tmp);
624 : }
625 : } else {
626 100 : norm_Ec += scs_norm_squared(work->c, n);
627 : }
628 28371 : residuals->res_unbdd = -SQRTF(norm_Ec) * norm_D_A_x_plus_s / tmp__c_times_x;
629 28371 : residuals->res_unbdd /= (work->stgs->normalize ? work->stgs->scale : 1);
630 : } else {
631 1530 : residuals->res_unbdd = NAN; /* not unbounded */
632 : }
633 :
634 :
635 : /* INFEASIBILITY */
636 29901 : if (tmp__b_times_yb < 0) {
637 119 : scs_float norm_Db_squared = 0;
638 119 : if (work->stgs->normalize) {
639 432503 : for (i = 0; i < m; ++i) {
640 432503 : scs_float tmp = work->scal->D[i] * work->b[i];
641 432503 : norm_Db_squared += (tmp * tmp);
642 : }
643 : } else {
644 0 : norm_Db_squared += scs_norm_squared(work->b, m);
645 : }
646 119 : residuals->res_infeas = -SQRTF(norm_Db_squared) * norm_E_Atran_yb / tmp__b_times_yb;
647 119 : residuals->res_infeas /= (work->stgs->normalize ? work->stgs->scale : 1);
648 : } else {
649 29782 : residuals->res_infeas = NAN; /* not infeasible */
650 : }
651 29901 : residuals->rel_gap = ABS(cTx + bTy) / (1 + ABS(cTx) + ABS(bTy));
652 : }
653 :
654 327 : static void scs_cold_start_vars(ScsWork * RESTRICT work) {
655 327 : scs_int l = work->l;
656 327 : memset(work->u, 0, l * sizeof (scs_float));
657 327 : work->u[l - 1] = SQRTF((scs_float) l);
658 327 : if (work->stgs->do_super_scs == 0) {
659 2 : memset(work->v, 0, l * sizeof (scs_float));
660 2 : work->v[l - 1] = SQRTF((scs_float) l);
661 : }
662 327 : }
663 :
664 : /* status < 0 indicates failure */
665 200 : static scs_int scs_project_lin_sys(
666 : ScsWork * RESTRICT work,
667 : scs_int iter) {
668 : /* ut = u + v */
669 200 : scs_int n = work->n, m = work->m, l = n + m + 1, status;
670 200 : memcpy(work->u_t, work->u, l * sizeof (scs_float));
671 200 : scs_add_scaled_array(work->u_t, work->v, l, 1.0);
672 :
673 200 : scs_scale_array(work->u_t, work->stgs->rho_x, n);
674 :
675 200 : scs_add_scaled_array(work->u_t, work->h, l - 1, -work->u_t[l - 1]);
676 200 : scs_add_scaled_array(work->u_t, work->h, l - 1,
677 200 : -scs_inner_product(work->u_t, work->g, l - 1) / (work->gTh + 1));
678 200 : scs_scale_array(&(work->u_t[n]), -1, m);
679 :
680 200 : status = scs_solve_lin_sys(work->A, work->stgs, work->p, work->u_t, work->u, iter);
681 :
682 200 : work->u_t[l - 1] += scs_inner_product(work->u_t, work->h, l - 1);
683 :
684 200 : return status;
685 : }
686 :
687 : /* status < 0 indicates failure */
688 48779 : static scs_int superscs_project_lin_sys(
689 : scs_float * RESTRICT u_t,
690 : scs_float * RESTRICT u,
691 : ScsWork * RESTRICT work,
692 : scs_int iter) {
693 : scs_int status;
694 48779 : const scs_int l = work->l;
695 :
696 : /* x_t = rho_x * x_t */
697 48779 : memcpy(u_t + work->n, u + work->n, (work->m + 1) * sizeof (scs_float));
698 48779 : scs_set_as_scaled_array(u_t, u, work->stgs->rho_x, work->n);
699 :
700 : /* (x_t, y_t) -= tau_t * h */
701 48779 : scs_add_scaled_array(u_t, work->h, l - 1, -u_t[l - 1]);
702 :
703 : /* u_t -= scalar * h */
704 48779 : scs_add_scaled_array(u_t, work->h, l - 1,
705 48779 : -scs_inner_product(u_t, work->g, l - 1) / (work->gTh + 1));
706 :
707 : /* y_t *= (-1) */
708 48779 : scs_scale_array(u_t + work->n, -1, work->m);
709 :
710 : /* call `scs_solve_lin_sys` to update (x_t, y_t) */
711 48779 : status = scs_solve_lin_sys(work->A, work->stgs, work->p, u_t, u, iter);
712 :
713 : /* tau_t += h'*(x_t, y_t) */
714 48779 : u_t[l - 1] += scs_inner_product(u_t, work->h, l - 1);
715 :
716 48779 : return status;
717 : }
718 :
719 : /* LCOV_EXCL_START */
720 : void scs_print_sol(
721 : ScsWork * RESTRICT work,
722 : ScsSolution * RESTRICT sol,
723 : ScsInfo * RESTRICT info) {
724 : scs_int i;
725 : FILE * RESTRICT stream = work->stgs->output_stream;
726 : scs_int print_mode = work->stgs->do_override_streams;
727 : scs_special_print(print_mode, stream, "%s\n", info->status);
728 : if (sol->x != SCS_NULL) {
729 : for (i = 0; i < work->n; ++i) {
730 : scs_special_print(print_mode, stream, "x[%i] = %4f\n", (int) i, sol->x[i]);
731 : }
732 : }
733 : if (sol->y != SCS_NULL) {
734 : for (i = 0; i < work->m; ++i) {
735 : scs_special_print(print_mode, stream, "y[%i] = %4f\n", (int) i, sol->y[i]);
736 : }
737 : }
738 : }
739 :
740 : /* LCOV_EXCL_STOP */
741 :
742 200 : static void scs_update_dual_vars(ScsWork * RESTRICT work) {
743 200 : scs_int n = work->n, l = n + work->m + 1;
744 : /* this does not relax 'x' variable */
745 200 : scs_add_array(work->v, work->u, l);
746 200 : scs_add_scaled_array(work->v, work->u_t, l, -work->stgs->alpha);
747 200 : scs_add_scaled_array(work->v, work->u_prev, l, -1.0 + work->stgs->alpha);
748 200 : }
749 :
750 : /* Calculates the fixed point residual R */
751 : static void scs_calc_FPR(
752 : scs_float * RESTRICT fpr,
753 : scs_float * RESTRICT u_t,
754 : scs_float * RESTRICT u_b,
755 : scs_int l) {
756 79458 : scs_axpy(fpr, u_t, u_b, 1.0, -1.0, l);
757 : }
758 :
759 : /* status < 0 indicates failure */
760 200 : static scs_int scs_project_cones(
761 : ScsWork * RESTRICT work,
762 : const ScsCone * RESTRICT cone,
763 : scs_int iter) {
764 200 : scs_int i, n = work->n, l = n + work->m + 1, status;
765 : /* this does not relax 'x' variable */
766 960666 : for (i = 0; i < n; ++i) {
767 960466 : work->u[i] = work->u_t[i] - work->v[i];
768 : }
769 1921203 : for (i = n; i < l; ++i) {
770 5763609 : work->u[i] = work->stgs->alpha * work->u_t[i] +
771 3842406 : (1 - work->stgs->alpha) * work->u_prev[i] - work->v[i];
772 : }
773 : /* u = [x;y;tau] */
774 200 : status = scs_project_dual_cone(&(work->u[n]), cone, work->coneWork, &(work->u_prev[n]), iter);
775 200 : if (work->u[l - 1] < 0.0)
776 29 : work->u[l - 1] = 0.0;
777 :
778 200 : return status;
779 : }
780 :
781 : /* status < 0 indicates failure */
782 79458 : static scs_int superscs_project_cones(
783 : scs_float * RESTRICT u_b,
784 : scs_float * RESTRICT u_t,
785 : scs_float * RESTRICT u,
786 : ScsWork * RESTRICT work,
787 : const ScsCone * RESTRICT cone,
788 : scs_int iter) {
789 79458 : scs_int n = work->n;
790 79458 : scs_int l = n + work->m + 1;
791 : scs_int status;
792 : /* this does not relax 'x' variable */
793 79458 : scs_axpy(u_b, u_t, u, 2.0, -1.0, l);
794 :
795 : /* u = [x;y;tau] */
796 79458 : status = scs_project_dual_cone(&(u_b[n]), cone, work->coneWork, &(work->u_prev[n]), iter);
797 79458 : if (u_b[l - 1] < 0.0) {
798 3565 : u_b[l - 1] = 0.0;
799 : }
800 79458 : return status;
801 : }
802 :
803 : /* LCOV_EXCL_START */
804 : static scs_int scs_indeterminate(
805 : ScsWork * RESTRICT w,
806 : ScsSolution * RESTRICT sol,
807 : ScsInfo * RESTRICT info) {
808 : strncpy(info->status, "Indeterminate", SCS_INFO_STATUS_MSG_LENGTH);
809 : scs_scale_array(sol->x, NAN, w->n);
810 : scs_scale_array(sol->y, NAN, w->m);
811 : scs_scale_array(sol->s, NAN, w->m);
812 : return SCS_INDETERMINATE;
813 : }
814 :
815 : /* LCOV_EXCL_STOP */
816 :
817 323 : static scs_int scs_solved(
818 : ScsWork * RESTRICT work,
819 : ScsSolution * RESTRICT sol,
820 : ScsInfo * RESTRICT info,
821 : scs_float tau) {
822 323 : scs_scale_array(sol->x, 1.0 / tau, work->n);
823 323 : scs_scale_array(sol->y, 1.0 / tau, work->m);
824 323 : scs_scale_array(sol->s, 1.0 / tau, work->m);
825 323 : if (info->statusVal == 0) {
826 11 : strncpy(info->status, "Solved/Inaccurate", SCS_INFO_STATUS_MSG_LENGTH);
827 : return SCS_SOLVED_INACCURATE;
828 : }
829 312 : strncpy(info->status, "Solved", SCS_INFO_STATUS_MSG_LENGTH);
830 : return SCS_SOLVED;
831 : }
832 :
833 1 : static scs_int scs_infeasible(
834 : ScsWork * RESTRICT work,
835 : ScsSolution * RESTRICT sol,
836 : ScsInfo * RESTRICT info,
837 : scs_float bTy) {
838 1 : scs_scale_array(sol->y, -1 / bTy, work->m);
839 1 : scs_scale_array(sol->x, NAN, work->n);
840 1 : scs_scale_array(sol->s, NAN, work->m);
841 1 : if (info->statusVal == 0) {
842 0 : strncpy(info->status, "Infeasible/Inaccurate", SCS_INFO_STATUS_MSG_LENGTH);
843 : return SCS_INFEASIBLE_INACCURATE;
844 : }
845 1 : strncpy(info->status, "Infeasible", SCS_INFO_STATUS_MSG_LENGTH);
846 : return SCS_INFEASIBLE;
847 : }
848 :
849 5 : static scs_int scs_unbounded(
850 : ScsWork * RESTRICT work,
851 : ScsSolution * RESTRICT sol,
852 : ScsInfo * RESTRICT info,
853 : scs_float cTx) {
854 5 : scs_scale_array(sol->x, -1 / cTx, work->n);
855 5 : scs_scale_array(sol->s, -1 / cTx, work->m);
856 5 : scs_scale_array(sol->y, NAN, work->m);
857 5 : if (info->statusVal == 0) {
858 3 : strncpy(info->status, "Unbounded/Inaccurate", SCS_INFO_STATUS_MSG_LENGTH);
859 : return SCS_UNBOUNDED_INACCURATE;
860 : }
861 2 : strncpy(info->status, "Unbounded", SCS_INFO_STATUS_MSG_LENGTH);
862 : return SCS_UNBOUNDED;
863 : }
864 :
865 329 : static void scs_set_y(ScsWork * RESTRICT work, ScsSolution * RESTRICT sol) {
866 329 : if (sol->y == SCS_NULL) {
867 18 : sol->y = scs_malloc(sizeof (scs_float) * work->m);
868 : }
869 329 : if (work->stgs->do_super_scs == 0) {
870 3 : memcpy(sol->y, &(work->u[work->n]), work->m * sizeof (scs_float));
871 : } else {
872 326 : memcpy(sol->y, &(work->u_b[work->n]), work->m * sizeof (scs_float));
873 : }
874 329 : }
875 :
876 329 : static void scs_set_s(ScsWork * RESTRICT work, ScsSolution * RESTRICT sol) {
877 329 : if (sol->s == SCS_NULL) {
878 18 : sol->s = scs_malloc(sizeof (scs_float) * work->m);
879 : }
880 329 : if (work->stgs->do_super_scs == 0) {
881 3 : memcpy(sol->s, &(work->v[work->n]), work->m * sizeof (scs_float));
882 : } else {
883 326 : memcpy(sol->s, work->s_b, work->m * sizeof (scs_float));
884 : }
885 329 : }
886 :
887 329 : static void scs_set_x(ScsWork * RESTRICT work, ScsSolution * RESTRICT sol) {
888 329 : if (sol->x == SCS_NULL)
889 18 : sol->x = scs_malloc(sizeof (scs_float) * work->n);
890 329 : if (work->stgs->do_super_scs == 0) {
891 3 : memcpy(sol->x, work->u, work->n * sizeof (scs_float));
892 : } else {
893 326 : memcpy(sol->x, work->u_b, work->n * sizeof (scs_float));
894 : }
895 329 : }
896 :
897 644 : scs_int scs_is_solved_status(scs_int status) {
898 644 : return status == SCS_SOLVED || status == SCS_SOLVED_INACCURATE;
899 : }
900 :
901 6 : scs_int scs_is_infeasible_status(scs_int status) {
902 6 : return status == SCS_INFEASIBLE || status == SCS_INFEASIBLE_INACCURATE;
903 : }
904 :
905 8 : scs_int scs_is_unbounded_status(scs_int status) {
906 8 : return status == SCS_UNBOUNDED || status == SCS_UNBOUNDED_INACCURATE;
907 : }
908 :
909 329 : static void scs_get_info(
910 : ScsWork * RESTRICT work,
911 : ScsSolution * RESTRICT sol,
912 : ScsInfo * RESTRICT info,
913 : struct scs_residuals * RESTRICT residuals,
914 : scs_int iter) {
915 329 : info->iter = iter;
916 329 : info->resInfeas = residuals->res_infeas;
917 329 : info->resUnbdd = residuals->res_unbdd;
918 329 : if (scs_is_solved_status(info->statusVal)) {
919 323 : info->relGap = residuals->rel_gap;
920 323 : info->resPri = residuals->res_pri;
921 323 : info->resDual = residuals->res_dual;
922 323 : info->pobj = residuals->cTx_by_tau / residuals->tau;
923 323 : info->dobj = -residuals->bTy_by_tau / residuals->tau;
924 6 : } else if (scs_is_unbounded_status(info->statusVal)) {
925 5 : info->relGap = NAN;
926 5 : info->resPri = NAN;
927 5 : info->resDual = NAN;
928 5 : info->pobj = -INFINITY;
929 5 : info->dobj = -INFINITY;
930 1 : } else if (scs_is_infeasible_status(info->statusVal)) {
931 1 : info->relGap = NAN;
932 1 : info->resPri = NAN;
933 1 : info->resDual = NAN;
934 1 : info->pobj = INFINITY;
935 1 : info->dobj = INFINITY;
936 : }
937 329 : }
938 :
939 : /* sets solutions, re-scales by inner prods if infeasible or unbounded */
940 329 : static void scs_get_solution(
941 328 : ScsWork * RESTRICT work,
942 : ScsSolution * RESTRICT sol,
943 : ScsInfo * RESTRICT info,
944 : struct scs_residuals * RESTRICT r,
945 : scs_int iter) {
946 329 : scs_int l = work->l;
947 329 : if (work->stgs->do_super_scs == 0) {
948 3 : scs_calc_residuals(work, r, iter);
949 : } else {
950 326 : scs_calc_residuals_superscs(work, r, iter);
951 652 : r->kap = ABS(work->kap_b) /
952 326 : (work->stgs->normalize ? (work->stgs->scale * work->sc_c * work->sc_b) : 1.0);
953 : }
954 329 : scs_set_x(work, sol);
955 329 : scs_set_y(work, sol);
956 329 : scs_set_s(work, sol);
957 329 : if (info->statusVal == SCS_UNFINISHED) {
958 : /* not yet converged, take best guess */
959 14 : if (r->tau > SCS_INDETERMINATE_TOL && r->tau > r->kap) {
960 22 : info->statusVal = scs_solved(work, sol, info, r->tau);
961 6 : } else if (scs_norm(work->u, l) <
962 3 : SCS_INDETERMINATE_TOL * SQRTF((scs_float) l)) {
963 0 : info->statusVal = scs_indeterminate(work, sol, info);
964 3 : } else if (r->bTy_by_tau < r->cTx_by_tau) {
965 0 : info->statusVal = scs_infeasible(work, sol, info, r->bTy_by_tau);
966 : } else {
967 6 : info->statusVal = scs_unbounded(work, sol, info, r->cTx_by_tau);
968 : }
969 315 : } else if (scs_is_solved_status(info->statusVal)) {
970 624 : info->statusVal = scs_solved(work, sol, info, r->tau);
971 3 : } else if (scs_is_infeasible_status(info->statusVal)) {
972 1 : info->statusVal = scs_infeasible(work, sol, info, r->bTy_by_tau);
973 : } else {
974 4 : info->statusVal = scs_unbounded(work, sol, info, r->cTx_by_tau);
975 : }
976 329 : if (work->stgs->normalize) {
977 326 : scs_unnormalize_sol(work, sol);
978 : }
979 329 : scs_get_info(work, sol, info, r, iter);
980 329 : }
981 :
982 : /* LCOV_EXCL_START */
983 : static void scs_print_summary(
984 : ScsWork * RESTRICT work,
985 : scs_int i,
986 : struct scs_residuals * RESTRICT residuals,
987 : ScsTimer * RESTRICT solveTimer) {
988 : FILE * RESTRICT stream = work->stgs->output_stream;
989 : scs_int print_mode = work->stgs->do_override_streams;
990 : scs_special_print(print_mode, stream, "%*i|", (int) strlen(SCS_HEADER[0]), (int) i);
991 : scs_special_print(print_mode, stream, "%*.2e ", (int) scs_hspace, residuals->res_pri);
992 : scs_special_print(print_mode, stream, "%*.2e ", (int) scs_hspace, residuals->res_dual);
993 : scs_special_print(print_mode, stream, "%*.2e ", (int) scs_hspace, residuals->rel_gap);
994 : scs_special_print(print_mode, stream, "%*.2e ", (int) scs_hspace, residuals->cTx_by_tau / residuals->tau);
995 : scs_special_print(print_mode, stream, "%*.2e ", (int) scs_hspace, -residuals->bTy_by_tau / residuals->tau);
996 : scs_special_print(print_mode, stream, "%*.2e ", (int) scs_hspace, residuals->kap / residuals->tau);
997 : if (work->stgs->do_super_scs) {
998 : scs_special_print(print_mode, stream, "%*.2e ", (int) scs_hspace, work->nrmR_con);
999 : }
1000 : scs_special_print(print_mode, stream, "%*.2e ", (int) scs_hspace, scs_toc_quiet(solveTimer) / 1e3);
1001 : scs_special_print(print_mode, stream, "\n");
1002 :
1003 : #ifdef MATLAB_MEX_FILE
1004 : mexEvalString("drawnow;");
1005 : #endif
1006 : }
1007 :
1008 : static void scs_print_header(
1009 : ScsWork * RESTRICT work,
1010 : const ScsCone * RESTRICT cone) {
1011 : scs_int i;
1012 : FILE * RESTRICT stream = work->stgs->output_stream;
1013 : scs_int print_mode = work->stgs->do_override_streams;
1014 : if (work->stgs->warm_start)
1015 : scs_special_print(print_mode, stream, "SCS using variable warm-starting\n");
1016 : for (i = 0; i < scs_header_line_length; ++i) {
1017 : scs_special_print(print_mode, stream, "-");
1018 : }
1019 : scs_special_print(print_mode, stream, "\n");
1020 : for (i = 0; i < scs_header_length - 2; ++i) {
1021 : scs_special_print(print_mode, stream, "%s|", SCS_HEADER[i]);
1022 : }
1023 : if (work->stgs->do_super_scs) {
1024 : scs_special_print(print_mode, stream, "%s|", SCS_HEADER[scs_header_length - 2]);
1025 : }
1026 : scs_special_print(print_mode, stream, "%s\n", SCS_HEADER[scs_header_length - 1]);
1027 : for (i = 0; i < scs_header_line_length; ++i) {
1028 : scs_special_print(print_mode, stream, "-");
1029 : }
1030 : scs_special_print(print_mode, stream, "\n");
1031 : #ifdef MATLAB_MEX_FILE
1032 : mexEvalString("drawnow;");
1033 : #endif
1034 : }
1035 :
1036 : /* LCOV_EXCL_STOP */
1037 :
1038 2 : scs_float scs_get_dual_cone_dist(
1039 : const scs_float * RESTRICT y,
1040 : const ScsCone * RESTRICT k,
1041 : ScsConeWork * RESTRICT cone_workspace,
1042 : scs_int m) {
1043 : scs_float dist;
1044 2 : scs_float * RESTRICT t = scs_malloc(sizeof (scs_float) * m);
1045 2 : memcpy(t, y, m * sizeof (scs_float));
1046 2 : scs_project_dual_cone(t, k, cone_workspace, SCS_NULL, -1);
1047 2 : dist = scs_norm_infinity_difference(t, y, m);
1048 2 : scs_free(t);
1049 2 : return dist;
1050 : }
1051 :
1052 : /* via moreau */
1053 2 : scs_float scs_get_primal_cone_dist(
1054 : const scs_float * RESTRICT s,
1055 : const ScsCone * RESTRICT cone,
1056 : ScsConeWork * RESTRICT cone_workspace,
1057 : scs_int m) {
1058 : scs_float dist;
1059 2 : scs_float * RESTRICT t = scs_malloc(sizeof (scs_float) * m);
1060 2 : memcpy(t, s, m * sizeof (scs_float));
1061 2 : scs_scale_array(t, -1.0, m);
1062 2 : scs_project_dual_cone(t, cone, cone_workspace, SCS_NULL, -1);
1063 2 : dist = scs_norm_infinity(t, m); /* ||s - Pi_c(s)|| = ||Pi_c*(-s)|| */
1064 2 : scs_free(t);
1065 2 : return dist;
1066 : }
1067 :
1068 6 : void scs_millis_to_time(scs_float time,
1069 : scs_int * hours,
1070 : scs_int * minutes,
1071 : scs_int * secs,
1072 : scs_float * sec_rest) {
1073 : scs_float integral;
1074 6 : scs_float time_seconds = time / 1e3;
1075 6 : *sec_rest = (scs_float) modf((double) time_seconds, &integral);
1076 6 : *secs = (scs_int) time_seconds;
1077 6 : *minutes = *secs / 60;
1078 6 : *secs %= 60;
1079 6 : if (*minutes >= 60) {
1080 2 : *hours = *minutes / 60;
1081 2 : *minutes %= 60;
1082 : } else {
1083 4 : *hours = 0;
1084 : }
1085 6 : }
1086 :
1087 : /* LCOV_EXCL_START */
1088 : static void scs_print_footer(
1089 : const ScsData * RESTRICT data,
1090 : const ScsCone * RESTRICT cone,
1091 : ScsSolution * RESTRICT sol,
1092 : ScsWork * RESTRICT work,
1093 : ScsInfo * RESTRICT info) {
1094 : scs_int i;
1095 : scs_int hours, minutes, seconds;
1096 : scs_float millis;
1097 :
1098 : char *linSysStr = scs_get_linsys_summary(work->p, info);
1099 : char *coneStr = scs_get_cone_summary(info, work->coneWork);
1100 : FILE * stream = work->stgs->output_stream;
1101 : scs_int print_mode = work->stgs->do_override_streams;
1102 : for (i = 0; i < scs_header_line_length; ++i) {
1103 : scs_special_print(print_mode, stream, "-");
1104 : }
1105 : scs_special_print(print_mode, stream, "\nStatus: %s\n", info->status);
1106 : if (info->iter == work->stgs->max_iters) {
1107 : scs_special_print(print_mode, stream, "Hit max_iters, solution may be inaccurate\n");
1108 : }
1109 :
1110 : scs_millis_to_time(info->solveTime, &hours, &minutes, &seconds, &millis);
1111 : scs_special_print(print_mode, stream, "Timing: Solve time: %02d:%02d:%02d.%d\n",
1112 : hours, minutes, seconds, (scs_int) round(1e4 * millis) / 10);
1113 :
1114 : if (linSysStr) {
1115 : scs_special_print(print_mode, stream, "%s", linSysStr);
1116 : scs_free(linSysStr);
1117 : }
1118 :
1119 : if (coneStr) {
1120 : scs_special_print(print_mode, stream, "%s", coneStr);
1121 : scs_free(coneStr);
1122 : }
1123 :
1124 : for (i = 0; i < scs_header_line_length; ++i) {
1125 : scs_special_print(print_mode, stream, "-");
1126 : }
1127 : scs_special_print(print_mode, stream, "\n");
1128 :
1129 : if (scs_is_infeasible_status(info->statusVal)) {
1130 : scs_special_print(print_mode, stream, "Certificate of primal infeasibility:\n");
1131 : scs_special_print(print_mode, stream, "dist(y, K*) = %.4e\n",
1132 : scs_get_dual_cone_dist(sol->y, cone, work->coneWork, data->m));
1133 : scs_special_print(print_mode, stream, "|A'y|_2 * |b|_2 = %.4e\n", info->resInfeas);
1134 : scs_special_print(print_mode, stream, "b'y = %.4f\n", scs_inner_product(data->b, sol->y, data->m));
1135 : } else if (scs_is_unbounded_status(info->statusVal)) {
1136 : scs_special_print(print_mode, stream, "Certificate of dual infeasibility:\n");
1137 : scs_special_print(print_mode, stream, "dist(s, K) = %.4e\n",
1138 : scs_get_primal_cone_dist(sol->s, cone, work->coneWork, data->m));
1139 : scs_special_print(print_mode, stream, "|Ax + s|_2 * |c|_2 = %.4e\n", info->resUnbdd);
1140 : scs_special_print(print_mode, stream, "c'x = %.4f\n", scs_inner_product(data->c, sol->x, data->n));
1141 : } else {
1142 : scs_special_print(print_mode, stream, "Error metrics:\n");
1143 : scs_special_print(print_mode, stream, "dist(s, K) = %.4e, dist(y, K*) = %.4e, s'y/|s||y| = %.4e\n",
1144 : scs_get_primal_cone_dist(sol->s, cone, work->coneWork, data->m),
1145 : scs_get_dual_cone_dist(sol->y, cone, work->coneWork, data->m),
1146 : scs_inner_product(sol->s, sol->y, data->m) / scs_norm(sol->s, data->m) /
1147 : scs_norm(sol->y, data->m));
1148 : scs_special_print(print_mode, stream, "|Ax + s - b|_2 / (1 + |b|_2) = %.4e\n", info->resPri);
1149 : scs_special_print(print_mode, stream, "|A'y + c|_2 / (1 + |c|_2) = %.4e\n", info->resDual);
1150 : scs_special_print(print_mode, stream, "|c'x + b'y| / (1 + |c'x| + |b'y|) = %.4e\n", info->relGap);
1151 : for (i = 0; i < scs_header_line_length; ++i) {
1152 : scs_special_print(print_mode, stream, "-");
1153 : }
1154 : scs_special_print(print_mode, stream, "\n");
1155 : scs_special_print(print_mode, stream, "c'x = %.4f, -b'y = %.4f\n", info->pobj, info->dobj);
1156 : }
1157 : for (i = 0; i < scs_header_line_length; ++i) {
1158 : scs_special_print(print_mode, stream, "=");
1159 : }
1160 : scs_special_print(print_mode, stream, "\n");
1161 : #ifdef MATLAB_MEX_FILE
1162 : mexEvalString("drawnow;");
1163 : #endif
1164 : }
1165 :
1166 : /* LCOV_EXCL_STOP */
1167 :
1168 : static scs_int scs_has_converged(
1169 : ScsWork * RESTRICT work,
1170 : struct scs_residuals * RESTRICT residuals,
1171 : scs_int iter) {
1172 30088 : scs_float eps = work->stgs->eps;
1173 30088 : if (residuals->res_pri < eps && residuals->res_dual < eps && residuals->rel_gap < eps) {
1174 : return SCS_SOLVED;
1175 : }
1176 29776 : if (residuals->res_unbdd < eps) {
1177 : return SCS_UNBOUNDED;
1178 : }
1179 29774 : if (residuals->res_infeas < eps) {
1180 : return SCS_INFEASIBLE;
1181 : }
1182 : return 0;
1183 : }
1184 :
1185 339 : static scs_int scs_validate_superscs_settings(const ScsData *data) {
1186 339 : ScsSettings * stgs = data->stgs;
1187 339 : scs_int print_mode = stgs->do_override_streams;
1188 :
1189 339 : if (stgs->do_super_scs == 0) return 0;
1190 :
1191 336 : if (stgs->thetabar < 0 || stgs->thetabar > 1) {
1192 : /* LCOV_EXCL_START */
1193 : scs_special_print(print_mode, stderr, "Parameters `thetabar` must "
1194 : "be a scalar between 0 and 1 (thetabar=%g)\n", stgs->thetabar);
1195 : return SCS_FAILED;
1196 : /* LCOV_EXCL_STOP */
1197 : }
1198 672 : if ((stgs->direction == restarted_broyden
1199 336 : || stgs->direction == anderson_acceleration)
1200 331 : && stgs->memory <= 1) {
1201 : /* LCOV_EXCL_START */
1202 : scs_special_print(print_mode, stderr, "Quasi-Newton memory length "
1203 : "(mem=%ld) is too low; choose an integer at least equal to 2.\n",
1204 : (long) stgs->memory);
1205 : return SCS_FAILED;
1206 : /* LCOV_EXCL_STOP */
1207 : }
1208 336 : if (stgs->direction == anderson_acceleration && stgs->memory > data->m + data->n + 1) {
1209 0 : scs_special_print(print_mode, stderr, "Quasi-Newton memory length (mem=%ld) "
1210 : "is too high for Anderson's method (l=%d).\n",
1211 : (long) stgs->memory, data->m + data->n + 1);
1212 : return SCS_FAILED;
1213 : }
1214 336 : if (stgs->beta >= 1 || stgs->beta <= 0) {
1215 : /* LCOV_EXCL_START */
1216 : scs_special_print(print_mode, stderr, "Stepsize reduction factor (beta=%g) out of bounds.\n", stgs->beta);
1217 : return SCS_FAILED;
1218 : /* LCOV_EXCL_STOP */
1219 : }
1220 334 : if (stgs->ls < 0) {
1221 : /* LCOV_EXCL_START */
1222 : scs_special_print(print_mode, stderr, "Illegal maximum number of line search iterations (ls=%ld).\n", (long) stgs->ls);
1223 : return SCS_FAILED;
1224 : /* LCOV_EXCL_STOP */
1225 : }
1226 333 : if (stgs->ls >= 40) {
1227 : /* LCOV_EXCL_START */
1228 : scs_special_print(print_mode, stderr, "WARNING! The value ls=%ld is too high. The maximum allowed "
1229 : "number of line search iteration is 40. We recommend a value about 10.\n", (long) stgs->ls);
1230 : return SCS_FAILED;
1231 : /* LCOV_EXCL_STOP */
1232 : }
1233 332 : if (stgs->ls > 10) {
1234 : /* LCOV_EXCL_START */
1235 : scs_special_print(print_mode, stderr, "WARNING! The value ls=%ld is too high. We highly recommend"
1236 : "the maximum number of line search iterations to be at most 10.\n", (long) stgs->ls);
1237 : /* LCOV_EXCL_STOP */
1238 : }
1239 332 : if (stgs->sigma < 0) {
1240 : /* LCOV_EXCL_START */
1241 : scs_special_print(print_mode, stderr, "Parameter sigma of the line search (sigma=%g) cannot be negative.\n", stgs->sigma);
1242 : return SCS_FAILED;
1243 : /* LCOV_EXCL_STOP */
1244 : }
1245 331 : if (stgs->c_bl < 0 || stgs->c_bl >= 1) {
1246 : /* LCOV_EXCL_START */
1247 : scs_special_print(print_mode, stderr, "Parameter (c_0=%g) for blind updates out of bounds.\n", stgs->c_bl);
1248 : return SCS_FAILED;
1249 : /* LCOV_EXCL_STOP */
1250 : }
1251 330 : if (stgs->c1 < 0 || stgs->c1 >= 1) {
1252 : /* LCOV_EXCL_START */
1253 : scs_special_print(print_mode, stderr, "Parameter (c1=%g) for step K1 out of bounds.\n", stgs->c1);
1254 : return SCS_FAILED;
1255 : /* LCOV_EXCL_STOP */
1256 : }
1257 329 : if (stgs->sse < 0 || stgs->sse >= 1) {
1258 : /* LCOV_EXCL_START */
1259 : scs_special_print(print_mode, stderr, "Parameter (sse=%g) for step K1 out of bounds.\n", stgs->sse);
1260 : return SCS_FAILED;
1261 : /* LCOV_EXCL_STOP */
1262 : }
1263 329 : if (stgs->k0 != 0 && stgs->k0 != 1) {
1264 : /* LCOV_EXCL_START */
1265 : scs_special_print(print_mode, stderr, "Parameter (k0=%d) can be eiter 0 (k0: off) or 1 (k0: on).\n", (int) stgs->k0);
1266 : return SCS_FAILED;
1267 : /* LCOV_EXCL_STOP */
1268 : }
1269 328 : if (stgs->k1 != 0 && stgs->k1 != 1) {
1270 : /* LCOV_EXCL_START */
1271 : scs_special_print(print_mode, stderr, "Parameter (k1=%d) can be eiter 0 (k1: off) or 1 (k1: on).\n", (int) stgs->k1);
1272 : return SCS_FAILED;
1273 : /* LCOV_EXCL_STOP */
1274 : }
1275 327 : if (stgs->k2 != 0 && stgs->k2 != 1) {
1276 : /* LCOV_EXCL_START */
1277 : scs_special_print(print_mode, stderr, "Parameter (k2=%d) can be eiter 0 (k2: off) or 1 (k2: on).\n", (int) stgs->k2);
1278 : return SCS_FAILED;
1279 : /* LCOV_EXCL_STOP */
1280 : }
1281 652 : if (stgs->direction != restarted_broyden
1282 326 : && stgs->direction != fixed_point_residual
1283 3 : && stgs->direction != anderson_acceleration
1284 0 : && stgs->direction != full_broyden) {
1285 : /* LCOV_EXCL_START */
1286 : scs_special_print(print_mode, stderr, "Invalid direction (%ld).\n",
1287 : (long) stgs->direction);
1288 : return SCS_FAILED;
1289 : /* LCOV_EXCL_STOP */
1290 : }
1291 : return 0;
1292 : }
1293 :
1294 341 : static scs_int scs_validate_general_settings(const ScsData *data) {
1295 341 : ScsSettings * stgs = data->stgs;
1296 341 : scs_int print_mode = stgs->do_override_streams;
1297 341 : if (stgs->max_iters <= 0) {
1298 : /* LCOV_EXCL_START */
1299 : scs_special_print(print_mode, stderr, "max_iters must be positive (max_iters=%ld)\n", (long) stgs->max_iters);
1300 : return SCS_FAILED;
1301 : /* LCOV_EXCL_STOP */
1302 : }
1303 341 : if (stgs->eps <= 0) {
1304 : /* LCOV_EXCL_START */
1305 : scs_special_print(print_mode, stderr, "eps tolerance must be positive (eps=%g)\n", stgs->eps);
1306 : return SCS_FAILED;
1307 : /* LCOV_EXCL_STOP */
1308 : }
1309 341 : if (stgs->alpha <= 0 || stgs->alpha >= 2) {
1310 : /* LCOV_EXCL_START */
1311 : scs_special_print(print_mode, stderr, "alpha must be in (0,2) (alpha=%g)\n", stgs->alpha);
1312 : return SCS_FAILED;
1313 : /* LCOV_EXCL_STOP */
1314 : }
1315 339 : if (stgs->rho_x <= 0) {
1316 : /* LCOV_EXCL_START */
1317 : scs_special_print(print_mode, stderr, "rho_x must be positive (1e-3 works well) (rho_x=%g).\n", stgs->rho_x);
1318 : return SCS_FAILED;
1319 : /* LCOV_EXCL_STOP */
1320 : }
1321 339 : if (stgs->scale <= 0) {
1322 : /* LCOV_EXCL_START */
1323 : scs_special_print(print_mode, stderr, "Parameter `scale` must be positive (1 works well).\n");
1324 : return SCS_FAILED;
1325 : /* LCOV_EXCL_STOP */
1326 : }
1327 339 : if (stgs->do_super_scs != 0 && stgs->do_super_scs != 1) {
1328 : /* LCOV_EXCL_START */
1329 : scs_special_print(print_mode, stderr, "do_super_scs (=%d) can be either 0 or 1.\n", (int) stgs->do_super_scs);
1330 : return SCS_FAILED;
1331 : /* LCOV_EXCL_STOP */
1332 : }
1333 :
1334 : return 0;
1335 : }
1336 :
1337 341 : static scs_int scs_validate(
1338 1358 : const ScsData * RESTRICT data,
1339 : const ScsCone * RESTRICT cone) {
1340 341 : ScsSettings *stgs = data->stgs;
1341 341 : scs_int print_mode = stgs->do_override_streams;
1342 341 : if (data->m <= 0 || data->n <= 0) {
1343 : /* LCOV_EXCL_START */
1344 : scs_special_print(print_mode, stderr, "m and n must both be greater than 0; m = %li, n = %li\n",
1345 : (long) data->m, (long) data->n);
1346 : return -1;
1347 : /* LCOV_EXCL_STOP */
1348 : }
1349 341 : if (data->m < data->n) {
1350 : /* LCOV_EXCL_START */
1351 : scs_special_print(print_mode, stderr, "WARN: m less than n, problem likely degenerate\n");
1352 : /* LCOV_EXCL_STOP */
1353 : }
1354 341 : if (scs_validate_linsys(data->A) < 0) {
1355 : /* LCOV_EXCL_START */
1356 : scs_special_print(print_mode, stderr, "invalid linear system input data\n");
1357 : return SCS_FAILED;
1358 : /* LCOV_EXCL_STOP */
1359 : }
1360 341 : if (scs_validate_cones(data, cone) < 0) {
1361 : /* LCOV_EXCL_START */
1362 : scs_special_print(print_mode, stderr, "cone validation error\n");
1363 : return SCS_FAILED;
1364 : /* LCOV_EXCL_STOP */
1365 : }
1366 :
1367 341 : if (scs_validate_general_settings(data) != 0
1368 1017 : || scs_validate_superscs_settings(data) != 0)
1369 : return SCS_FAILED;
1370 :
1371 : return 0;
1372 : }
1373 :
1374 329 : static ScsWork * scs_init_work(
1375 : const ScsData * RESTRICT data,
1376 : const ScsCone * RESTRICT cone) {
1377 329 : ScsWork *work = scs_calloc(1, sizeof (*work));
1378 329 : const scs_int l = data->n + data->m + 1;
1379 329 : scs_int print_mode = data->stgs->do_override_streams;
1380 329 : if (data->stgs->verbose) {
1381 2 : scs_print_init_header(data, cone);
1382 : }
1383 329 : if (work == SCS_NULL) {
1384 0 : scs_special_print(print_mode, stderr, "ERROR: allocating work failure\n");
1385 0 : return SCS_NULL;
1386 : }
1387 :
1388 : /* get settings and dims from data struct */
1389 329 : work->stgs = data->stgs;
1390 329 : work->m = data->m;
1391 329 : work->n = data->n;
1392 329 : work->l = l; /* total dimension */
1393 :
1394 : /* -------------------------------------
1395 : * Workspace allocation:
1396 : *
1397 : * After every call to scs_malloc or scs_calloc
1398 : * we check whether the allocation has been
1399 : * successful.
1400 : * ------------------------------------- */
1401 329 : work->u = scs_calloc(l, sizeof (scs_float));
1402 329 : if (work->u == SCS_NULL) {
1403 : /* LCOV_EXCL_START */
1404 : scs_special_print(print_mode, stderr, "ERROR: `u` memory allocation failure\n");
1405 : return SCS_NULL;
1406 : /* LCOV_EXCL_STOP */
1407 : }
1408 329 : work->u_b = scs_calloc(l, sizeof (scs_float));
1409 329 : if (work->u_b == SCS_NULL) {
1410 : /* LCOV_EXCL_START */
1411 : scs_special_print(print_mode, stderr, "ERROR: `u_b` memory allocation failure\n");
1412 : return SCS_NULL;
1413 : /* LCOV_EXCL_STOP */
1414 : }
1415 329 : if (work->stgs->do_super_scs == 0) {
1416 3 : work->v = scs_calloc(l, sizeof (scs_float));
1417 3 : if (work->v == SCS_NULL) {
1418 : /* LCOV_EXCL_START */
1419 : scs_special_print(print_mode, stderr, "ERROR: `v` memory allocation failure\n");
1420 : return SCS_NULL;
1421 : /* LCOV_EXCL_STOP */
1422 : }
1423 : }
1424 329 : work->u_t = scs_malloc(l * sizeof (scs_float));
1425 329 : if (work->u_t == SCS_NULL) {
1426 : /* LCOV_EXCL_START */
1427 : scs_special_print(print_mode, stderr, "ERROR: `u_t` memory allocation failure\n");
1428 : return SCS_NULL;
1429 : /* LCOV_EXCL_STOP */
1430 : }
1431 329 : work->u_prev = scs_malloc(l * sizeof (scs_float));
1432 329 : if (work->u_prev == SCS_NULL) {
1433 : /* LCOV_EXCL_START */
1434 : scs_special_print(print_mode, stderr, "ERROR: `u_prev` memory allocation failure\n");
1435 : return SCS_NULL;
1436 : /* LCOV_EXCL_STOP */
1437 : }
1438 329 : work->h = scs_malloc((l - 1) * sizeof (scs_float));
1439 329 : if (work->h == SCS_NULL) {
1440 : /* LCOV_EXCL_START */
1441 : scs_special_print(print_mode, stderr, "ERROR: `h` memory allocation failure\n");
1442 : return SCS_NULL;
1443 : /* LCOV_EXCL_STOP */
1444 : }
1445 329 : work->g = scs_malloc((l - 1) * sizeof (scs_float));
1446 329 : if (work->g == SCS_NULL) {
1447 : /* LCOV_EXCL_START */
1448 : scs_special_print(print_mode, stderr, "ERROR: `g` memory allocation failure\n");
1449 : return SCS_NULL;
1450 : /* LCOV_EXCL_STOP */
1451 : }
1452 329 : work->pr = scs_malloc(data->m * sizeof (scs_float));
1453 329 : if (work->pr == SCS_NULL) {
1454 : /* LCOV_EXCL_START */
1455 : scs_special_print(print_mode, stderr, "ERROR: `pr` memory allocation failure\n");
1456 : return SCS_NULL;
1457 : /* LCOV_EXCL_STOP */
1458 : }
1459 329 : work->dr = scs_malloc(data->n * sizeof (scs_float));
1460 329 : if (work->dr == SCS_NULL) {
1461 : /* LCOV_EXCL_START */
1462 : scs_special_print(print_mode, stderr, "ERROR: `dr` memory allocation failure\n");
1463 : return SCS_NULL;
1464 : /* LCOV_EXCL_STOP */
1465 : }
1466 329 : work->b = scs_malloc(data->m * sizeof (scs_float));
1467 329 : if (work->b == SCS_NULL) {
1468 : /* LCOV_EXCL_START */
1469 : scs_special_print(print_mode, stderr, "ERROR: `b` memory allocation failure\n");
1470 : return SCS_NULL;
1471 : /* LCOV_EXCL_STOP */
1472 : }
1473 329 : work->c = scs_malloc(data->n * sizeof (scs_float));
1474 329 : if (work->c == SCS_NULL) {
1475 : /* LCOV_EXCL_START */
1476 : scs_special_print(print_mode, stderr, "ERROR: `c` memory allocation failure\n");
1477 : return SCS_NULL;
1478 : /* LCOV_EXCL_STOP */
1479 : }
1480 :
1481 :
1482 :
1483 329 : if (work->stgs->do_super_scs == 1) {
1484 : /* -------------------------------------
1485 : * Additional memory needs to be allocated
1486 : * in SuperSCS
1487 : * ------------------------------------- */
1488 326 : work->R = scs_calloc(l, sizeof (scs_float));
1489 326 : if (work->R == SCS_NULL) {
1490 : /* LCOV_EXCL_START */
1491 : scs_special_print(print_mode, stderr, "ERROR: `R` memory allocation failure\n");
1492 : return SCS_NULL;
1493 : /* LCOV_EXCL_STOP */
1494 : }
1495 326 : work->R_prev = scs_calloc(l, sizeof (scs_float));
1496 326 : if (work->R_prev == SCS_NULL) {
1497 : /* LCOV_EXCL_START */
1498 : scs_special_print(print_mode, stderr, "ERROR: `R_prev` memory allocation failure\n");
1499 : return SCS_NULL;
1500 : /* LCOV_EXCL_STOP */
1501 : }
1502 326 : work->dir = scs_malloc(l * sizeof (scs_float));
1503 326 : if (work->dir == SCS_NULL) {
1504 : /* LCOV_EXCL_START */
1505 : scs_special_print(print_mode, stderr, "ERROR: `dir` memory allocation failure\n");
1506 : return SCS_NULL;
1507 : /* LCOV_EXCL_STOP */
1508 : }
1509 326 : work->dut = scs_malloc(l * sizeof (scs_float));
1510 326 : if (work->dut == SCS_NULL) {
1511 : /* LCOV_EXCL_START */
1512 : scs_special_print(print_mode, stderr, "ERROR: `dut` memory allocation failure\n");
1513 : return SCS_NULL;
1514 : /* LCOV_EXCL_STOP */
1515 : }
1516 326 : work->s_b = scs_malloc(data->m * sizeof (scs_float));
1517 326 : if (work->s_b == SCS_NULL) {
1518 : /* LCOV_EXCL_START */
1519 : scs_special_print(print_mode, stderr, "ERROR: `s_b` memory allocation failure\n");
1520 : return SCS_NULL;
1521 : /* LCOV_EXCL_STOP */
1522 : }
1523 :
1524 326 : work->stepsize = 1.0;
1525 :
1526 : /* -------------------------------------
1527 : * Restarted Broyden requires the allocation
1528 : * of an (S,U)-cache.
1529 : * ------------------------------------- */
1530 652 : if ((work->stgs->direction == restarted_broyden
1531 326 : || work->stgs->direction == anderson_acceleration)
1532 321 : && work->stgs->memory > 0) {
1533 321 : work->direction_cache = scs_init_direction_cache(work->stgs->memory, l, print_mode, work->stgs->direction);
1534 321 : if (work->direction_cache == SCS_NULL) {
1535 : /* LCOV_EXCL_START */
1536 : scs_special_print(print_mode, stderr,
1537 : "ERROR: `direction_cache` memory allocation failure\n");
1538 : return SCS_NULL;
1539 : /* LCOV_EXCL_STOP */
1540 : }
1541 : } else {
1542 5 : work->direction_cache = SCS_NULL;
1543 : }
1544 :
1545 : /* -------------------------------------
1546 : * Allocate memory for the full Broyden
1547 : * method
1548 : * ------------------------------------- */
1549 326 : if (work->stgs->direction == full_broyden) {
1550 : scs_int i;
1551 0 : work->H = scs_malloc(l * l * sizeof (scs_float));
1552 0 : if (work->H == SCS_NULL) {
1553 : /* LCOV_EXCL_START */
1554 : scs_special_print(print_mode, stderr, "ERROR: `H` memory allocation failure\n");
1555 : return SCS_NULL;
1556 : /* LCOV_EXCL_STOP */
1557 : }
1558 : /* H = I */
1559 0 : for (i = 0; i < l; ++i) {
1560 0 : work->H[i * (l + 1)] = 1.0;
1561 : }
1562 : } else {
1563 326 : work->H = SCS_NULL;
1564 : }
1565 :
1566 326 : work->Sk = scs_malloc(l * sizeof (scs_float));
1567 326 : if (work->Sk == SCS_NULL) {
1568 : /* LCOV_EXCL_START */
1569 : scs_special_print(print_mode, stderr, "ERROR: `Sk` memory allocation failure\n");
1570 : return SCS_NULL;
1571 : /* LCOV_EXCL_STOP */
1572 : }
1573 326 : work->Yk = scs_malloc(l * sizeof (scs_float));
1574 326 : if (work->Yk == SCS_NULL) {
1575 : /* LCOV_EXCL_START */
1576 : scs_special_print(print_mode, stderr, "ERROR: `Yk` memory allocation failure\n");
1577 : return SCS_NULL;
1578 : /* LCOV_EXCL_STOP */
1579 : }
1580 :
1581 326 : if (work->stgs->ls > 0) {
1582 322 : work->wu = scs_malloc(l * sizeof (scs_float));
1583 322 : if (work->wu == SCS_NULL) {
1584 : /* LCOV_EXCL_START */
1585 : scs_special_print(print_mode, stderr, "ERROR: `wu` memory allocation failure\n");
1586 : return SCS_NULL;
1587 : /* LCOV_EXCL_STOP */
1588 : }
1589 322 : work->Rwu = scs_malloc(l * sizeof (scs_float));
1590 322 : if (work->Rwu == SCS_NULL) {
1591 : /* LCOV_EXCL_START */
1592 : scs_special_print(print_mode, stderr, "ERROR: `Rwu` memory allocation failure\n");
1593 : return SCS_NULL;
1594 : /* LCOV_EXCL_STOP */
1595 : }
1596 322 : work->wu_t = scs_malloc(l * sizeof (scs_float));
1597 322 : if (work->wu_t == SCS_NULL) {
1598 : /* LCOV_EXCL_START */
1599 : scs_special_print(print_mode, stderr, "ERROR: `wu_t` memory allocation failure\n");
1600 : return SCS_NULL;
1601 : /* LCOV_EXCL_STOP */
1602 : }
1603 322 : work->wu_b = scs_malloc(l * sizeof (scs_float));
1604 322 : if (work->wu_b == SCS_NULL) {
1605 : /* LCOV_EXCL_START */
1606 : scs_special_print(print_mode, stderr, "ERROR: `wu_b` memory allocation failure\n");
1607 : return SCS_NULL;
1608 : /* LCOV_EXCL_STOP */
1609 : }
1610 : }
1611 : } else {
1612 : /* -------------------------------------
1613 : * In SCS, the pointers which correspond to
1614 : * SuperSCS are set to SCS_NULL and are
1615 : * inactive.
1616 : * ------------------------------------- */
1617 3 : work->R = SCS_NULL;
1618 3 : work->R_prev = SCS_NULL;
1619 3 : work->dir = SCS_NULL;
1620 3 : work->dut = SCS_NULL;
1621 3 : work->s_b = SCS_NULL;
1622 3 : work->direction_cache = SCS_NULL;
1623 3 : work->Yk = SCS_NULL;
1624 3 : work->Sk = SCS_NULL;
1625 3 : work->wu = SCS_NULL;
1626 3 : work->Rwu = SCS_NULL;
1627 3 : work->wu_t = SCS_NULL;
1628 3 : work->wu_b = SCS_NULL;
1629 : }
1630 :
1631 329 : work->A = data->A;
1632 329 : if (work->stgs->normalize) {
1633 : #ifdef COPYAMATRIX
1634 326 : if (scs_copy_a_matrix(&(work->A), data->A) == 0) {
1635 : /* LCOV_EXCL_START */
1636 : scs_special_print(print_mode, stderr, "ERROR: copy A matrix failed\n");
1637 : return SCS_NULL;
1638 : /* LCOV_EXCL_STOP */
1639 : }
1640 : #endif
1641 326 : work->scal = scs_malloc(sizeof (ScsScaling));
1642 326 : scs_normalize_a(work->A, work->stgs, cone, work->scal);
1643 : } else {
1644 3 : work->scal = SCS_NULL;
1645 : }
1646 329 : if ((work->coneWork = scs_init_conework(cone)) == SCS_NULL) {
1647 : /* LCOV_EXCL_START */
1648 : scs_special_print(print_mode, stderr, "ERROR: initCone failure\n");
1649 : return SCS_NULL;
1650 : /* LCOV_EXCL_STOP */
1651 : }
1652 329 : work->p = scs_init_priv(work->A, work->stgs);
1653 329 : if (work->p == SCS_NULL) {
1654 : /* LCOV_EXCL_START */
1655 : scs_special_print(print_mode, stderr, "ERROR: scs_init_priv failure\n");
1656 : return SCS_NULL;
1657 : /* LCOV_EXCL_STOP */
1658 : }
1659 : return work;
1660 : }
1661 :
1662 329 : static scs_int scs_update_work(
1663 : const ScsData * RESTRICT data,
1664 : ScsWork * RESTRICT work,
1665 : const ScsSolution * RESTRICT sol) {
1666 : /* before normalization */
1667 329 : scs_int n = data->n;
1668 329 : scs_int m = data->m;
1669 :
1670 329 : work->nm_b = scs_norm(data->b, m);
1671 329 : work->nm_c = scs_norm(data->c, n);
1672 329 : memcpy(work->b, data->b, data->m * sizeof (scs_float));
1673 329 : memcpy(work->c, data->c, data->n * sizeof (scs_float));
1674 :
1675 329 : if (work->stgs->normalize) {
1676 326 : scs_normalize_bc(work);
1677 : }
1678 329 : if (work->stgs->warm_start) {
1679 2 : scs_warm_start_vars(work, sol);
1680 : } else {
1681 327 : scs_cold_start_vars(work);
1682 : }
1683 329 : memcpy(work->h, work->c, n * sizeof (scs_float));
1684 329 : memcpy(&(work->h[n]), work->b, m * sizeof (scs_float));
1685 329 : memcpy(work->g, work->h, (n + m) * sizeof (scs_float));
1686 329 : scs_solve_lin_sys(work->A, work->stgs, work->p, work->g, SCS_NULL, -1);
1687 329 : scs_scale_array(&(work->g[n]), -1, m);
1688 329 : work->gTh = scs_inner_product(work->h, work->g, n + m);
1689 329 : return 0;
1690 : }
1691 :
1692 : static scs_int scs_validate_solve_input(
1693 : const ScsWork * RESTRICT work,
1694 : const ScsData * RESTRICT data,
1695 : const ScsCone * RESTRICT cone,
1696 : const ScsSolution * RESTRICT sol,
1697 : const ScsInfo * RESTRICT info) {
1698 329 : return ((data == SCS_NULL
1699 329 : || cone == SCS_NULL
1700 329 : || sol == SCS_NULL
1701 329 : || info == SCS_NULL
1702 329 : || work == SCS_NULL
1703 329 : || data->b == SCS_NULL
1704 658 : || data->c == SCS_NULL));
1705 : }
1706 :
1707 329 : static scs_int scs_init_progress_data(
1708 : ScsInfo * RESTRICT info,
1709 : ScsWork * RESTRICT work) {
1710 : /* --------------------------------------------------------------
1711 : * As the might be successive calls to superSCS (superscs_solve)
1712 : * with different values of `do_record_progress`, we nee to make
1713 : * sure that we neither allocate memory for the same array twice,
1714 : * nor do we set any arrays to SCS_NULL unnecessarily. Note that,
1715 : * provided that `info` has been initialized with scs_init_info (this
1716 : * is highly recommended), all pointers are initially set to SCS_NULL.
1717 : * -------------------------------------------------------------- */
1718 329 : if (work->stgs->do_record_progress) {
1719 8 : const scs_int max_history_alloc = work->stgs->max_iters;
1720 :
1721 : /* ----------------------------------------
1722 : * If a pointer is SCS_NULL, it means that
1723 : * no memory has been allocated for that
1724 : * previously.
1725 : * ---------------------------------------- */
1726 8 : if (info->progress_relgap == SCS_NULL) {
1727 3 : info->progress_relgap = scs_malloc(sizeof (scs_float) * max_history_alloc);
1728 3 : if (info->progress_relgap == SCS_NULL) return -1;
1729 : }
1730 8 : if (info->progress_respri == SCS_NULL) {
1731 3 : info->progress_respri = scs_malloc(sizeof (scs_float) * max_history_alloc);
1732 3 : if (info->progress_respri == SCS_NULL) return -2;
1733 : }
1734 8 : if (info->progress_resdual == SCS_NULL) {
1735 3 : info->progress_resdual = scs_malloc(sizeof (scs_float) * max_history_alloc);
1736 3 : if (info->progress_resdual == SCS_NULL) return -3;
1737 : }
1738 8 : if (info->progress_pcost == SCS_NULL) {
1739 3 : info->progress_pcost = scs_malloc(sizeof (scs_float) * max_history_alloc);
1740 3 : if (info->progress_pcost == SCS_NULL) return -4;
1741 : }
1742 8 : if (info->progress_dcost == SCS_NULL) {
1743 3 : info->progress_dcost = scs_malloc(sizeof (scs_float) * max_history_alloc);
1744 3 : if (info->progress_dcost == SCS_NULL) return -5;
1745 : }
1746 8 : if (info->progress_iter == SCS_NULL) {
1747 3 : info->progress_iter = scs_malloc(sizeof (scs_int) * max_history_alloc);
1748 3 : if (info->progress_iter == SCS_NULL) return -6;
1749 : }
1750 8 : if (info->progress_norm_fpr == SCS_NULL) {
1751 3 : info->progress_norm_fpr = scs_malloc(sizeof (scs_float) * max_history_alloc);
1752 3 : if (info->progress_norm_fpr == SCS_NULL) return -7;
1753 : }
1754 8 : if (info->progress_time == SCS_NULL) {
1755 3 : info->progress_time = scs_malloc(sizeof (scs_float) * max_history_alloc);
1756 3 : if (info->progress_time == SCS_NULL) return -8;
1757 : }
1758 8 : if (info->progress_mode == SCS_NULL) {
1759 3 : info->progress_mode = scs_malloc(sizeof (scs_int) * max_history_alloc);
1760 3 : if (info->progress_mode == SCS_NULL) return -9;
1761 : }
1762 8 : if (info->progress_ls == SCS_NULL) {
1763 3 : info->progress_ls = scs_malloc(sizeof (scs_int) * max_history_alloc);
1764 3 : if (info->progress_ls == SCS_NULL) return -10;
1765 : }
1766 :
1767 : /* ---------------------------------------------------------
1768 : * If `do_record_progress` is true, and there has
1769 : * been a previous allocation, but now the maximum
1770 : * number of iterations has increased
1771 : *
1772 : * Note: if the current `max_iters` is smaller than
1773 : * the previous value, it means that more than adequate
1774 : * memory space has been allocated for the progress arrays.
1775 : * However, we will not use `realloc` to size it down.
1776 : * --------------------------------------------------------- */
1777 8 : if (work->stgs->previous_max_iters != -1
1778 6 : && max_history_alloc > work->stgs->previous_max_iters) {
1779 : /* ------------------------------------
1780 : * We don't check for NULL values here
1781 : * because `realloc` on NULL pointers
1782 : * behaves like `malloc`
1783 : * ------------------------------------
1784 : */
1785 5 : info->progress_relgap = realloc(info->progress_relgap, sizeof (scs_float) * max_history_alloc);
1786 5 : if (info->progress_relgap == SCS_NULL) return -100;
1787 :
1788 5 : info->progress_respri = realloc(info->progress_respri, sizeof (scs_float) * max_history_alloc);
1789 5 : if (info->progress_respri == SCS_NULL) return -101;
1790 :
1791 5 : info->progress_resdual = realloc(info->progress_resdual, sizeof (scs_float) * max_history_alloc);
1792 5 : if (info->progress_resdual == SCS_NULL) return -102;
1793 :
1794 5 : info->progress_pcost = realloc(info->progress_pcost, sizeof (scs_float) * max_history_alloc);
1795 5 : if (info->progress_pcost == SCS_NULL) return -103;
1796 :
1797 5 : info->progress_dcost = realloc(info->progress_dcost, sizeof (scs_float) * max_history_alloc);
1798 5 : if (info->progress_dcost == SCS_NULL) return -104;
1799 :
1800 5 : info->progress_iter = realloc(info->progress_iter, sizeof (scs_int) * max_history_alloc);
1801 5 : if (info->progress_iter == SCS_NULL) return -105;
1802 :
1803 5 : info->progress_norm_fpr = realloc(info->progress_norm_fpr, sizeof (scs_float) * max_history_alloc);
1804 5 : if (info->progress_norm_fpr == SCS_NULL) return -106;
1805 :
1806 5 : info->progress_time = realloc(info->progress_time, sizeof (scs_float) * max_history_alloc);
1807 5 : if (info->progress_time == SCS_NULL) return -107;
1808 :
1809 5 : info->progress_mode = realloc(info->progress_mode, sizeof (scs_int) * max_history_alloc);
1810 5 : if (info->progress_mode == SCS_NULL) return -108;
1811 :
1812 5 : info->progress_ls = realloc(info->progress_ls, sizeof (scs_int) * max_history_alloc);
1813 5 : if (info->progress_ls == SCS_NULL) return -109;
1814 : }
1815 : }
1816 : return 0;
1817 : }
1818 :
1819 987 : static void scs_record_progress_data(
1820 : ScsInfo * info,
1821 : struct scs_residuals * res,
1822 : const ScsWork * work,
1823 : ScsTimer * solveTimer,
1824 : scs_int iter) {
1825 987 : scs_int idx_progress = iter / SCS_CONVERGED_INTERVAL;
1826 987 : info->progress_iter[idx_progress] = iter;
1827 987 : info->progress_relgap[idx_progress] = res->rel_gap;
1828 987 : info->progress_respri[idx_progress] = res->res_pri;
1829 987 : info->progress_resdual[idx_progress] = res->res_dual;
1830 987 : info->progress_pcost[idx_progress] = res->cTx_by_tau / res->tau;
1831 987 : info->progress_dcost[idx_progress] = -res->bTy_by_tau / res->tau;
1832 987 : info->progress_norm_fpr[idx_progress] = work->nrmR_con;
1833 987 : info->progress_time[idx_progress] = scs_toc_quiet(solveTimer);
1834 987 : }
1835 :
1836 3 : scs_int scs_solve(
1837 200 : ScsWork * RESTRICT work,
1838 : const ScsData * RESTRICT data,
1839 : const ScsCone * RESTRICT cone,
1840 : ScsSolution * RESTRICT sol,
1841 : ScsInfo * RESTRICT info) {
1842 : scs_int i;
1843 : ScsTimer solveTimer;
1844 : struct scs_residuals r;
1845 3 : scs_int print_mode = work->stgs->do_override_streams;
1846 3 : const scs_float max_runtime_millis = work->stgs->max_time_milliseconds;
1847 :
1848 3 : if ((i = scs_init_progress_data(info, work)) < 0) {
1849 : /* LCOV_EXCL_START */
1850 : scs_special_print(print_mode, stderr, "Memory allocation error (progress arrays), code: %d\n", (int) i);
1851 : return SCS_FAILED;
1852 : /* LCOV_EXCL_STOP */
1853 : }
1854 :
1855 3 : if (scs_validate_solve_input(work, data, cone, sol, info)) {
1856 0 : scs_special_print(print_mode, stderr, "ERROR: SCS_NULL input\n");
1857 0 : return SCS_FAILED;
1858 : }
1859 : /* initialize ctrl-c support */
1860 3 : startInterruptListener();
1861 3 : scs_tic(&solveTimer);
1862 3 : info->statusVal = SCS_UNFINISHED; /* not yet converged */
1863 3 : r.last_iter = -1;
1864 3 : scs_update_work(data, work, sol);
1865 :
1866 3 : if (work->stgs->verbose)
1867 0 : scs_print_header(work, cone);
1868 : /* scs: */
1869 198 : for (i = 0; i < work->stgs->max_iters && scs_toc_quiet(&solveTimer) < max_runtime_millis; ++i) {
1870 200 : memcpy(work->u_prev, work->u, work->l * sizeof (scs_float));
1871 :
1872 200 : if (scs_project_lin_sys(work, i) < 0) {
1873 0 : return scs_failure(work, work->m, work->n, sol, info, SCS_FAILED,
1874 : "error in projectLinSys", "Failure", print_mode);
1875 : }
1876 200 : if (scs_project_cones(work, cone, i) < 0) {
1877 0 : return scs_failure(work, work->m, work->n, sol, info, SCS_FAILED,
1878 : "error in projectCones", "Failure", print_mode);
1879 : }
1880 :
1881 200 : scs_update_dual_vars(work);
1882 :
1883 200 : if (isInterrupted()) {
1884 0 : return scs_failure(work, work->m, work->n, sol, info, SCS_SIGINT, "Interrupted",
1885 : "Interrupted", print_mode);
1886 : }
1887 : if (i % SCS_CONVERGED_INTERVAL == 0) {
1888 200 : scs_calc_residuals(work, &r, i);
1889 200 : if (work->stgs->do_record_progress) scs_record_progress_data(info, &r, work, &solveTimer, i);
1890 400 : if ((info->statusVal = scs_has_converged(work, &r, i)) != 0) break;
1891 : }
1892 :
1893 198 : if (work->stgs->verbose && i % SCS_PRINT_INTERVAL == 0) {
1894 0 : scs_calc_residuals(work, &r, i);
1895 0 : scs_print_summary(work, i, &r, &solveTimer);
1896 : }
1897 : }
1898 3 : if (work->stgs->verbose) {
1899 0 : scs_calc_residuals(work, &r, i);
1900 0 : scs_print_summary(work, i, &r, &solveTimer);
1901 : }
1902 : /* populate solution vectors (unnormalized) and info */
1903 3 : scs_get_solution(work, sol, info, &r, i);
1904 3 : info->solveTime = scs_toc_quiet(&solveTimer);
1905 :
1906 3 : if (work->stgs->verbose)
1907 : scs_print_footer(data, cone, sol, work, info); /* LCOV_EXCL_LINE */
1908 3 : endInterruptListener();
1909 3 : info->history_length = i / SCS_CONVERGED_INTERVAL;
1910 :
1911 3 : return info->statusVal;
1912 : }
1913 :
1914 29901 : static void scs_compute_sb_kapb(
1915 : ScsWork * RESTRICT work) {
1916 29901 : scs_axpy(work->s_b, work->u_b + work->n, work->u_t + work->n, 1.0, -2.0, work->m);
1917 29901 : scs_add_array(work->s_b, work->u + work->n, work->m);
1918 29901 : work->kap_b = work->u_b[work->l - 1] - 2.0 * work->u_t[work->l - 1] + work->u[work->l - 1];
1919 29901 : }
1920 :
1921 5097 : static void scs_step_k1(
1922 : ScsWork * RESTRICT work,
1923 : scs_float nrmRw_con_current,
1924 : scs_float * r_safe,
1925 : scs_float q_to_power_iter_times_nrm_R_init,
1926 : scs_int * how) {
1927 5097 : memcpy(work->u, work->wu, work->l * sizeof (scs_float)); /* u = wu */
1928 5097 : memcpy(work->u_t, work->wu_t, work->l * sizeof (scs_float)); /* u_t = wu_t */
1929 5097 : memcpy(work->u_b, work->wu_b, work->l * sizeof (scs_float)); /* u_b = wu_b */
1930 5097 : memcpy(work->R, work->Rwu, work->l * sizeof (scs_float)); /* R = Rw */
1931 5097 : scs_compute_sb_kapb(work);
1932 5097 : work->nrmR_con = nrmRw_con_current;
1933 5097 : *r_safe = work->nrmR_con + q_to_power_iter_times_nrm_R_init;
1934 5097 : *how = (scs_int) 1;
1935 5097 : }
1936 :
1937 24944 : static scs_int scs_step_k2(
1938 : ScsWork * RESTRICT work,
1939 : scs_float nrmRw_con,
1940 : scs_int * how) {
1941 :
1942 24944 : scs_int do_break_loop = 0;
1943 : scs_float slack;
1944 : scs_float rhs;
1945 49888 : slack = nrmRw_con * nrmRw_con - work->stepsize * (
1946 24944 : scs_inner_product(work->dir + work->n, work->Rwu + work->n, work->m + 1)
1947 24944 : + work->stgs->rho_x * scs_inner_product(work->dir, work->Rwu, work->n));
1948 24944 : rhs = work->stgs->sigma * work->nrmR_con * nrmRw_con;
1949 24944 : if (slack >= rhs) {
1950 : scs_float stepsize2;
1951 12563 : stepsize2 = (work->stgs->alpha * (slack / (nrmRw_con * nrmRw_con)));
1952 12563 : scs_add_scaled_array(work->u, work->Rwu, work->l, -stepsize2);
1953 12563 : *how = 2;
1954 12563 : do_break_loop = 1;
1955 : }
1956 24944 : return do_break_loop;
1957 : }
1958 :
1959 29120 : static void scs_update_caches(
1960 : ScsWork * RESTRICT work,
1961 : scs_float sqrt_rhox,
1962 : scs_int how) {
1963 29120 : if (how == 0 || work->stgs->ls == 0) {
1964 5310 : scs_axpy(work->Sk, work->u, work->u_prev, 1.0, -1.0, work->l); /* Sk = u - u_prev */
1965 5310 : scs_axpy(work->Yk, work->R, work->R_prev, sqrt_rhox, -1.0, work->n); /* Yk = sqrt_rhox * R - R_prev */
1966 5310 : scs_axpy(work->Yk + work->n, work->R + work->n, work->R_prev + work->n, 1.0, -1.0, work->m + 1);
1967 5310 : scs_scale_array(work->Sk, sqrt_rhox, work->n); /* Sk *= sqrt_rhox */
1968 : } else {
1969 23810 : scs_axpy(work->Sk, work->wu, work->u_prev, sqrt_rhox, -sqrt_rhox, work->n);
1970 23810 : scs_axpy(work->Sk + work->n, work->wu + work->n, work->u_prev + work->n, 1.0, -1.0, work->m + 1);
1971 23810 : scs_axpy(work->Yk, work->Rwu, work->R_prev, sqrt_rhox, -1.0, work->n);
1972 23810 : scs_axpy(work->Yk + work->n, work->Rwu + work->n, work->R_prev + work->n, 1.0, -1.0, work->m + 1);
1973 : }
1974 29120 : }
1975 :
1976 24478 : static scs_int scs_exit_loop_without_k1(
1977 : ScsWork * RESTRICT work,
1978 : ScsSolution * RESTRICT sol,
1979 : ScsInfo * RESTRICT info,
1980 : const ScsCone * RESTRICT cone,
1981 : scs_int i,
1982 : scs_int print_mode) {
1983 24478 : if (superscs_project_lin_sys(work->u_t, work->u, work, i) < 0) {
1984 0 : return scs_failure(work, work->m, work->n, sol, info, SCS_FAILED,
1985 : "error in projectLinSysv2", "Failure", print_mode);
1986 : }
1987 : /* u_bar = proj_C(2u_t - u) */
1988 24478 : if (superscs_project_cones(work->u_b, work->u_t, work->u, work, cone, i) < 0) {
1989 0 : return scs_failure(work, work->m, work->n, sol, info, SCS_FAILED,
1990 : "error in projectConesv2", "Failure", print_mode);
1991 : }
1992 24478 : scs_compute_sb_kapb(work);
1993 24478 : scs_calc_FPR(work->R, work->u_t, work->u_b, work->l);
1994 24478 : work->nrmR_con = SQRTF(
1995 24478 : work->stgs->rho_x * scs_norm_squared(work->R, work->n)
1996 24478 : + scs_norm_squared(work->R + work->n, work->m + 1)
1997 : );
1998 24478 : return 0;
1999 : }
2000 :
2001 326 : scs_int superscs_solve(
2002 29888 : ScsWork * RESTRICT work,
2003 : const ScsData * RESTRICT data,
2004 : const ScsCone * RESTRICT cone,
2005 : ScsSolution * RESTRICT sol,
2006 : ScsInfo * RESTRICT info) {
2007 : scs_int i; /* i indexes the (outer) iterations */
2008 326 : scs_int how = 0; /* -1:unsuccessful backtracking, 0:K0, 1:K1, 2:K2 */
2009 : scs_float eta;
2010 : scs_float nrm_R_0;
2011 : scs_float r_safe;
2012 : scs_float nrmRw_con; /* norm of FP res at line-search */
2013 : scs_float nrmR_con_old; /* keeps previous FP res */
2014 326 : scs_float q_to_power_iter_times_nrmR0 = work->stgs->sse; /* sse^(iter) */
2015 326 : const scs_float rhox = work->stgs->rho_x;
2016 326 : const scs_float sqrt_rhox = SQRTF(rhox);
2017 326 : const scs_float q0 = work->stgs->sse;
2018 326 : const scs_int n = work->n;
2019 326 : const scs_int m = work->m;
2020 326 : const scs_int l = work->l;
2021 : ScsTimer solveTimer;
2022 : struct scs_residuals r;
2023 326 : scs_int print_mode = work->stgs->do_override_streams;
2024 : /* ------------------------------------
2025 : * Store some pointers in function-scope
2026 : * variables for performance.
2027 : * ------------------------------------ */
2028 326 : ScsSettings * RESTRICT settings = work->stgs;
2029 326 : scs_float alpha = settings->alpha;
2030 326 : scs_float * RESTRICT dir = work->dir;
2031 326 : scs_float * RESTRICT R = work->R;
2032 326 : scs_float * RESTRICT R_prev = work->R_prev;
2033 326 : scs_float * RESTRICT Rwu = work->Rwu;
2034 326 : scs_float * RESTRICT u = work->u;
2035 326 : scs_float * RESTRICT u_t = work->u_t;
2036 326 : scs_float * RESTRICT u_b = work->u_b;
2037 326 : scs_float * RESTRICT u_prev = work->u_prev;
2038 326 : scs_float * RESTRICT wu = work->wu;
2039 326 : scs_float * RESTRICT wu_t = work->wu_t;
2040 326 : scs_float * RESTRICT wu_b = work->wu_b;
2041 326 : scs_float * RESTRICT dut = work->dut;
2042 :
2043 326 : if ((i = scs_init_progress_data(info, work)) < 0) {
2044 : /* LCOV_EXCL_START */
2045 : scs_special_print(print_mode, stderr, "Memory allocation error (progress arrays), code: %d\n", (int) i);
2046 : return SCS_FAILED;
2047 : /* LCOV_EXCL_STOP */
2048 : }
2049 326 : settings->previous_max_iters = settings->max_iters;
2050 :
2051 326 : if (scs_validate_solve_input(work, data, cone, sol, info)) {
2052 0 : scs_special_print(print_mode, stderr, "ERROR: SCS_NULL input\n");
2053 0 : return SCS_FAILED;
2054 : }
2055 :
2056 : /* initialize ctrl-c support */
2057 326 : startInterruptListener();
2058 326 : scs_tic(&solveTimer);
2059 326 : info->statusVal = SCS_UNFINISHED; /* not yet converged */
2060 326 : r.last_iter = -1;
2061 326 : scs_update_work(data, work, sol);
2062 :
2063 326 : if (settings->verbose > 0) scs_print_header(work, cone);
2064 :
2065 : /* Initialize: */
2066 326 : i = 0; /* Needed for the next two functions */
2067 326 : if (superscs_project_lin_sys(u_t, u, work, i) < 0) { /* u_t = (I+Q)^{-1} u*/
2068 0 : return scs_failure(work, m, n, sol, info, SCS_FAILED,
2069 : "error in projectLinSysv2", "Failure", print_mode);
2070 : }
2071 326 : if (superscs_project_cones(u_b, u_t, u, work, cone, i) < 0) { /* u_bar = proj_C(2u_t - u) */
2072 0 : return scs_failure(work, m, n, sol, info, SCS_FAILED,
2073 : "error in projectConesv2", "Failure", print_mode);
2074 : }
2075 326 : scs_compute_sb_kapb(work); /* compute s_b and kappa_b */
2076 : scs_calc_FPR(R, u_t, u_b, l); /* compute Ru */
2077 : /* initialize eta = |Ru^0| (norm of R using rho_x)... */
2078 326 : eta = SQRTF(rhox * scs_norm_squared(R, n) + scs_norm_squared(R + n, m + 1));
2079 326 : r_safe = eta;
2080 326 : work->nrmR_con = eta;
2081 326 : nrm_R_0 = MIN(1.0, eta);
2082 326 : q_to_power_iter_times_nrmR0 *= nrm_R_0;
2083 :
2084 : /* +++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2085 : *
2086 : * MAIN SUPER SCS LOOP
2087 : *
2088 : * ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ */
2089 30227 : for (i = 0; i < settings->max_iters
2090 59464 : && scs_toc_quiet(&solveTimer) < work->stgs->max_time_milliseconds; ++i) {
2091 29888 : scs_int j_iter_ls = 0; /* j indexes the line search iterations */
2092 :
2093 29888 : if (isInterrupted()) {
2094 0 : return scs_failure(work, m, n, sol, info, SCS_SIGINT, "Interrupted",
2095 : "Interrupted", print_mode);
2096 : }
2097 :
2098 : /* Convergence checks */
2099 : if (i % SCS_CONVERGED_INTERVAL == 0) {
2100 29888 : scs_calc_residuals_superscs(work, &r, i);
2101 29888 : if (settings->do_record_progress) scs_record_progress_data(info, &r, work, &solveTimer, i);
2102 59776 : if ((info->statusVal = scs_has_converged(work, &r, i)) != 0) break;
2103 : }
2104 :
2105 : /* Prints results every PRINT_INTERVAL iterations */
2106 29575 : if (settings->verbose && i % SCS_PRINT_INTERVAL == 0) {
2107 2 : scs_calc_residuals_superscs(work, &r, i);
2108 2 : scs_print_summary(work, i, &r, &solveTimer);
2109 : }
2110 :
2111 29575 : if (settings->ls > 0 || settings->k0 == 1) {
2112 29441 : q_to_power_iter_times_nrmR0 *= q0; /* q = q0^i = sse^i */
2113 29441 : if (i == 0) {
2114 : /* -------------------------------------------
2115 : * At i=0, the direction is defined using the
2116 : * FPR: dir^0 = -R
2117 : * -------------------------------------------- */
2118 321 : scs_set_as_scaled_array(dir, R, -sqrt_rhox, n);
2119 321 : scs_set_as_scaled_array(dir + n, R + n, -1, m + 1);
2120 :
2121 : } else {
2122 29120 : scs_update_caches(work, sqrt_rhox, how);
2123 29120 : scs_scale_array(R, sqrt_rhox, n); /* R *= sqrt_rhox */
2124 29120 : if (scs_compute_direction(work, i) < 0) { /* compute direction */
2125 : /*~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2126 : * Function `scs_compute_direction` is invoked at iterations i>=1.
2127 : * At i=1, the first pair (S,Y) has been computed (see work).
2128 : *~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~*/
2129 0 : return scs_failure(work, m, n, sol, info, SCS_FAILED,
2130 : "error in scs_compute_direction", "Failure", print_mode);
2131 : }
2132 29120 : scs_scale_array(R, 1 / sqrt_rhox, n); /* R = R/sqrt_rhox */
2133 : }
2134 : /* -------------------------------------------
2135 : * Scale the x-part of dir using sqrt_rhox
2136 : * -------------------------------------------- */
2137 29441 : scs_scale_array(dir, 1 / sqrt_rhox, n);
2138 : }
2139 :
2140 29575 : memcpy(u_prev, u, l * sizeof (scs_float)); /* u_prev = u */
2141 29575 : memcpy(R_prev, R, l * sizeof (scs_float)); /* R_prev = R */
2142 29575 : scs_scale_array(R_prev, sqrt_rhox, n);
2143 29575 : how = -1; /* no backtracking (yet) */
2144 29575 : nrmR_con_old = work->nrmR_con;
2145 :
2146 29575 : if (i >= settings->warm_start) {
2147 : /* ------------------------------------------------
2148 : * Blind updates (K0)
2149 : * ------------------------------------------------ */
2150 29575 : if (settings->k0 == 1 && work->nrmR_con <= settings->c_bl * eta) {
2151 5466 : scs_add_array(u, dir, l); /* u += dir */
2152 5466 : how = 0;
2153 5466 : eta = work->nrmR_con;
2154 5466 : work->stepsize = 1.0;
2155 24109 : } else if (settings->ls > 0) {
2156 23975 : if (superscs_project_lin_sys(dut, dir, work, i) < 0) {
2157 0 : return scs_failure(work, m, n, sol, info, SCS_FAILED, "error in superscs_project_lin_sys", "Failure", print_mode);
2158 : }
2159 23975 : work->stepsize = 2.0;
2160 :
2161 : /* ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2162 : * Line search
2163 : * Main computational burden: 1 projection
2164 : * ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ */
2165 60969 : for (j_iter_ls = 0; j_iter_ls < settings->ls; ++j_iter_ls) {
2166 54654 : work->stepsize *= settings->beta;
2167 54654 : scs_axpy(wu, u, dir, 1.0, work->stepsize, l); /* wu = u + step * dir */
2168 54654 : scs_axpy(wu_t, u_t, dut, 1.0, work->stepsize, l); /* wut = u_t + step * dut */
2169 :
2170 54654 : if (superscs_project_cones(wu_b, wu_t, wu, work, cone, i) < 0) {
2171 0 : return scs_failure(work, m, n, sol, info, SCS_FAILED, "error in superscs_project_cones", "Failure", print_mode);
2172 : }
2173 : scs_calc_FPR(Rwu, wu_t, wu_b, l); /* calculate FPR on scaled vectors */
2174 :
2175 54654 : nrmRw_con = SQRTF(scs_norm_squared(Rwu + n, m + 1) + rhox * scs_norm_squared(Rwu, n));
2176 :
2177 : /* K1 */
2178 54654 : if (settings->k1
2179 21108 : && nrmRw_con <= settings->c1 * nrmR_con_old
2180 7050 : && work->nrmR_con <= r_safe) {
2181 5097 : scs_step_k1(work, nrmRw_con, &r_safe, q_to_power_iter_times_nrmR0, &how);
2182 5097 : break;
2183 : }
2184 :
2185 : /* K2 */
2186 49557 : if (settings->k2 && scs_step_k2(work, nrmRw_con, &how)) break;
2187 : } /* end of line-search */
2188 23975 : j_iter_ls++; /* to get the number of LS iterations */
2189 : } /* end of `else if` block (when !K1 OR no blind update) */
2190 : } /* IF-block: iterated after warm start */
2191 :
2192 29575 : if (how == -1) { /* means that R didn't change */
2193 : /* x -= alpha*Rx */
2194 6449 : scs_add_scaled_array(u, R, l, -alpha);
2195 : } /* how == -1 */
2196 29575 : if (how != 1) { /* exited with other than K1 */
2197 24478 : scs_int status = scs_exit_loop_without_k1(work, sol, info, cone, i, print_mode);
2198 24478 : if (status < 0)
2199 : return status;
2200 : } /* how != 1 */
2201 :
2202 : /* -------------------------------------------
2203 : * Record some more progress information
2204 : * -------------------------------------------*/
2205 29575 : if (settings->do_record_progress && i % SCS_CONVERGED_INTERVAL == 0) {
2206 905 : scs_int idx_progress = i / SCS_CONVERGED_INTERVAL;
2207 905 : info->progress_mode[idx_progress] = how;
2208 905 : info->progress_ls[idx_progress] = j_iter_ls;
2209 : }
2210 :
2211 : } /* main for loop */
2212 :
2213 :
2214 326 : info->cg_total_iters = scs_linsys_total_cg_iters(work->p);
2215 326 : info->linsys_total_solve_time_ms = scs_linsys_total_solve_time_ms(work->p);
2216 :
2217 326 : scs_calc_residuals_superscs(work, &r, i);
2218 :
2219 : /* prints summary of last iteration */
2220 326 : if (settings->verbose) scs_print_summary(work, i, &r, &solveTimer);
2221 :
2222 : /* populate solution vectors (unnormalized) and info */
2223 : /* update u, denormalize, etc */
2224 326 : scs_get_solution(work, sol, info, &r, i);
2225 326 : info->iter = i;
2226 326 : info->solveTime = scs_toc_quiet(&solveTimer);
2227 326 : info->history_length = i / SCS_CONVERGED_INTERVAL;
2228 :
2229 326 : if (settings->verbose)
2230 : scs_print_footer(data, cone, sol, work, info); /* LCOV_EXCL_LINE */
2231 326 : endInterruptListener();
2232 :
2233 326 : return info->statusVal;
2234 : }
2235 :
2236 341 : void scs_finish(ScsWork * RESTRICT w) {
2237 341 : if (w != SCS_NULL) {
2238 329 : scs_finish_cone(w->coneWork);
2239 329 : if (w->stgs != SCS_NULL && w->stgs->normalize != 0) {
2240 : #ifndef COPYAMATRIX
2241 : scs_unnormalize_a(w->A, w->stgs, w->scal);
2242 : #else
2243 326 : scs_free_a_matrix(w->A);
2244 : #endif
2245 : }
2246 329 : if (w->p != SCS_NULL)
2247 329 : scs_free_priv(w->p);
2248 329 : scs_free_work(w);
2249 : }
2250 341 : }
2251 :
2252 341 : ScsWork * scs_init(
2253 : const ScsData * RESTRICT data,
2254 : const ScsCone * RESTRICT cone,
2255 : ScsInfo * RESTRICT info) {
2256 : ScsWork * RESTRICT work;
2257 : ScsTimer initTimer;
2258 341 : startInterruptListener();
2259 682 : if (data == SCS_NULL
2260 341 : || cone == SCS_NULL
2261 341 : || info == SCS_NULL) {
2262 : /* LCOV_EXCL_START */
2263 : scs_special_print(data->stgs->do_override_streams, stderr, "ERROR: Missing Data, Cone or Info input\n");
2264 : return SCS_NULL;
2265 : /* LCOV_EXCL_STOP */
2266 : }
2267 : #ifndef NOVALIDATE
2268 341 : if (scs_validate(data, cone) < 0) {
2269 12 : scs_special_print(data->stgs->do_override_streams, stderr, "ERROR: Validation returned failure\n");
2270 12 : return SCS_NULL;
2271 : }
2272 : #endif
2273 329 : scs_tic(&initTimer);
2274 329 : work = scs_init_work(data, cone);
2275 : /* strtoc("init", &initTimer); */
2276 329 : info->setupTime = scs_toc_quiet(&initTimer);
2277 329 : if (data->stgs->verbose) {
2278 4 : scs_special_print(work->stgs->do_override_streams,
2279 2 : work->stgs->output_stream, "Setup time: %1.2es\n", info->setupTime / 1e3);
2280 : }
2281 329 : endInterruptListener();
2282 329 : return work;
2283 : }
2284 :
2285 2 : static void scs_compute_allocated_memory(
2286 : const ScsWork * RESTRICT work,
2287 : const ScsCone * RESTRICT k,
2288 : const ScsData * RESTRICT data,
2289 : ScsInfo * RESTRICT info) {
2290 2 : blasint temp_n_max = 0;
2291 : unsigned long long allocated_memory;
2292 : size_t i;
2293 2 : size_t float_size = (scs_int) sizeof (scs_float);
2294 2 : size_t int_size = (scs_int) sizeof (scs_int);
2295 2 : size_t l = data->m + data->n + 1;
2296 2 : long mem = work->stgs->memory;
2297 :
2298 2 : for (i = 0; i < k->ssize; ++i) {
2299 0 : if (k->s[i] > temp_n_max) {
2300 0 : temp_n_max = (blasint) k->s[i];
2301 : }
2302 : }
2303 :
2304 2 : allocated_memory =
2305 4 : float_size * (2 * data->A->p[data->A->n]
2306 2 : + 6 * data->m
2307 2 : + 9 * data->n
2308 2 : + k->qsize
2309 2 : + k->psize
2310 2 : + k->ssize
2311 2 : + 2 * temp_n_max * temp_n_max
2312 2 : + temp_n_max
2313 2 : + work->coneWork->lwork
2314 2 : + 10 * l)
2315 4 : + int_size * (2 * data->A->p[data->A->n]
2316 2 : + data->n
2317 2 : + work->coneWork->liwork
2318 2 : + data->m + 2);
2319 :
2320 2 : if (work->stgs->ls > 0) {
2321 2 : allocated_memory += float_size * 4 * l;
2322 : }
2323 2 : if (work->stgs->direction == restarted_broyden && mem > 0) {
2324 2 : allocated_memory += float_size * 2 * l * (mem + 1);
2325 : }
2326 2 : if (work->stgs->direction == anderson_acceleration) {
2327 0 : allocated_memory += float_size * (4 * l * mem + l
2328 0 : + mem + scs_svd_workspace_size(l, mem));
2329 : }
2330 :
2331 2 : if (work->stgs->normalize) {
2332 1 : allocated_memory += float_size * (data->m + data->n);
2333 : }
2334 2 : info->allocated_memory = allocated_memory;
2335 2 : }
2336 :
2337 : #define SCS_BYTE_UNIT_STR_LEN 3
2338 :
2339 2 : static void scs_units_allocated_mem(
2340 : const ScsInfo * RESTRICT info,
2341 : scs_float * memory,
2342 : char *unit_of_measurement) {
2343 2 : if (info->allocated_memory > 1e9) {
2344 : strncpy(unit_of_measurement, "GB", SCS_BYTE_UNIT_STR_LEN);
2345 0 : *memory = info->allocated_memory / 1e9;
2346 2 : } else if (info->allocated_memory > 1e6) {
2347 : strncpy(unit_of_measurement, "MB", SCS_BYTE_UNIT_STR_LEN);
2348 0 : *memory = info->allocated_memory / 1e6;
2349 2 : } else if (info->allocated_memory > 1e3) {
2350 : strncpy(unit_of_measurement, "kB", SCS_BYTE_UNIT_STR_LEN);
2351 2 : *memory = info->allocated_memory / 1e3;
2352 : } else {
2353 : strncpy(unit_of_measurement, "B", SCS_BYTE_UNIT_STR_LEN);
2354 0 : *memory = info->allocated_memory;
2355 : }
2356 2 : unit_of_measurement[SCS_BYTE_UNIT_STR_LEN - 1] = '\0';
2357 2 : }
2358 :
2359 : #define SCS_DIRECTION_STRING_LENGTH 24
2360 :
2361 1 : static void scs_print_parameter_details(const ScsData * RESTRICT data) {
2362 1 : scs_int print_mode = data->stgs->do_override_streams;
2363 1 : FILE * stream = data->stgs->output_stream;
2364 : char dir_string[SCS_DIRECTION_STRING_LENGTH];
2365 1 : switch (data->stgs->direction) {
2366 : case anderson_acceleration:
2367 : strncpy(dir_string, "anderson", SCS_DIRECTION_STRING_LENGTH);
2368 : break;
2369 : case restarted_broyden:
2370 : strncpy(dir_string, "restarted broyden", SCS_DIRECTION_STRING_LENGTH);
2371 : break;
2372 : case fixed_point_residual:
2373 : strncpy(dir_string, "fixed point residual", SCS_DIRECTION_STRING_LENGTH);
2374 : break;
2375 : case full_broyden:
2376 : strncpy(dir_string, "full broyden", SCS_DIRECTION_STRING_LENGTH);
2377 : break;
2378 : default:
2379 : strncpy(dir_string, "unknown", SCS_DIRECTION_STRING_LENGTH);
2380 : }
2381 3 : scs_special_print(
2382 : print_mode,
2383 : stream,
2384 : "\nSettings:\n"
2385 : ".....................................................................\n"
2386 : "alpha : %2.1f\t\t"
2387 : "beta : %2.1f\n"
2388 : "c1 : %2.1f\t\t"
2389 : "c_bl : %g\n"
2390 : "cg_rate : %g\t\t"
2391 : "dir : %s\n"
2392 : "superscs : %s\t\t"
2393 : "eps : %g\n"
2394 : "(k0, k1, k2) : (%d, %d, %d)\t"
2395 : "ls : %d\n"
2396 : "max_iters : %d\t\t"
2397 : "max_time (min) : %g\n"
2398 : "memory : %d\t\t"
2399 : "normalize : %d\n"
2400 : "rho_x : %g\t\t"
2401 : "scale : %g\n"
2402 : "sigma : %g\t\t"
2403 : "sse : %g\n"
2404 : "thetabar : %g\t\t"
2405 : "warm_start : %d\n"
2406 : ".....................................................................\n\n",
2407 : (double) data->stgs->alpha,
2408 : (double) data->stgs->beta,
2409 : (double) data->stgs->c1,
2410 : (double) data->stgs->c_bl,
2411 : (double) data->stgs->cg_rate,
2412 : dir_string,
2413 1 : data->stgs->do_super_scs == 1 ? "yes" : "no",
2414 : (double) data->stgs->eps,
2415 : (int) data->stgs->k0,
2416 : (int) data->stgs->k1,
2417 : (int) data->stgs->k2,
2418 : (int) data->stgs->ls,
2419 : (int) data->stgs->max_iters,
2420 1 : (double) data->stgs->max_time_milliseconds / 60000,
2421 : (int) data->stgs->memory,
2422 : (int) data->stgs->normalize,
2423 : (double) data->stgs->rho_x,
2424 : (double) data->stgs->scale,
2425 : (double) data->stgs->sigma,
2426 : (double) data->stgs->sse,
2427 : (double) data->stgs->thetabar,
2428 : (int) data->stgs->warm_start);
2429 1 : }
2430 :
2431 2 : static void scs_print_allocated_memory(
2432 : const ScsData * RESTRICT data,
2433 : const ScsInfo * RESTRICT info) {
2434 2 : scs_int print_mode = data->stgs->do_override_streams;
2435 2 : FILE * stream = data->stgs->output_stream;
2436 2 : scs_float mem = 0;
2437 : char unit_mem[SCS_BYTE_UNIT_STR_LEN];
2438 2 : scs_units_allocated_mem(info, &mem, unit_mem);
2439 2 : scs_special_print(print_mode, stream, "Allocated Memory: %5.2f%s\n", (double) mem, unit_mem);
2440 2 : }
2441 :
2442 : /* this just calls scs_init, scs_solve, and scs_finish */
2443 341 : scs_int scs(
2444 2 : const ScsData * RESTRICT data,
2445 : const ScsCone * RESTRICT cone,
2446 : ScsSolution * RESTRICT sol,
2447 : ScsInfo * RESTRICT info) {
2448 : scs_int status;
2449 : /* --------------------------------------------------
2450 : * Create a Work object. It may happen that scs_init
2451 : * returns SCS_NULL (e.g., if some parameter is out
2452 : * of range, or memory could not be allocated).
2453 : * -------------------------------------------------- */
2454 341 : if (data->stgs->verbose >= 2) scs_print_parameter_details(data);
2455 341 : ScsWork *work = scs_init(data, cone, info);
2456 341 : scs_int print_mode = data->stgs->do_override_streams;
2457 341 : FILE * stream = data->stgs->output_stream;
2458 :
2459 341 : if (work != SCS_NULL) {
2460 329 : if (work->stgs->do_super_scs) {
2461 : /* solve with SuperSCS*/
2462 326 : if (work->stgs->verbose > 0) {
2463 2 : scs_special_print(print_mode, stream, "\nRunning SuperSCS...\n");
2464 2 : scs_compute_allocated_memory(work, cone, data, info);
2465 2 : scs_print_allocated_memory(data, info);
2466 : }
2467 326 : superscs_solve(work, data, cone, sol, info);
2468 : } else {
2469 : /* solve with SCS */
2470 3 : if (work->stgs->verbose > 0)
2471 0 : scs_special_print(print_mode, stream, "Running Standard SCS...\n");
2472 3 : scs_solve(work, data, cone, sol, info);
2473 : }
2474 329 : status = info->statusVal;
2475 : } else {
2476 12 : status = scs_failure(SCS_NULL, data ? data->m : -1, data ? data->n : -1, sol, info,
2477 : SCS_FAILED, "could not initialize work", "Failure", print_mode);
2478 : }
2479 341 : scs_finish(work);
2480 341 : return status;
2481 : }
2482 :
2483 19 : ScsSolution * scs_init_sol() {
2484 19 : ScsSolution * RESTRICT sol = scs_calloc(1, sizeof (* sol));
2485 19 : if (sol == SCS_NULL) {
2486 : /* LCOV_EXCL_START */
2487 : printf("ERROR: allocating sol failure\n");
2488 : return SCS_NULL;
2489 : /* LCOV_EXCL_STOP */
2490 : }
2491 19 : sol->s = SCS_NULL;
2492 19 : sol->x = SCS_NULL;
2493 19 : sol->y = SCS_NULL;
2494 19 : return sol;
2495 : }
2496 :
2497 19 : ScsInfo * scs_init_info() {
2498 19 : ScsInfo * RESTRICT info = scs_calloc(1, sizeof (*info));
2499 19 : if (info == SCS_NULL) {
2500 : /* LCOV_EXCL_START */
2501 : printf("ERROR: allocating info failure\n");
2502 : return SCS_NULL;
2503 : /* LCOV_EXCL_STOP */
2504 : }
2505 19 : info->pobj = NAN;
2506 19 : info->dobj = NAN;
2507 19 : info->iter = -1;
2508 19 : info->cg_total_iters = -1;
2509 19 : info->linsys_total_solve_time_ms = -1;
2510 19 : info->relGap = NAN;
2511 19 : info->resDual = NAN;
2512 19 : info->resInfeas = NAN;
2513 19 : info->resPri = NAN;
2514 19 : info->resUnbdd = NAN;
2515 19 : info->setupTime = NAN;
2516 19 : info->solveTime = NAN;
2517 19 : info->statusVal = SCS_INDETERMINATE;
2518 19 : info->progress_iter = SCS_NULL;
2519 19 : info->progress_dcost = SCS_NULL;
2520 19 : info->progress_pcost = SCS_NULL;
2521 19 : info->progress_relgap = SCS_NULL;
2522 19 : info->progress_resdual = SCS_NULL;
2523 19 : info->progress_respri = SCS_NULL;
2524 19 : info->progress_norm_fpr = SCS_NULL;
2525 19 : info->progress_time = SCS_NULL;
2526 19 : info->progress_mode = SCS_NULL;
2527 19 : info->progress_ls = SCS_NULL;
2528 19 : return info;
2529 : }
2530 :
2531 45 : ScsData * scs_init_data() {
2532 45 : ScsData * RESTRICT data = scs_malloc(sizeof (*data));
2533 :
2534 45 : if (data == SCS_NULL) {
2535 : /* LCOV_EXCL_START */
2536 : printf("ERROR: allocating data failure\n");
2537 : return SCS_NULL;
2538 : /* LCOV_EXCL_STOP */
2539 : }
2540 :
2541 45 : data->A = SCS_NULL;
2542 45 : data->b = SCS_NULL;
2543 45 : data->c = SCS_NULL;
2544 45 : data->m = 0;
2545 45 : data->n = 0;
2546 :
2547 45 : data->stgs = scs_malloc(sizeof (ScsSettings));
2548 45 : scs_set_default_settings(data);
2549 :
2550 45 : return data;
2551 : }
2552 :
|
