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 "scs_blas.h" /* contains BLAS(X) macros and type info */
30 :
31 : #define CONE_RATE (2)
32 : #define CONE_TOL (1e-8)
33 : #define CONE_THRESH (1e-6)
34 : #define EXP_CONE_MAX_ITERS (100)
35 : #define POW_CONE_MAX_ITERS (20)
36 :
37 : #ifdef LAPACK_LIB_FOUND
38 : extern void BLAS(syevr)(const char *jobz, const char *range, const char *uplo,
39 : blasint *n, scs_float *a, blasint *lda, scs_float *vl,
40 : scs_float *vu, blasint *il, blasint *iu, scs_float *abstol,
41 : blasint *m, scs_float *w, scs_float *z, blasint *ldz,
42 : blasint *isuppz, scs_float *work, blasint *lwork,
43 : blasint *iwork, blasint *liwork, blasint *info);
44 : extern void BLAS(syr)(const char *uplo, const blasint *n, const scs_float *alpha,
45 : const scs_float *x, const blasint *incx, scs_float *a,
46 : const blasint *lda);
47 : extern void BLAS(scal)(const blasint *n, const scs_float *sa, scs_float *sx,
48 : const blasint *incx);
49 : extern scs_float BLAS(nrm2)(const blasint *n, scs_float *x, const blasint *incx);
50 : #endif
51 :
52 : static scs_int getSdConeSize(scs_int s) {
53 660 : return (s * (s + 1)) / 2;
54 : }
55 :
56 : /*
57 : * boundaries will contain array of indices of rows of A corresponding to
58 : * cone boundaries, boundaries[0] is starting index for cones of size strictly
59 : * larger than 1
60 : * returns length of boundaries array, boundaries malloc-ed here so should be
61 : * freed
62 : */
63 326 : scs_int scs_get_cone_boundaries(
64 : const ScsCone * RESTRICT k,
65 : scs_int * * RESTRICT boundaries) {
66 326 : scs_int i, count = 0;
67 326 : scs_int len = 1 + k->qsize + k->ssize + k->ed + k->ep + k->psize;
68 326 : scs_int *RESTRICT b = scs_malloc(sizeof (scs_int) * len);
69 326 : b[count] = k->f + k->l;
70 326 : count += 1;
71 326 : if (k->qsize > 0) {
72 323 : memcpy(&b[count], k->q, k->qsize * sizeof (scs_int));
73 : }
74 326 : count += k->qsize;
75 329 : for (i = 0; i < k->ssize; ++i) {
76 6 : b[count + i] = getSdConeSize(k->s[i]);
77 : }
78 : count += k->ssize;
79 122 : for (i = 0; i < k->ep + k->ed; ++i) {
80 122 : b[count + i] = 3;
81 : }
82 326 : count += k->ep + k->ed;
83 330 : for (i = 0; i < k->psize; ++i) {
84 4 : b[count + i] = 3;
85 : }
86 : /* count += k->psize; */
87 326 : *boundaries = b;
88 326 : return len;
89 : }
90 :
91 341 : static scs_int getFullConeDims(const ScsCone *RESTRICT k) {
92 341 : scs_int i, c = 0;
93 341 : if (k->f)
94 0 : c += k->f;
95 341 : if (k->l)
96 7 : c += k->l;
97 341 : if (k->qsize && k->q) {
98 338 : for (i = 0; i < k->qsize; ++i) {
99 338 : c += k->q[i];
100 : }
101 : }
102 341 : if (k->ssize && k->s) {
103 3 : for (i = 0; i < k->ssize; ++i) {
104 6 : c += getSdConeSize(k->s[i]);
105 : }
106 : }
107 341 : if (k->ed)
108 1 : c += 3 * k->ed;
109 341 : if (k->ep)
110 2 : c += 3 * k->ep;
111 341 : if (k->p)
112 1 : c += 3 * k->psize;
113 341 : return c;
114 : }
115 :
116 341 : scs_int scs_validate_cones(const ScsData * RESTRICT d, const ScsCone * RESTRICT k) {
117 : scs_int i;
118 341 : if (getFullConeDims(k) != d->m) {
119 0 : scs_printf("cone dimensions %li not equal to num rows in A = m = %li\n",
120 : (long) getFullConeDims(k), (long) d->m);
121 0 : return -1;
122 : }
123 341 : if (k->f && k->f < 0) {
124 : scs_printf("free cone error\n");
125 0 : return -1;
126 : }
127 341 : if (k->l && k->l < 0) {
128 : scs_printf("lp cone error\n");
129 0 : return -1;
130 : }
131 341 : if (k->qsize && k->q) {
132 338 : if (k->qsize < 0) {
133 : scs_printf("soc cone error\n");
134 0 : return -1;
135 : }
136 338 : for (i = 0; i < k->qsize; ++i) {
137 338 : if (k->q[i] < 0) {
138 : scs_printf("soc cone error\n");
139 0 : return -1;
140 : }
141 : }
142 : }
143 341 : if (k->ssize && k->s) {
144 1 : if (k->ssize < 0) {
145 : scs_printf("sd cone error\n");
146 0 : return -1;
147 : }
148 3 : for (i = 0; i < k->ssize; ++i) {
149 3 : if (k->s[i] < 0) {
150 : scs_printf("sd cone error\n");
151 0 : return -1;
152 : }
153 : }
154 : }
155 341 : if (k->ed && k->ed < 0) {
156 : scs_printf("ep cone error\n");
157 0 : return -1;
158 : }
159 341 : if (k->ep && k->ep < 0) {
160 : scs_printf("ed cone error\n");
161 0 : return -1;
162 : }
163 341 : if (k->psize && k->p) {
164 1 : if (k->psize < 0) {
165 : scs_printf("power cone error\n");
166 0 : return -1;
167 : }
168 4 : for (i = 0; i < k->psize; ++i) {
169 4 : if (k->p[i] < -1 || k->p[i] > 1) {
170 : scs_printf("power cone error, values must be in [-1,1]\n");
171 0 : return -1;
172 : }
173 : }
174 : }
175 : return 0;
176 : }
177 :
178 2 : char *scs_get_cone_summary(const ScsInfo * RESTRICT info, ScsConeWork * RESTRICT c) {
179 2 : char *str = scs_malloc(sizeof (char) * 64);
180 2 : sprintf(str, "\tCones: avg projection time: %1.2es\n",
181 2 : c->total_cone_time / (info->iter + 1) / 1e3);
182 2 : c->total_cone_time = 0.0;
183 2 : return str;
184 : }
185 :
186 329 : void scs_finish_cone(ScsConeWork * RESTRICT c) {
187 : #ifdef LAPACK_LIB_FOUND
188 329 : scs_free(c->Xs);
189 329 : scs_free(c->Z);
190 329 : scs_free(c->e);
191 329 : scs_free(c->work);
192 329 : scs_free(c->iwork);
193 : #endif
194 329 : scs_free(c);
195 329 : }
196 :
197 2 : char *scs_get_cone_header(const ScsCone * RESTRICT k) {
198 2 : char *tmp = scs_malloc(sizeof (char) * 512);
199 : scs_int i, socVars, socBlks, sdVars, sdBlks;
200 : sprintf(tmp, "Cones:");
201 2 : if (k->f) {
202 0 : sprintf(tmp + strlen(tmp), "\tprimal zero / dual free vars: %li\n",
203 : (long) k->f);
204 : }
205 2 : if (k->l) {
206 0 : sprintf(tmp + strlen(tmp), "\tlinear vars: %li\n", (long) k->l);
207 : }
208 2 : socVars = 0;
209 2 : socBlks = 0;
210 2 : if (k->qsize && k->q) {
211 : socBlks = k->qsize;
212 2 : for (i = 0; i < k->qsize; i++) {
213 2 : socVars += k->q[i];
214 : }
215 2 : sprintf(tmp + strlen(tmp), "\tsoc vars: %li, soc blks: %li\n",
216 : (long) socVars, (long) socBlks);
217 : }
218 2 : sdVars = 0;
219 2 : sdBlks = 0;
220 2 : if (k->ssize && k->s) {
221 : sdBlks = k->ssize;
222 0 : for (i = 0; i < k->ssize; i++) {
223 0 : sdVars += getSdConeSize(k->s[i]);
224 : }
225 0 : sprintf(tmp + strlen(tmp), "\tsd vars: %li, sd blks: %li\n",
226 : (long) sdVars, (long) sdBlks);
227 : }
228 2 : if (k->ep || k->ed) {
229 0 : sprintf(tmp + strlen(tmp), "\texp vars: %li, dual exp vars: %li\n",
230 0 : (long) 3 * k->ep, (long) 3 * k->ed);
231 : }
232 2 : if (k->psize && k->p) {
233 0 : sprintf(tmp + strlen(tmp), "\tprimal + dual power vars: %li\n",
234 : (long) 3 * k->psize);
235 : }
236 2 : return tmp;
237 : }
238 :
239 : static scs_int isSimpleSemiDefiniteCone(scs_int * RESTRICT s, scs_int ssize) {
240 : scs_int i;
241 0 : for (i = 0; i < ssize; i++) {
242 1 : if (s[i] > 2) {
243 : return 0; /* false */
244 : }
245 : }
246 : return 1; /* true */
247 : }
248 :
249 50435 : static scs_float expNewtonOneD(scs_float rho, scs_float y_hat, scs_float z_hat) {
250 50435 : scs_float t = MAX(-z_hat, 1e-6);
251 : scs_float f, fp;
252 : scs_int i;
253 271525 : for (i = 0; i < EXP_CONE_MAX_ITERS; ++i) {
254 271525 : f = t * (t + z_hat) / rho / rho - y_hat / rho + log(t / rho) + 1;
255 271525 : fp = (2 * t + z_hat) / rho / rho + 1 / t;
256 :
257 271525 : t = t - f / fp;
258 :
259 271525 : if (t <= -z_hat) {
260 : return 0;
261 271525 : } else if (t <= 0) {
262 : return z_hat;
263 271525 : } else if (ABS(f) < CONE_TOL) {
264 : break;
265 : }
266 : }
267 50435 : return t + z_hat;
268 : }
269 :
270 50435 : static void expSolveForXWithRho(scs_float * RESTRICT v, scs_float * RESTRICT x, scs_float rho) {
271 50435 : x[2] = expNewtonOneD(rho, v[1], v[2]);
272 50435 : x[1] = (x[2] - v[2]) * x[2] / rho;
273 50435 : x[0] = v[0] - rho;
274 50435 : }
275 :
276 50435 : static scs_float expCalcGrad(scs_float * RESTRICT v, scs_float * RESTRICT x, scs_float rho) {
277 50435 : expSolveForXWithRho(v, x, rho);
278 50435 : if (x[1] <= 1e-12) {
279 0 : return x[0];
280 : }
281 50435 : return x[0] + x[1] * log(x[1] / x[2]);
282 : }
283 :
284 1761 : static void expGetRhoUb(scs_float * RESTRICT v, scs_float * RESTRICT x, scs_float * RESTRICT ub, scs_float *lb) {
285 1761 : *lb = 0;
286 1761 : *ub = 0.125;
287 7590 : while (expCalcGrad(v, x, *ub) > 0) {
288 4068 : *lb = *ub;
289 4068 : (*ub) *= 2;
290 : }
291 1761 : }
292 :
293 : /* project onto the exponential cone, v has dimension *exactly* 3 */
294 1893 : static scs_int projExpCone(scs_float * RESTRICT v, scs_int iter) {
295 : scs_int i;
296 : scs_float ub, lb, rho, g, x[3];
297 1893 : scs_float r = v[0], s = v[1], t = v[2];
298 1893 : scs_float tol = CONE_TOL; /* iter < 0 ? CONE_TOL : MAX(CONE_TOL, 1 /
299 : POWF((iter + 1), CONE_RATE)); */
300 :
301 : /* v in cl(Kexp) */
302 3689 : if ((s > 0 && s * exp(r / s) - t <= CONE_THRESH) ||
303 1796 : (r <= 0 && s == 0 && t >= 0)) {
304 : return 0;
305 : }
306 :
307 : /* -v in Kexp^* */
308 3573 : if ((-r < 0 && r * exp(s / r) + exp(1) * t <= CONE_THRESH) ||
309 1777 : (-r == 0 && -s >= 0 && -t >= 0)) {
310 : memset(v, 0, 3 * sizeof (scs_float));
311 : return 0;
312 : }
313 :
314 : /* special case with analytical solution */
315 1777 : if (r < 0 && s < 0) {
316 16 : v[1] = 0.0;
317 16 : v[2] = MAX(v[2], 0);
318 : return 0;
319 : }
320 :
321 : /* iterative procedure to find projection, bisects on dual variable: */
322 1761 : expGetRhoUb(v, x, &ub, &lb); /* get starting upper and lower bounds */
323 44606 : for (i = 0; i < EXP_CONE_MAX_ITERS; ++i) {
324 44606 : rho = (ub + lb) / 2; /* halfway between upper and lower bounds */
325 44606 : g = expCalcGrad(v, x, rho); /* calculates gradient wrt dual var */
326 44606 : if (g > 0) {
327 22517 : lb = rho;
328 : } else {
329 22089 : ub = rho;
330 : }
331 44606 : if (ub - lb < tol) {
332 : break;
333 : }
334 : }
335 1761 : v[0] = x[0];
336 1761 : v[1] = x[1];
337 1761 : v[2] = x[2];
338 : return 0;
339 : }
340 :
341 1 : static scs_int setUpSdScsConeWorkSpace(ScsConeWork * RESTRICT c, const ScsCone * RESTRICT k) {
342 : #ifdef LAPACK_LIB_FOUND
343 : scs_int i;
344 1 : blasint nMax = 0;
345 1 : scs_float eigTol = 1e-8;
346 1 : blasint negOne = -1;
347 1 : blasint m = 0;
348 : blasint info;
349 : scs_float wkopt;
350 :
351 : /* eigenvector decomp workspace */
352 4 : for (i = 0; i < k->ssize; ++i) {
353 3 : if (k->s[i] > nMax) {
354 1 : nMax = (blasint) k->s[i];
355 : }
356 : }
357 1 : c->Xs = scs_calloc(nMax * nMax, sizeof (scs_float));
358 1 : c->Z = scs_calloc(nMax * nMax, sizeof (scs_float));
359 1 : c->e = scs_calloc(nMax, sizeof (scs_float));
360 :
361 1 : BLAS(syevr)("Vectors", "All", "Lower", &nMax, c->Xs, &nMax, SCS_NULL,
362 : SCS_NULL, SCS_NULL, SCS_NULL, &eigTol, &m, c->e, c->Z, &nMax,
363 : SCS_NULL, &wkopt, &negOne, &(c->liwork), &negOne, &info);
364 :
365 1 : if (info != 0) {
366 0 : scs_printf("FATAL: syevr failure, info = %li\n", (long) info);
367 : return -1;
368 : }
369 1 : c->lwork = (blasint) (wkopt + 0.01); /* 0.01 for int casting safety */
370 1 : c->work = scs_malloc(c->lwork * sizeof (scs_float));
371 1 : c->iwork = scs_malloc(c->liwork * sizeof (blasint));
372 :
373 1 : if (c->Xs == SCS_NULL || c->Z == SCS_NULL
374 1 : || c->e == SCS_NULL || c->work == SCS_NULL
375 1 : || c->iwork == SCS_NULL) {
376 : return -1;
377 : }
378 : return 0;
379 : #else
380 : scs_printf("FATAL: Cannot solve SDPs with > 2x2 matrices without linked "
381 : "blas+lapack libraries\n");
382 : scs_printf("Install blas+lapack and re-compile SCS with blas+lapack libray "
383 : "locations\n");
384 : return -1;
385 : #endif
386 : }
387 :
388 329 : ScsConeWork *scs_init_conework(const ScsCone * RESTRICT k) {
389 329 : ScsConeWork * RESTRICT coneWork = scs_calloc(1, sizeof (ScsConeWork));
390 329 : coneWork->total_cone_time = 0.0;
391 329 : if (k->ssize && k->s) {
392 2 : if (isSimpleSemiDefiniteCone(k->s, k->ssize) == 0 &&
393 1 : setUpSdScsConeWorkSpace(coneWork, k) < 0) {
394 0 : scs_finish_cone(coneWork);
395 0 : return SCS_NULL;
396 : }
397 : }
398 329 : return coneWork;
399 : }
400 :
401 109 : scs_int project2By2Sdc(scs_float *X) {
402 : scs_float a, b, d, l1, l2, x1, x2, rad;
403 109 : scs_float sqrt2 = SQRTF(2.0);
404 109 : a = X[0];
405 109 : b = X[1] / sqrt2;
406 109 : d = X[2];
407 :
408 109 : if (ABS(b) < 1e-6) { /* diagonal matrix */
409 0 : X[0] = MAX(a, 0);
410 0 : X[1] = 0;
411 0 : X[2] = MAX(d, 0);
412 0 : return 0;
413 : }
414 :
415 109 : rad = SQRTF((a - d) * (a - d) + 4 * b * b);
416 : /* l1 >= l2 always, since rad >= 0 */
417 109 : l1 = 0.5 * (a + d + rad);
418 109 : l2 = 0.5 * (a + d - rad);
419 :
420 :
421 109 : if (l2 >= 0) { /* both eigs positive already */
422 : return 0;
423 : }
424 91 : if (l1 <= 0) { /* both eigs negative, set to 0 */
425 91 : X[0] = 0;
426 91 : X[1] = 0;
427 91 : X[2] = 0;
428 91 : return 0;
429 : }
430 :
431 : /* l1 pos, l2 neg */
432 0 : x1 = 1 / SQRTF(1 + (l1 - a) * (l1 - a) / b / b);
433 0 : x2 = x1 * (l1 - a) / b;
434 :
435 0 : X[0] = l1 * x1 * x1;
436 0 : X[1] = (l1 * x1 * x2) * sqrt2;
437 0 : X[2] = l1 * x2 * x2;
438 0 : return 0;
439 : }
440 :
441 : /* size of X is getSdConeSize(n) */
442 327 : static scs_int projSemiDefiniteCone(
443 : scs_float * RESTRICT X,
444 : const scs_int n,
445 : ScsConeWork * RESTRICT c,
446 : const scs_int iter) {
447 : /* project onto the positive semi-definite cone */
448 : #ifdef LAPACK_LIB_FOUND
449 : scs_int i;
450 327 : blasint one = 1;
451 327 : blasint m = 0;
452 327 : blasint nb = (blasint) n;
453 327 : blasint nbPlusOne = (blasint) (n + 1);
454 327 : blasint coneSz = (blasint) (getSdConeSize(n));
455 :
456 327 : scs_float sqrt2 = SQRTF(2.0);
457 327 : scs_float sqrt2Inv = 1.0 / sqrt2;
458 327 : scs_float *RESTRICT Xs = c->Xs;
459 327 : scs_float *RESTRICT Z = c->Z;
460 327 : scs_float *RESTRICT e = c->e;
461 327 : scs_float *RESTRICT work = c->work;
462 327 : blasint *RESTRICT iwork = c->iwork;
463 327 : blasint lwork = c->lwork;
464 327 : blasint liwork = c->liwork;
465 :
466 327 : scs_float eigTol = CONE_TOL; /* iter < 0 ? CONE_TOL : MAX(CONE_TOL, 1 /
467 : POWF(iter + 1, CONE_RATE)); */
468 327 : scs_float zero = 0.0;
469 : blasint info;
470 : scs_float vupper;
471 : #endif /* LAPACK_LIB_FOUND */
472 327 : if (n == 0) {
473 : return 0;
474 : }
475 327 : if (n == 1) {
476 0 : if (X[0] < 0.0) {
477 0 : X[0] = 0.0;
478 : }
479 : return 0;
480 : }
481 327 : if (n == 2) {
482 109 : return project2By2Sdc(X);
483 : }
484 : #ifdef LAPACK_LIB_FOUND
485 : /* expand lower triangular matrix to full matrix */
486 4360 : for (i = 0; i < n; ++i) {
487 4360 : memcpy(&(Xs[i * (n + 1)]), &(X[i * n - ((i - 1) * i) / 2]),
488 4360 : (n - i) * sizeof (scs_float));
489 : }
490 : /*
491 : rescale so projection works, and matrix norm preserved
492 : see http://www.seas.ucla.edu/~vandenbe/publications/mlbook.pdf pg 3
493 : */
494 : /* scale diags by sqrt(2) */
495 218 : BLAS(scal)(&nb, &sqrt2, Xs, &nbPlusOne); /* not nSquared */
496 :
497 : /* max-eig upper bounded by frobenius norm */
498 218 : vupper = 1.1 * sqrt2 *
499 218 : BLAS(nrm2)(&coneSz, X,
500 : &one); /* mult by factor to make sure is upper bound */
501 218 : vupper = MAX(vupper, 0.01);
502 :
503 : /* Solve eigenproblem, reuse workspaces */
504 218 : BLAS(syevr)("Vectors", "VInterval", "Lower", &nb, Xs, &nb, &zero, &vupper,
505 : SCS_NULL, SCS_NULL, &eigTol, &m, e, Z, &nb, SCS_NULL, work,
506 : &lwork, iwork, &liwork, &info);
507 :
508 218 : if (info < 0)
509 : return -1;
510 :
511 218 : memset(Xs, 0, n * n * sizeof (scs_float));
512 2800 : for (i = 0; i < m; ++i) {
513 2582 : scs_float a = e[i];
514 2582 : BLAS(syr)("Lower", &nb, &a, &(Z[i * n]), &one, Xs, &nb);
515 : }
516 : /* scale diags by 1/sqrt(2) */
517 218 : BLAS(scal)(&nb, &sqrt2Inv, Xs, &nbPlusOne); /* not nSquared */
518 : /* extract just lower triangular matrix */
519 4578 : for (i = 0; i < n; ++i) {
520 4360 : memcpy(&(X[i * n - ((i - 1) * i) / 2]), &(Xs[i * (n + 1)]),
521 4360 : (n - i) * sizeof (scs_float));
522 : }
523 :
524 : #else /* LAPACK_LIB_FOUND */
525 : scs_printf("FAILURE: solving SDP with > 2x2 matrices, but no blas/lapack "
526 : "libraries were linked!\n");
527 : scs_printf("SCS will return nonsense!\n");
528 : scs_scale_array(X, NAN, n);
529 : return -1;
530 : #endif /* LAPACK_LIB_FOUND */
531 : return 0;
532 : }
533 :
534 2790 : static scs_float powCalcX(scs_float r, scs_float xh, scs_float rh, scs_float a) {
535 2790 : scs_float x = 0.5 * (xh + SQRTF(xh * xh + 4 * a * (rh - r) * r));
536 2790 : return MAX(x, 1e-12);
537 : }
538 :
539 : static scs_float powCalcdxdr(scs_float x, scs_float xh, scs_float rh, scs_float r,
540 : scs_float a) {
541 2324 : return a * (rh - 2 * r) / (2 * x - xh);
542 : }
543 :
544 1395 : static scs_float powCalcF(scs_float x, scs_float y, scs_float r, scs_float a) {
545 1395 : return POWF(x, a) * POWF(y, (1 - a)) - r;
546 : }
547 :
548 1162 : static scs_float powCalcFp(scs_float x, scs_float y, scs_float dxdr, scs_float dydr,
549 : scs_float a) {
550 1162 : return POWF(x, a) * POWF(y, (1 - a)) * (a * dxdr / x + (1 - a) * dydr / y) -
551 : 1;
552 : }
553 :
554 368 : static void projPowerCone(scs_float *RESTRICT v, scs_float a) {
555 368 : scs_float xh = v[0], yh = v[1], rh = ABS(v[2]);
556 : scs_float x, y, r;
557 : scs_int i;
558 : /* v in K_a */
559 442 : if (xh >= 0 && yh >= 0 &&
560 74 : CONE_THRESH + POWF(xh, a) * POWF(yh, (1 - a)) >= rh)
561 : return;
562 :
563 : /* -v in K_a^* */
564 544 : if (xh <= 0 && yh <= 0 &&
565 199 : CONE_THRESH + POWF(-xh, a) * POWF(-yh, 1 - a) >=
566 199 : rh * POWF(a, a) * POWF(1 - a, 1 - a)) {
567 112 : v[0] = v[1] = v[2] = 0;
568 112 : return;
569 : }
570 :
571 233 : r = rh / 2;
572 1395 : for (i = 0; i < POW_CONE_MAX_ITERS; ++i) {
573 : scs_float f, fp, dxdr, dydr;
574 1395 : x = powCalcX(r, xh, rh, a);
575 1395 : y = powCalcX(r, yh, rh, 1 - a);
576 :
577 1395 : f = powCalcF(x, y, r, a);
578 1395 : if (ABS(f) < CONE_TOL)
579 : break;
580 :
581 1162 : dxdr = powCalcdxdr(x, xh, rh, r, a);
582 2324 : dydr = powCalcdxdr(y, yh, rh, r, (1 - a));
583 1162 : fp = powCalcFp(x, y, dxdr, dydr, a);
584 :
585 1162 : r = MAX(r - f / fp, 0);
586 1162 : r = MIN(r, rh);
587 : }
588 233 : v[0] = x;
589 233 : v[1] = y;
590 233 : v[2] = (v[2] < 0) ? -(r) : (r);
591 : }
592 :
593 : /* outward facing cone projection routine, iter is outer algorithm iteration, if
594 : iter < 0 then iter is ignored
595 : warm_start contains guess of projection (can be set to SCS_NULL) */
596 79662 : scs_int scs_project_dual_cone(
597 : scs_float * RESTRICT x,
598 : const ScsCone * RESTRICT k,
599 : ScsConeWork * RESTRICT c,
600 : const scs_float * RESTRICT warm_start,
601 : scs_int iter) {
602 : scs_int i;
603 79662 : scs_int count = (k->f ? k->f : 0);
604 : ScsTimer coneTimer;
605 79662 : scs_tic(&coneTimer);
606 :
607 79662 : if (k->l) {
608 : /* project onto positive orthant */
609 1144300 : for (i = count; i < count + k->l; ++i) {
610 1144300 : if (x[i] < 0.0)
611 1108700 : x[i] = 0.0;
612 : /* x[i] = (x[i] < 0.0) ? 0.0 : x[i]; */
613 : }
614 : count += k->l;
615 : }
616 :
617 79662 : if (k->qsize && k->q) {
618 : /* project onto SOC */
619 79596 : for (i = 0; i < k->qsize; ++i) {
620 79596 : if (k->q[i] == 0) {
621 0 : continue;
622 : }
623 79596 : if (k->q[i] == 1) {
624 0 : if (x[count] < 0.0)
625 0 : x[count] = 0.0;
626 : } else {
627 79596 : scs_float v1 = x[count];
628 79596 : scs_float s = scs_norm(&(x[count + 1]), k->q[i] - 1);
629 79596 : scs_float alpha = (s + v1) / 2.0;
630 :
631 79596 : if (s <= v1) { /* do nothing */
632 79253 : } else if (s <= -v1) {
633 2898 : memset(&(x[count]), 0, k->q[i] * sizeof (scs_float));
634 : } else {
635 76355 : x[count] = alpha;
636 76355 : scs_scale_array(&(x[count + 1]), alpha / s, k->q[i] - 1);
637 : }
638 : }
639 79596 : count += k->q[i];
640 : }
641 : }
642 :
643 79662 : if (k->ssize && k->s) {
644 : /* project onto PSD cone */
645 327 : for (i = 0; i < k->ssize; ++i) {
646 327 : if (k->s[i] == 0) {
647 0 : continue;
648 : }
649 327 : if (projSemiDefiniteCone(&(x[count]), k->s[i], c, iter) < 0)
650 : return -1;
651 654 : count += getSdConeSize(k->s[i]);
652 : }
653 : }
654 :
655 79662 : if (k->ep) {
656 : scs_float r, s, t;
657 : scs_int idx;
658 : /*
659 : * exponential cone is not self dual, if s \in K
660 : * then y \in K^* and so if K is the primal cone
661 : * here we project onto K^*, via Moreau
662 : * \Pi_C^*(y) = y + \Pi_C(-y)
663 : */
664 51 : scs_scale_array(&(x[count]), -1, 3 * k->ep); /* x = -x; */
665 : #ifdef _OPENMP
666 : #pragma omp parallel for private(r, s, t, idx)
667 : #endif
668 204 : for (i = 0; i < k->ep; ++i) {
669 153 : idx = count + 3 * i;
670 153 : r = x[idx];
671 153 : s = x[idx + 1];
672 153 : t = x[idx + 2];
673 :
674 153 : projExpCone(&(x[idx]), iter);
675 :
676 153 : x[idx] -= r;
677 153 : x[idx + 1] -= s;
678 153 : x[idx + 2] -= t;
679 : }
680 51 : count += 3 * k->ep;
681 : }
682 :
683 79662 : if (k->ed) {
684 : /* exponential cone: */
685 : #ifdef _OPENMP
686 : #pragma omp parallel for
687 : #endif
688 1740 : for (i = 0; i < k->ed; ++i) {
689 1740 : projExpCone(&(x[count + 3 * i]), iter);
690 : }
691 15 : count += 3 * k->ed;
692 : }
693 :
694 79662 : if (k->psize && k->p) {
695 : scs_float v[3];
696 : scs_int idx;
697 : /* don't use openmp for power cone
698 : ifdef _OPENMP
699 : pragma omp parallel for private(v, idx)
700 : endif
701 : */
702 368 : for (i = 0; i < k->psize; ++i) {
703 368 : idx = count + 3 * i;
704 368 : if (k->p[i] <= 0) {
705 : /* dual power cone */
706 184 : projPowerCone(&(x[idx]), -k->p[i]);
707 : } else {
708 : /* primal power cone, using Moreau */
709 184 : v[0] = -x[idx];
710 184 : v[1] = -x[idx + 1];
711 184 : v[2] = -x[idx + 2];
712 :
713 184 : projPowerCone(v, k->p[i]);
714 :
715 184 : x[idx] += v[0];
716 184 : x[idx + 1] += v[1];
717 184 : x[idx + 2] += v[2];
718 : }
719 : }
720 : /* count += 3 * k->psize; */
721 : }
722 : /* project onto OTHER cones */
723 79662 : if (c) {
724 79662 : c->total_cone_time += scs_toc_quiet(&coneTimer);
725 : }
726 : return 0;
727 : }
|
