1 /*
2 * ECOS - Embedded Conic Solver.
3 * Copyright (C) 2012-2015 A. Domahidi [domahidi@embotech.com],
4 * Automatic Control Lab, ETH Zurich & embotech GmbH, Zurich, Switzerland.
5 *
6 * This program is free software: you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation, either version 3 of the License, or
9 * (at your option) any later version.
10 *
11 * This program is distributed in the hope that it will be useful,
12 * but WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
14 * GNU General Public License for more details.
15 *
16 * You should have received a copy of the GNU General Public License
17 * along with this program. If not, see <http://www.gnu.org/licenses/>.
18 */
19
20 /* Main solver module */
21
22
23 /* ECOS HEADER FILE ---------------------------------------------------- */
24 #include "ecos.h"
25 #include "splamm.h"
26 #include "equil.h"
27
28 /* NEEDED FOR SQRT ----------------------------------------------------- */
29 #include <math.h>
30
31 #if (defined _WIN32 || defined _WIN64 )
32 /* include stdio.h for _set_output_format */
33 #include <stdio.h>
34 #endif
35
36
37
38 /* Some internal defines */
39 #define ECOS_NOT_CONVERGED_YET (-87) /* indicates no convergence yet */
40
41 /**
42 * Private static const char * for version numbering.
43 * All versions point to this string.
44 */
45 static const char* thisVersion = ECOS_VERSION;
46
47
48 /**
49 * Version: returns the current version number
50 * Returns the version number in the format X.Y.Z
51 */
ECOS_ver(void)52 const char* ECOS_ver(void)
53 {
54 return thisVersion;
55 }
56
57
58 /* Compares stats of two iterates with each other.
59 * Returns 1 if infoA is better than infoB, zero otherwise.
60 */
compareStatistics(stats * infoA,stats * infoB)61 idxint compareStatistics(stats* infoA, stats* infoB)
62 {
63
64 if ( infoA->pinfres != ECOS_NAN && infoA->kapovert > 1){
65 if( infoB->pinfres != ECOS_NAN) {
66 /* A->pinfres != NAN, B->pinfres!=NAN */
67 if ( ( infoA->gap > 0 && infoB->gap > 0 && infoA->gap < infoB->gap ) &&
68 ( infoA->pinfres > 0 && infoA->pinfres < infoB->pres ) &&
69 ( infoA->mu > 0 && infoA->mu < infoB->mu ) ){
70 /* PRINTTEXT("BRANCH 1 "); */
71 return 1;
72 } else {
73 /* PRINTTEXT("BRANCH 1 not OK"); */
74 return 0;
75 }
76 } else {
77 /* A->pinfres != NAN, B->pinfres==NAN */
78 if ( ( infoA->gap > 0 && infoB->gap > 0 && infoA->gap < infoB->gap ) &&
79 ( infoA->mu > 0 && infoA->mu < infoB->mu ) ){
80 /* PRINTTEXT("BRANCH 2 "); */
81 return 1;
82 } else {
83 /* PRINTTEXT("BRANCH 2 not OK"); */
84 return 0;
85 }
86 }
87 } else {
88 /* A->pinfres == NAN or pinfres too large */
89 if ( ( infoA->gap > 0 && infoB->gap > 0 && infoA->gap < infoB->gap ) &&
90 ( infoA->pres > 0 && infoA->pres < infoB->pres ) &&
91 ( infoA->dres > 0 && infoA->dres < infoB->dres ) &&
92 ( infoA->kapovert > 0 && infoA->kapovert < infoB->kapovert) &&
93 ( infoA->mu > 0 && infoA->mu < infoB->mu ) ){
94 /* PRINTTEXT("BRANCH 3 OK "); */
95 return 1;
96 } else {
97 /* PRINTTEXT("BRANCH 3 not OK"); */
98 return 0;
99 }
100 }
101 }
102
103
104 /* Copy variables from current to best iterate */
saveIterateAsBest(pwork * w)105 void saveIterateAsBest(pwork* w)
106 {
107 idxint i;
108 for (i=0; i<w->n; i++) { w->best_x[i] = w->x[i]; }
109 for (i=0; i<w->p; i++) { w->best_y[i] = w->y[i]; }
110 for (i=0; i<w->m; i++) { w->best_z[i] = w->z[i]; }
111 for (i=0; i<w->m; i++) { w->best_s[i] = w->s[i]; }
112 w->best_kap = w->kap;
113 w->best_tau = w->tau;
114 w->best_cx = w->cx;
115 w->best_by = w->by;
116 w->best_hz = w->hz;
117 w->best_info->pcost = w->info->pcost;
118 w->best_info->dcost = w->info->dcost;
119 w->best_info->pres = w->info->pres;
120 w->best_info->dres = w->info->dres;
121 w->best_info->pinfres = w->info->pinfres;
122 w->best_info->dinfres = w->info->dinfres;
123 w->best_info->gap = w->info->gap;
124 w->best_info->relgap = w->info->relgap;
125 w->best_info->mu = w->info->mu;
126 w->best_info->kapovert = w->info->kapovert;
127 w->best_info->iter = w->info->iter;
128 }
129
130
131 /* Copy variables from best iterate to current
132 * Comapred to saveIterateAsBest, only the iteration number
133 * is not restored.
134 */
restoreBestIterate(pwork * w)135 void restoreBestIterate(pwork* w)
136 {
137 idxint i;
138 for (i=0; i<w->n; i++) { w->x[i] = w->best_x[i]; }
139 for (i=0; i<w->p; i++) { w->y[i] = w->best_y[i]; }
140 for (i=0; i<w->m; i++) { w->z[i] = w->best_z[i]; }
141 for (i=0; i<w->m; i++) { w->s[i] = w->best_s[i]; }
142 w->kap = w->best_kap;
143 w->tau = w->best_tau;
144 w->cx = w->best_cx;
145 w->by = w->best_by;
146 w->hz = w->best_hz;
147 w->info->pcost = w->best_info->pcost;
148 w->info->dcost = w->best_info->dcost;
149 w->info->pres = w->best_info->pres;
150 w->info->dres = w->best_info->dres;
151 w->info->pinfres = w->best_info->pinfres;
152 w->info->dinfres = w->best_info->dinfres;
153 w->info->gap = w->best_info->gap;
154 w->info->relgap = w->best_info->relgap;
155 w->info->mu = w->best_info->mu;
156 w->info->kapovert = w->best_info->kapovert;
157 }
158
159
160 /*
161 * This function is reponsible for checking the exit/convergence conditions of ECOS.
162 * If one of the exit conditions is met, ECOS displays an exit message and returns
163 * the corresponding exit code. The calling function must then make sure that ECOS
164 * is indeed correctly exited, so a call to this function should always be followed
165 * by a break statement.
166 *
167 * If mode == 0, normal precisions are checked.
168 *
169 * If mode != 0, reduced precisions are checked, and the exit display is augmented
170 * by "Close to". The exitcodes returned are increased by the value
171 * of mode.
172 *
173 * The primal and dual infeasibility flags w->info->pinf and w->info->dinf are raised
174 * according to the outcome of the test.
175 *
176 * If none of the exit tests are met, the function returns ECOS_NOT_CONVERGED_YET.
177 * This should not be an exitflag that is ever returned to the outside world.
178 */
checkExitConditions(pwork * w,idxint mode)179 idxint checkExitConditions(pwork* w, idxint mode)
180 {
181 pfloat feastol;
182 pfloat abstol;
183 pfloat reltol;
184
185 /* set accuracy against which to check */
186 if( mode == 0) {
187 /* check convergence against normal precisions */
188 feastol = w->stgs->feastol;
189 abstol = w->stgs->abstol;
190 reltol = w->stgs->reltol;
191 } else {
192 /* check convergence against reduced precisions */
193 feastol = w->stgs->feastol_inacc;
194 abstol = w->stgs->abstol_inacc;
195 reltol = w->stgs->reltol_inacc;
196 }
197
198 /* Optimal? */
199 if( ( -w->cx > 0 || -w->by - w->hz >= -abstol ) &&
200 ( w->info->pres < feastol && w->info->dres < feastol ) &&
201 ( w->info->gap < abstol || w->info->relgap < reltol )){
202 #if PRINTLEVEL > 0
203 if( w->stgs->verbose ) {
204 if( mode == 0) {
205 PRINTTEXT("\nOPTIMAL (within feastol=%3.1e, reltol=%3.1e, abstol=%3.1e).", MAX(w->info->dres, w->info->pres), w->info->relgap, w->info->gap);
206 } else {
207 PRINTTEXT("\nClose to OPTIMAL (within feastol=%3.1e, reltol=%3.1e, abstol=%3.1e).", MAX(w->info->dres, w->info->pres), w->info->relgap, w->info->gap);
208 }
209 }
210 #endif
211 w->info->pinf = 0;
212 w->info->dinf = 0;
213 return ECOS_OPTIMAL + mode;
214 }
215
216 /* Dual infeasible? */
217 else if( (w->info->dinfres != ECOS_NAN) && (w->info->dinfres < feastol) && (w->tau < w->kap) ){
218 #if PRINTLEVEL > 0
219 if( w->stgs->verbose ) {
220 if( mode == 0) {
221 PRINTTEXT("\nUNBOUNDED (within feastol=%3.1e).", w->info->dinfres );
222 } else {
223 PRINTTEXT("\nClose to UNBOUNDED (within feastol=%3.1e).", w->info->dinfres );
224 }
225 }
226 #endif
227 w->info->pinf = 0;
228 w->info->dinf = 1;
229 return ECOS_DINF + mode;
230 }
231
232 /* Primal infeasible? */
233 else if( ((w->info->pinfres != ECOS_NAN && w->info->pinfres < feastol) && (w->tau < w->kap)) ||
234 ( w->tau < w->stgs->feastol && w->kap < w->stgs->feastol && w->info->pinfres < w->stgs->feastol) ){
235 #if PRINTLEVEL > 0
236 if( w->stgs->verbose ) {
237 if( mode == 0) {
238 PRINTTEXT("\nPRIMAL INFEASIBLE (within feastol=%3.1e).", w->info->pinfres );
239 } else {
240 PRINTTEXT("\nClose to PRIMAL INFEASIBLE (within feastol=%3.1e).", w->info->pinfres );
241 }
242 }
243 #endif
244 w->info->pinf = 1;
245 w->info->dinf = 0;
246 return ECOS_PINF + mode;
247 }
248
249
250 /* Indicate if none of the above criteria are met */
251 else {
252 return ECOS_NOT_CONVERGED_YET;
253 }
254 }
255
256
257 /*
258 * Initializes the solver.
259 */
init(pwork * w)260 idxint init(pwork* w)
261 {
262 idxint i, j, k, l, KKT_FACTOR_RETURN_CODE;
263 idxint* Pinv = w->KKT->Pinv;
264 pfloat rx, ry, rz;
265
266 #if PROFILING > 1
267 timer tfactor, tkktsolve;
268 #endif
269
270 /* set regularization parameter */
271 w->KKT->delta = w->stgs->delta;
272
273 /* Initialize KKT matrix */
274 kkt_init(w->KKT->PKPt, w->KKT->PK, w->C);
275
276 #if DEBUG > 0
277 dumpSparseMatrix(w->KKT->PKPt, "PKPt0.txt");
278 #endif
279
280
281 /* initialize RHS1 */
282 k = 0; j = 0;
283 for( i=0; i<w->n; i++ ){ w->KKT->RHS1[w->KKT->Pinv[k++]] = 0; }
284 for( i=0; i<w->p; i++ ){ w->KKT->RHS1[w->KKT->Pinv[k++]] = w->b[i]; }
285 for( i=0; i<w->C->lpc->p; i++ ){ w->KKT->RHS1[w->KKT->Pinv[k++]] = w->h[i]; j++; }
286 for( l=0; l<w->C->nsoc; l++ ){
287 for( i=0; i < w->C->soc[l].p; i++ ){ w->KKT->RHS1[w->KKT->Pinv[k++]] = w->h[j++]; }
288 #if CONEMODE == 0
289 w->KKT->RHS1[w->KKT->Pinv[k++]] = 0;
290 w->KKT->RHS1[w->KKT->Pinv[k++]] = 0;
291 #endif
292 }
293 #ifdef EXPCONE
294 for( l=0; l<w->C->nexc; l++)
295 {
296 for(i=0; i<3; i++){ w->KKT->RHS1[w->KKT->Pinv[k++]] = w->h[j++]; }
297 }
298 #endif
299 #if PRINTLEVEL > 2
300 PRINTTEXT("Written %i entries of RHS1\n", (int)k);
301 #endif
302
303 /* initialize RHS2 */
304 for( i=0; i<w->n; i++ ){ w->KKT->RHS2[w->KKT->Pinv[i]] = -w->c[i]; }
305 for( i=w->n; i<w->KKT->PKPt->n; i++ ){ w->KKT->RHS2[w->KKT->Pinv[i]] = 0; }
306
307 /* get scalings of problem data */
308 rx = norm2(w->c, w->n); w->resx0 = MAX(1, rx);
309 ry = norm2(w->b, w->p); w->resy0 = MAX(1, ry);
310 rz = norm2(w->h, w->m); w->resz0 = MAX(1, rz);
311
312 #ifdef EXPCONE
313 /* If expcone is enabled and there are exponential cones in
314 * the problem definition, then skip the Mehrotra initialization strategy
315 */
316 if(w->C->nexc == 0)
317 {
318 #endif
319
320 /* Factor KKT matrix - this is needed in all 3 linear system solves */
321 #if PROFILING > 1
322 tic(&tfactor);
323 KKT_FACTOR_RETURN_CODE = kkt_factor(w->KKT, w->stgs->eps, w->stgs->delta, &w->info->tfactor_t1, &w->info->tfactor_t2);
324 w->info->tfactor += toc(&tfactor);
325 #else
326 KKT_FACTOR_RETURN_CODE = kkt_factor(w->KKT, w->stgs->eps, w->stgs->delta);
327 #endif
328
329 /* check if factorization was successful, exit otherwise */
330 if( KKT_FACTOR_RETURN_CODE != KKT_OK ){
331 #if PRINTLEVEL > 0
332 if( w->stgs->verbose ) PRINTTEXT("\nProblem in factoring KKT system, aborting.");
333 #endif
334 return ECOS_FATAL;
335 }
336
337
338 /*
339 * PRIMAL VARIABLES:
340 * - solve xhat = arg min ||Gx-h||_2^2 such that Ax = b
341 * - r = h - G*xhat
342 * These two equations are solved by
343 *
344 * [ 0 A' G' ] [ xhat ] [ 0 ]
345 * [ A 0 0 ] [ y ] = [ b ]
346 * [ G 0 -I ] [ -r ] [ h ]
347 *
348 * and then take shat = r if alphap < 0, zbar + (1+alphap)*e otherwise
349 * where alphap = inf{ alpha | sbar + alpha*e >= 0 }
350 */
351
352 /* Solve for RHS [0; b; h] */
353 #if PROFILING > 1
354 tic(&tkktsolve);
355 #endif
356 w->info->nitref1 = kkt_solve(w->KKT, w->A, w->G, w->KKT->RHS1, w->KKT->dx1, w->KKT->dy1, w->KKT->dz1, w->n, w->p, w->m, w->C, 1, w->stgs->nitref);
357 #if PROFILING > 1
358 w->info->tkktsolve += toc(&tkktsolve);
359 #endif
360
361 #if DEBUG > 0
362 #if PRINTLEVEL > 3
363 printDenseMatrix(w->KKT->dx1, w->n, 1, "dx1_init");
364 printDenseMatrix(w->KKT->dy1, w->p, 1, "dy1_init");
365 printDenseMatrix(w->KKT->dz1, w->m, 1, "dz1_init");
366 #endif
367 #endif
368
369 /* Copy out initial value of x */
370 for( i=0; i<w->n; i++ ){ w->x[i] = w->KKT->dx1[i]; }
371
372 /* Copy out -r into temporary variable */
373 for( i=0; i<w->m; i++ ){ w->KKT->work1[i] = -w->KKT->dz1[i]; }
374
375 /* Bring variable to cone */
376 bring2cone(w->C, w->KKT->work1, w->s );
377
378 /*
379 * dual variables
380 * solve (yhat,zbar) = arg min ||z||_2^2 such that G'*z + A'*y + c = 0
381 *
382 * we can solve this by
383 *
384 * [ 0 A' G' ] [ x ] [ -c ]
385 * [ A 0 0 ] [ yhat ] = [ 0 ]
386 * [ G 0 -I ] [ zbar ] [ 0 ]
387 *
388 * and then take zhat = zbar if alphad < 0, zbar + (1+alphad)*e otherwise
389 * where alphad = inf{ alpha | zbar + alpha*e >= 0 }
390 */
391
392 /* Solve for RHS [-c; 0; 0] */
393 #if PROFILING > 1
394 tic(&tkktsolve);
395 #endif
396 w->info->nitref2 = kkt_solve(w->KKT, w->A, w->G, w->KKT->RHS2, w->KKT->dx2, w->KKT->dy2, w->KKT->dz2, w->n, w->p, w->m, w->C, 1, w->stgs->nitref);
397 #if PROFILING > 1
398 w->info->tkktsolve += toc(&tkktsolve);
399 #endif
400
401 #if DEBUG > 0
402 #if PRINTLEVEL > 3
403 printDenseMatrix(w->KKT->dx2, w->n, 1, "dx2_init");
404 printDenseMatrix(w->KKT->dy2, w->p, 1, "dy2_init");
405 printDenseMatrix(w->KKT->dz2, w->m, 1, "dz2_init");
406 #endif
407 #endif
408
409 /* Copy out initial value of y */
410 for( i=0; i<w->p; i++ ){ w->y[i] = w->KKT->dy2[i]; }
411
412 /* Bring variable to cone */
413 bring2cone(w->C, w->KKT->dz2, w->z );
414
415 #ifdef EXPCONE
416 }
417 else if(w->C->nexc > 0)
418 {
419 unitInitialization(w->C, w->s, w->z, 1.0);
420
421 for( i=0; i<w->p; i++ ){ w->y[i] = 0.0; }
422 for( i=0; i<w->n; i++ ){ w->x[i] = 0.0; }
423 w->info->nitref1 = 0;
424 w->info->nitref2 = 0;
425 }
426 #endif
427
428 /* Prepare RHS1 - before this line RHS1 = [0; b; h], after it holds [-c; b; h] */
429 for( i=0; i<w->n; i++){ w->KKT->RHS1[Pinv[i]] = -w->c[i]; }
430
431 /*
432 * other variables
433 */
434 w->kap = 1.0;
435 w->tau = 1.0;
436
437 w->info->step = 0;
438 w->info->step_aff = 0;
439 w->info->dinf = 0;
440 w->info->pinf = 0;
441
442 return 0;
443 }
444
445
446
447 /*
448 * Computes residuals.
449 *
450 * hrx = -A'*y - G'*z; rx = hrx - c.*tau; hresx = norm(rx,2);
451 * hry = A*x; ry = hry - b.*tau; hresy = norm(ry,2);
452 * hrz = s + G*x; rz = hrz - h.*tau; hresz = norm(rz,2);
453 * rt = kappa + c'*x + b'*y + h'*z;
454 */
computeResiduals(pwork * w)455 void computeResiduals(pwork *w)
456 {
457 /* rx = -A'*y - G'*z - c.*tau */
458 if( w->p > 0 ) {
459 sparseMtVm(w->A, w->y, w->rx, 1, 0);
460 sparseMtVm(w->G, w->z, w->rx, 0, 0);
461 } else {
462 sparseMtVm(w->G, w->z, w->rx, 1, 0);
463 }
464 w->hresx = norm2(w->rx, w->n);
465 vsubscale(w->n, w->tau, w->c, w->rx);
466
467 /* ry = A*x - b.*tau */
468 if( w->p > 0 ){
469 sparseMV(w->A, w->x, w->ry, 1, 1);
470 w->hresy = norm2(w->ry, w->p);
471 vsubscale(w->p, w->tau, w->b, w->ry);
472 } else {
473 w->hresy = 0;
474 w->ry = NULL;
475 }
476
477 /* rz = s + G*x - h.*tau */
478 sparseMV(w->G, w->x, w->rz, 1, 1);
479 vadd(w->m, w->s, w->rz);
480 w->hresz = norm2(w->rz, w->m);
481 vsubscale(w->m, w->tau, w->h, w->rz);
482
483 /* rt = kappa + c'*x + b'*y + h'*z; */
484 w->cx = eddot(w->n, w->c, w->x);
485 w->by = w->p > 0 ? eddot(w->p, w->b, w->y) : 0.0;
486 w->hz = eddot(w->m, w->h, w->z);
487 w->rt = w->kap + w->cx + w->by + w->hz;
488
489 /* Norms of x y z */
490 w->nx = norm2(w->x,w->n);
491 w->ny = norm2(w->y,w->p);
492 w->ns = norm2(w->s,w->m);
493 w->nz = norm2(w->z,w->m);
494
495 }
496
497
498
499 /*
500 * Updates statistics.
501 */
updateStatistics(pwork * w)502 void updateStatistics(pwork* w)
503 {
504 pfloat nry, nrz;
505
506 stats* info = w->info;
507
508 /* mu = (s'*z + kap*tau) / (D+1) where s'*z is the duality gap */
509 info->gap = eddot(w->m, w->s, w->z);
510 info->mu = (info->gap + w->kap*w->tau) / (w->D + 1);
511
512 info->kapovert = w->kap / w->tau;
513 info->pcost = w->cx / w->tau;
514 info->dcost = -(w->hz + w->by) / w->tau;
515
516 /* relative duality gap */
517 if( info->pcost < 0 ){ info->relgap = info->gap / (-info->pcost); }
518 else if( info->dcost > 0 ){ info->relgap = info->gap / info->dcost; }
519 else info->relgap = ECOS_NAN;
520
521 /* residuals */
522 nry = w->p > 0 ? norm2(w->ry, w->p)/MAX(w->resy0+w->nx,1) : 0.0;
523 nrz = norm2(w->rz, w->m)/MAX(w->resz0+w->nx+w->ns,1);
524 info->pres = MAX(nry, nrz) / w->tau;
525 info->dres = norm2(w->rx, w->n)/MAX(w->resx0+w->ny+w->nz,1) / w->tau;
526
527 /* infeasibility measures
528 *
529 * CVXOPT uses the following:
530 * info->pinfres = w->hz + w->by < 0 ? w->hresx / w->resx0 / (-w->hz - w->by) : NAN;
531 * info->dinfres = w->cx < 0 ? MAX(w->hresy/w->resy0, w->hresz/w->resz0) / (-w->cx) : NAN;
532 */
533 info->pinfres = (w->hz + w->by)/MAX(w->ny+w->nz,1) < -w->stgs->reltol ? w->hresx / MAX(w->ny+w->nz,1) : ECOS_NAN;
534 info->dinfres = w->cx/MAX(w->nx,1) < -w->stgs->reltol ? MAX(w->hresy/MAX(w->nx,1), w->hresz/MAX(w->nx+w->ns,1)) : ECOS_NAN;
535
536
537
538 #if PRINTLEVEL > 2
539 PRINTTEXT("TAU=%6.4e KAP=%6.4e PINFRES=%6.4e DINFRES=%6.4e\n",w->tau,w->kap,info->pinfres, info->dinfres );
540 #endif
541
542 }
543
544
545
546 #if PRINTLEVEL > 0
printProgress(stats * info)547 void printProgress(stats* info)
548 {
549 if( info->iter == 0 )
550 {
551 /* print header at very first iteration */
552 #if PRINTLEVEL == 2
553 PRINTTEXT("\nECOS %s - (C) embotech GmbH, Zurich Switzerland, 2012-15. Web: www.embotech.com/ECOS\n", ECOS_VERSION);
554 #if PRINTLEVEL > 2
555 PRINTTEXT("Contributors: Core interior point solver: A. Domahidi, embotech GmbH\n");
556 PRINTTEXT(" Python interface and equilibration: Eric Chu, Stanford University\n");
557 PRINTTEXT(" Branch-and-bound module: Han Wang, Stanford University\n");
558 PRINTTEXT(" Maths and Methods: Stephen Boyd, Stanford University\n");
559 #ifdef EXPCONE
560 PRINTTEXT(" Exponential cone support: Santiago Akle, Stanford University\n");
561 #endif
562 #endif
563 PRINTTEXT("\n");
564 #endif
565 #if defined _WIN32 || defined _WIN64
566 #if defined EXPCONE
567 PRINTTEXT("It pcost dcost gap pres dres k/t mu step sigma IR | BT\n");
568 PRINTTEXT("%2d %+5.3e %+5.3e %+2.0e %2.0e %2.0e %2.0e %2.0e --- --- %2d %2d - | - - \n",(int)info->iter, info->pcost, info->dcost, info->gap, info->pres, info->dres, info->kapovert, info->mu, (int)info->nitref1, (int)info->nitref2);
569 #else
570 PRINTTEXT("It pcost dcost gap pres dres k/t mu step sigma IR\n");
571 PRINTTEXT("%2d %+5.3e %+5.3e %+2.0e %2.0e %2.0e %2.0e %2.0e --- --- %2d %2d -\n",(int)info->iter, info->pcost, info->dcost, info->gap, info->pres, info->dres, info->kapovert, info->mu, (int)info->nitref1, (int)info->nitref2);
572 #endif //if defined expcone
573 #else
574 #if defined EXPCONE
575 PRINTTEXT("It pcost dcost gap pres dres k/t mu step sigma IR | BT\n");
576 PRINTTEXT("%2d %+5.3e %+5.3e %+2.0e %2.0e %2.0e %2.0e %2.0e --- --- %2d %2d - | - - \n",(int)info->iter, info->pcost, info->dcost, info->gap, info->pres, info->dres, info->kapovert, info->mu, (int)info->nitref1, (int)info->nitref2);
577 #else
578 PRINTTEXT("It pcost dcost gap pres dres k/t mu step sigma IR\n");
579 PRINTTEXT("%2d %+5.3e %+5.3e %+2.0e %2.0e %2.0e %2.0e %2.0e --- --- %2d %2d -\n",(int)info->iter, info->pcost, info->dcost, info->gap, info->pres, info->dres, info->kapovert, info->mu, (int)info->nitref1, (int)info->nitref2);
580 #endif /*End if defined EXPCONE for the branch of the first iteration non windows */
581 #endif
582 } else {
583 #if defined _WIN32 || defined _WIN64
584 #if defined EXPCONE
585 PRINTTEXT("%2d %+5.3e %+5.3e %+2.0e %2.0e %2.0e %2.0e %2.0e %6.4f %2.0e %2d %2d %2d | %2d %2d\n",(int)info->iter, info->pcost, info->dcost, info->gap, info->pres, info->dres, info->kapovert, info->mu, info->step, info->sigma,\
586 (int)info->nitref1,\
587 (int)info->nitref2,\
588 (int)info->nitref3,\
589 (int)info->affBack,\
590 (int)info->cmbBack);
591 #else
592 PRINTTEXT("%2d %+5.3e %+5.3e %+2.0e %2.0e %2.0e %2.0e %2.0e %6.4f %2.0e %2d %2d %2d\n",(int)info->iter, info->pcost, info->dcost, info->gap, info->pres, info->dres, info->kapovert, info->mu, info->step, info->sigma,\
593 (int)info->nitref1,\
594 (int)info->nitref2,\
595 (int)info->nitref3);
596 #endif //End of if defined EXPCONE
597 #else
598 #if defined EXPCONE
599 PRINTTEXT("%2d %+5.3e %+5.3e %+2.0e %2.0e %2.0e %2.0e %2.0e %6.4f %2.0e %2d %2d %2d | %2d %2d\n",(int)info->iter, info->pcost, info->dcost, info->gap, info->pres, info->dres, info->kapovert, info->mu, info->step, info->sigma,\
600 (int)info->nitref1,\
601 (int)info->nitref2,\
602 (int)info->nitref3,\
603 (int)info->affBack,\
604 (int)info->cmbBack);
605 #else
606 PRINTTEXT("%2d %+5.3e %+5.3e %+2.0e %2.0e %2.0e %2.0e %2.0e %6.4f %2.0e %2d %2d %2d\n",(int)info->iter, info->pcost, info->dcost, info->gap, info->pres, info->dres, info->kapovert, info->mu, info->step, info->sigma,\
607 (int)info->nitref1,\
608 (int)info->nitref2,\
609 (int)info->nitref3);
610 #endif /* End of if defined expcone */
611 #endif /* end of windows block */
612 }
613
614 /* enable to flush printf in Matlab immediately */
615 #ifdef MATLAB_MEX_FILE
616 #if defined MATLAB_FLUSH_PRINTS
617 mexEvalString("pause(0.0001);");
618 #endif
619 #endif
620
621
622 }
623
deleteLastProgressLine(stats * info)624 void deleteLastProgressLine( stats* info )
625 {
626 idxint i;
627 idxint offset = 0;
628
629 if( info->kapovert < 0 ) offset++;
630 if( info->mu < 0) offset++;
631 if( info->pres < 0 ) offset++;
632 if (info->dres < 0 ) offset++;
633
634 for (i=0; i<82+offset; i++) {
635 PRINTTEXT("%c",8);
636 }
637 }
638 #endif
639
640
641
642 /**
643 * Prepares the affine RHS for KKT system.
644 * Given the special way we store the KKT matrix (sparse representation
645 * of the scalings for the second-order cone), we need this to prepare
646 * the RHS before solving the KKT system in the special format.
647 */
RHS_affine(pwork * w)648 void RHS_affine(pwork* w)
649 {
650 pfloat* RHS = w->KKT->RHS2;
651 idxint n = w->n;
652 idxint p = w->p;
653 idxint i, j, k, l;
654 idxint* Pinv = w->KKT->Pinv;
655
656 j = 0;
657 for( i=0; i < n; i++ ){ RHS[Pinv[j++]] = w->rx[i]; }
658 for( i=0; i < p; i++ ){ RHS[Pinv[j++]] = -w->ry[i]; }
659 for( i=0; i < w->C->lpc->p; i++ ){ RHS[Pinv[j++]] = w->s[i] - w->rz[i]; }
660 k = w->C->lpc->p;
661 for( l=0; l < w->C->nsoc; l++ ){
662 for( i=0; i < w->C->soc[l].p; i++ ){
663 RHS[Pinv[j++]] = w->s[k] - w->rz[k]; k++;
664 }
665 #if CONEMODE == 0
666 RHS[Pinv[j++]] = 0;
667 RHS[Pinv[j++]] = 0;
668 #endif
669 }
670 #ifdef EXPCONE
671
672 /*Exponential cones*/
673 for(l=0;l<w->C->nexc; l++){
674 for(i=0;i<3;i++)
675 {
676 RHS[Pinv[j++]] = w->s[k] - w->rz[k]; k++;
677 }
678 }
679
680 #endif
681 }
682
683
684
685 /**
686 * Prepares the RHS for computing the combined search direction.
687 */
RHS_combined(pwork * w)688 void RHS_combined(pwork* w)
689 {
690 pfloat* ds1 = w->KKT->work1;
691 pfloat* ds2 = w->KKT->work2;
692 idxint i, j, k, l;
693 pfloat sigmamu = w->info->sigma * w->info->mu;
694 pfloat one_minus_sigma = 1.0 - w->info->sigma;
695 idxint* Pinv = w->KKT->Pinv;
696 #ifdef EXPCONE
697 pfloat* s;
698 #endif
699
700 /* ds = lambda o lambda + W\s o Wz - sigma*mu*e) */
701 conicProduct(w->lambda, w->lambda, w->C, ds1);
702 conicProduct(w->dsaff_by_W, w->W_times_dzaff, w->C, ds2);
703 for( i=0; i < w->C->lpc->p; i++ ){ ds1[i] += ds2[i] - sigmamu; }
704 k = w->C->lpc->p;
705 for( i=0; i < w->C->nsoc; i++ ){
706 ds1[k] += ds2[k] - sigmamu; k++;
707 for( j=1; j < w->C->soc[i].p; j++ ){ ds1[k] += ds2[k]; k++; }
708 }
709
710 /* dz = -(1-sigma)*rz + W*(lambda \ ds) */
711 conicDivision(w->lambda, ds1, w->C, w->dsaff_by_W);
712 scale(w->dsaff_by_W, w->C, ds1);
713
714 /* copy in RHS */
715 j = 0;
716 for( i=0; i < w->n; i++ ){ w->KKT->RHS2[Pinv[j++]] *= one_minus_sigma; }
717 for( i=0; i < w->p; i++ ){ w->KKT->RHS2[Pinv[j++]] *= one_minus_sigma; }
718 for( i=0; i < w->C->lpc->p; i++) { w->KKT->RHS2[Pinv[j++]] = -one_minus_sigma*w->rz[i] + ds1[i]; }
719 k = w->C->lpc->p;
720 for( l=0; l < w->C->nsoc; l++ ){
721 for( i=0; i < w->C->soc[l].p; i++ ){
722 w->KKT->RHS2[Pinv[j++]] = -one_minus_sigma*w->rz[k] + ds1[k];
723 k++;
724 }
725 #if CONEMODE == 0
726 w->KKT->RHS2[Pinv[j++]] = 0;
727 w->KKT->RHS2[Pinv[j++]] = 0;
728 #endif
729 }
730
731 #ifdef EXPCONE
732 s = w->s;
733 k = w->C->fexv;
734
735 /*Exponential cones*/
736 /* -(1-sigma)*rz +s + mu*sigma*g(x) */
737 for(l=0;l<w->C->nexc;l++)
738 {
739 w->C->expc[l].g[0] = s[k]+sigmamu*w->C->expc[l].g[0];
740 w->KKT->RHS2[Pinv[j++]] = -one_minus_sigma*w->rz[k] + w->C->expc[l].g[0];
741 k++;
742 w->C->expc[l].g[1] = s[k]+sigmamu*w->C->expc[l].g[1];
743 w->KKT->RHS2[Pinv[j++]] = -one_minus_sigma*w->rz[k] + w->C->expc[l].g[1];
744 k++;
745 w->C->expc[l].g[2] = s[k]+sigmamu*w->C->expc[l].g[2];
746 w->KKT->RHS2[Pinv[j++]] = -one_minus_sigma*w->rz[k] + w->C->expc[l].g[2];
747 k++;
748
749 }
750
751 #endif
752 }
753
754
755
756 #ifdef EXPCONE
757 /* Line search to backtrack into feasibility for the exponential cones
758 * When affine = 1 it starts from alpha=w->info->step_aff
759 * when affine = 0 it starts from alpha=w->info->step
760 */
expConeLineSearch(pwork * w,pfloat dtau,pfloat dkappa,idxint affine)761 pfloat expConeLineSearch(pwork* w, pfloat dtau, pfloat dkappa, idxint affine)
762 {
763 /* Iterates and workspace */
764 pfloat* ws = w->KKT->work1;
765 pfloat* wz = w->KKT->work2;
766 pfloat* s = w->s;
767 pfloat* z = w->z;
768 pfloat* ds = w->dsaff;
769 pfloat* dz = w->KKT->dz2;
770 pfloat tau = w->tau;
771 pfloat kap = w->kap;
772 /* pfloat* ds = w->dsaff_by_W; */
773 /* pfloat* dz = w->W_times_dzaff; */
774
775 /* Settings */
776 pfloat gamma = w->stgs->gamma;
777
778 /* Misc */
779 idxint bk_iter, j;
780
781 pfloat mui, mu;
782
783 /* Initial alpha */
784 pfloat alpha;
785
786 /* Pointers to counters */
787 idxint* pob;
788 idxint* cb;
789 idxint* cob;
790 idxint* pb;
791 idxint* db;
792
793 /* Initialize the other statistics */
794 pfloat* centrality = &w->info->centrality;
795 pfloat barrier = 0.0; /* Variable to hold the value of f_e(s_e)+f^\star_e(z_e) */
796
797 /* Constants */
798 const pfloat cent_constant = w->D+1;
799
800 /* Placeholder */
801 *(centrality) = 1e300; /* Placeholder */
802
803 /* Start from the alpha chosen for the symmetric cones */
804 if(affine ==1)
805 {
806 alpha = w->info->step_aff;
807 }
808 else /* Combined */
809 {
810 alpha = w->info->step;
811 }
812 pob = &w->info->pob;
813 cb = &w->info->cb;
814 cob = &w->info->cob;
815 pb = &w->info->pb;
816 db = &w->info->db;
817
818 /* Initialize the counters */
819 *pob = 0;
820 *cb = 0;
821 *cob = 0;
822 *pb = 0;
823 *db = 0;
824
825 for(bk_iter = 0;bk_iter<w->stgs->max_bk_iter;bk_iter++)
826 {
827 mu = 0.0;
828
829 for(j=0;j<w->m;j++)
830 {
831 ws[j] = s[j]+alpha*ds[j];
832 wz[j] = z[j]+alpha*dz[j];
833 mu += ws[j]*wz[j];
834 }
835
836 mu = mu+(tau+alpha*dtau)*(kap+alpha*dkappa);
837 mu = mu/(w->D+1);
838
839 /* Evaluate the feasibility */
840 if(evalExpDualFeas(wz+w->C->fexv,w->C->nexc)==1)
841 {
842 if(evalExpPrimalFeas(ws+w->C->fexv,w->C->nexc)==1)
843 {
844
845 /* Compute the mu_i for each cone
846 * Because the sum mu_i = mu*nexc,
847 * Dont let any mu_i get much smaller than mu
848 */
849
850 /* Evaluate the first mui */
851 j=w->C->fexv;
852 mui = ws[j]*wz[j]+ws[j+1]*wz[j+1]+ws[j+2]*wz[j+2];
853 mui /= 3.0;
854 while(mui>MIN_DISTANCE*mu && j<w->m-2)
855 {
856 /* Calculate mui */
857 j+=3;
858 if(j<w->m)
859 {
860 mui = ws[j]*wz[j]+ws[j+1]*wz[j+1]+ws[j+2]*wz[j+2];
861 mui /= 3.0;
862 }
863
864 }
865 /* If all the mui are larger than MIN_DISTANCE*exmpu then j == m */
866 if(j==w->m)
867 {
868
869 barrier = evalBarrierValue(ws,wz, w->C->fexv, w->C->nexc);
870
871 barrier += evalSymmetricBarrierValue(ws, wz, tau+alpha*dtau, kap+alpha*dkappa, w->C, w->D);
872
873 *centrality = barrier + cent_constant*log(mu) + cent_constant;
874
875 if(*centrality<w->stgs->centrality)
876 {
877 return alpha*gamma;
878 }
879 else /* Iterate violates centrality */
880 {
881
882 #if PRINTLEVEL > 2
883 PRINTTEXT("Centrality violation %e, step %e , mu %e, barrier %e\n",*centrality,alpha,mu,barrier);
884 #endif
885 /*Count a backtrack due to centrality violation*/
886 (*cb)++;
887 }
888 }
889 else /* At least one mui is too small*/
890 {
891 #if PRINTLEVEL > 2
892 PRINTTEXT("Cone too close to complementary %e, step %e , mu %e\n",mui,alpha,mu);
893 #endif
894 /*Count a backtrack due to primal infeasibility*/
895 (*cob)++;
896 }
897
898 }
899 else
900 {
901 #if PRINTLEVEL > 2
902 PRINTTEXT("Primal infeasible at backtrack %ld\n",bk_iter);
903 /* Something is fishy find the largest entry in the iterate */
904 double max = 0;
905 int mix = 0;
906 int ik = 0;
907 for(ik = 0; ik < w->m;ik++)
908 {
909 if(ds[ik]>max)
910 {
911 max = ds[ik];
912 mix = ik;
913 }
914 if(dz[ik]>max)
915 {
916 max = dz[ik];
917 mix = ik;
918 }
919 }
920 PRINTTEXT("Max dir %lf %ld \n",max,mix);
921 #endif
922 (*pb)++;
923 }
924 }
925 else
926 {
927
928 #if PRINTLEVEL > 2
929 PRINTTEXT("Dual infeasible at backtrack %ld\n",bk_iter);
930 #endif
931 /*Count a backtrack due to dual infeasibility*/
932 (*db)++;
933 }
934 alpha = alpha*w->stgs->bk_scale;
935
936 }
937 #if PRINTLEVEL > 2
938 PRINTTEXT("Max backtrack iters reached\n");
939 #endif
940
941 return -1;
942 }
943 #endif
944 /*
945 * Line search according to Vandenberghe (cf. �8 in his manual).
946 */
lineSearch(pfloat * lambda,pfloat * ds,pfloat * dz,pfloat tau,pfloat dtau,pfloat kap,pfloat dkap,cone * C,kkt * KKT)947 pfloat lineSearch(pfloat* lambda, pfloat* ds, pfloat* dz, pfloat tau, pfloat dtau, pfloat kap, pfloat dkap, cone* C, kkt* KKT)
948 {
949 idxint i, j, cone_start, conesize;
950 pfloat rhomin, sigmamin, alpha, lknorm2, lknorm, lknorminv, rhonorm, sigmanorm, conic_step, temp;
951 pfloat lkbar_times_dsk, lkbar_times_dzk, factor;
952 pfloat* lk;
953 pfloat* dsk;
954 pfloat* dzk;
955 pfloat* lkbar = KKT->work1;
956 pfloat* rho = KKT->work2;
957 pfloat* sigma = KKT->work2;
958 pfloat minus_tau_by_dtau = -tau/dtau;
959 pfloat minus_kap_by_dkap = -kap/dkap;
960
961
962 /* LP cone */
963 if( C->lpc->p > 0 ){
964 rhomin = ds[0] / lambda[0]; sigmamin = dz[0] / lambda[0];
965 for( i=1; i < C->lpc->p; i++ ){
966 rho[0] = ds[i] / lambda[i]; if( rho[0] < rhomin ){ rhomin = rho[0]; }
967 sigma[0] = dz[i] / lambda[i]; if( sigma[0] < sigmamin ){ sigmamin = sigma[0]; }
968 }
969
970 if( -sigmamin > -rhomin ){
971 alpha = sigmamin < 0 ? 1.0 / (-sigmamin) : 1.0 / EPS;
972 } else {
973 alpha = rhomin < 0 ? 1.0 / (-rhomin) : 1.0 / EPS;
974 }
975 } else {
976 alpha = 10;
977 }
978
979 /* tau and kappa */
980 if( minus_tau_by_dtau > 0 && minus_tau_by_dtau < alpha )
981 {
982 alpha = minus_tau_by_dtau;
983 }
984 if( minus_kap_by_dkap > 0 && minus_kap_by_dkap < alpha )
985 {
986 alpha = minus_kap_by_dkap;
987 }
988
989
990 /* Second-order cone */
991 cone_start = C->lpc->p;
992 for( i=0; i < C->nsoc; i++ ){
993
994 /* indices */
995 conesize = C->soc[i].p;
996 lk = lambda + cone_start; dsk = ds + cone_start; dzk = dz + cone_start;
997
998 /* normalize */
999 lknorm2 = lk[0]*lk[0] - eddot(conesize-1, lk+1, lk+1);
1000 if (lknorm2 <= 0.0)
1001 continue;
1002
1003 lknorm = sqrt( lknorm2 );
1004 for( j=0; j < conesize; j++ ){ lkbar[j] = lk[j] / lknorm; }
1005 lknorminv = 1.0 / lknorm;
1006
1007 /* calculate products */
1008 lkbar_times_dsk = lkbar[0]*dsk[0];
1009 for( j=1; j < conesize; j++ ){ lkbar_times_dsk -= lkbar[j]*dsk[j]; }
1010 lkbar_times_dzk = lkbar[0]*dzk[0];
1011 for( j=1; j < conesize; j++ ){ lkbar_times_dzk -= lkbar[j]*dzk[j]; }
1012
1013 /* now construct rhok and sigmak, the first element is different */
1014 rho[0] = lknorminv * lkbar_times_dsk;
1015 factor = (lkbar_times_dsk+dsk[0])/(lkbar[0]+1);
1016 for( j=1; j < conesize; j++ ){ rho[j] = lknorminv*(dsk[j] - factor*lkbar[j]); }
1017 rhonorm = norm2(rho+1, conesize-1) - rho[0];
1018
1019 sigma[0] = lknorminv * lkbar_times_dzk;
1020 factor = (lkbar_times_dzk+dzk[0])/(lkbar[0]+1);
1021 for( j=1; j < conesize; j++ ){ sigma[j] = lknorminv*(dzk[j] - factor*lkbar[j]); }
1022 sigmanorm = norm2(sigma+1, conesize-1) - sigma[0];
1023
1024 /* update alpha */
1025 conic_step = 0;
1026 if( rhonorm > conic_step ){ conic_step = rhonorm; }
1027 if( sigmanorm > conic_step ){ conic_step = sigmanorm; }
1028 if( conic_step != 0 ){
1029 temp = 1.0 / conic_step;
1030 if( temp < alpha ){ alpha = temp; }
1031 }
1032
1033 cone_start += C->soc[i].p;
1034
1035 }
1036
1037 /* saturate between STEPMIN and STEPMAX */
1038 if( alpha > STEPMAX ) alpha = STEPMAX;
1039 if( alpha < STEPMIN ) alpha = STEPMIN;
1040
1041 /* return alpha */
1042 return alpha;
1043 }
1044
1045
1046
1047 /**
1048 * Scales variables by 1.0/tau, i.e. computes
1049 * x = x./tau, y = y./tau, z = z./tau, s = s./tau
1050 */
backscale(pwork * w)1051 void backscale(pwork *w)
1052 {
1053 idxint i;
1054 #if defined EQUILIBRATE && EQUILIBRATE > 0
1055 /* We performed a change of variables on x so this unsets the equilibration */
1056 for( i=0; i < w->n; i++ ){ w->x[i] /= (w->xequil[i] * w->tau); }
1057 for( i=0; i < w->p; i++ ){ w->y[i] /= (w->Aequil[i] * w->tau); }
1058 for( i=0; i < w->m; i++ ){ w->z[i] /= (w->Gequil[i] * w->tau); }
1059 for( i=0; i < w->m; i++ ){ w->s[i] *= (w->Gequil[i] / w->tau); }
1060 for( i=0; i < w->n; i++ ){ w->c[i] *= w->xequil[i]; }
1061 #else
1062 /* standard back scaling without equilibration */
1063 for( i=0; i < w->n; i++ ){ w->x[i] /= w->tau; }
1064 for( i=0; i < w->p; i++ ){ w->y[i] /= w->tau; }
1065 for( i=0; i < w->m; i++ ){ w->z[i] /= w->tau; }
1066 for( i=0; i < w->m; i++ ){ w->s[i] /= w->tau; }
1067 #endif
1068 }
1069
1070
1071
1072 /*
1073 * Main solver routine.
1074 */
ECOS_solve(pwork * w)1075 idxint ECOS_solve(pwork* w)
1076 {
1077 idxint i, initcode, KKT_FACTOR_RETURN_CODE;
1078 pfloat dtau_denom, dtauaff, dkapaff, sigma, dtau, dkap, bkap;
1079 idxint exitcode = ECOS_FATAL, interrupted = 0;
1080 pfloat pres_prev = (pfloat)ECOS_NAN;
1081
1082 #ifdef EXPCONE
1083 idxint fc = w->C->fexv; /* First cone variable e */
1084 idxint k;
1085 #endif
1086
1087 #if DEBUG
1088 char fn[20];
1089 #endif
1090
1091 #if PROFILING > 0
1092 timer tsolve;
1093 #endif
1094 #if PROFILING > 1
1095 timer tfactor, tkktsolve;
1096 #endif
1097
1098 #if defined EQUILIBRATE && EQUILIBRATE > 0
1099 for(i = 0; i <w->n; i++) { w->c[i] /= w->xequil[i]; }
1100 #endif
1101
1102 #if PROFILING > 0
1103 /* start timer */
1104 tic(&tsolve);
1105 #endif
1106
1107 /* initialize ctrl-c support */
1108 #if CTRLC > 0
1109 init_ctrlc();
1110 #endif
1111
1112 /* Initialize solver */
1113 initcode = init(w);
1114 if( initcode == ECOS_FATAL ){
1115 #if PRINTLEVEL > 0
1116 if( w->stgs->verbose ) PRINTTEXT("\nFatal error during initialization, aborting.");
1117 #endif
1118 return ECOS_FATAL;
1119 }
1120
1121
1122
1123 /* MAIN INTERIOR POINT LOOP ---------------------------------------------------------------------- */
1124 for( w->info->iter = 0; w->info->iter <= w->stgs->maxit ; w->info->iter++ ){
1125
1126 /* Compute residuals */
1127 computeResiduals(w);
1128
1129 /* Update statistics */
1130 updateStatistics(w);
1131
1132 #if PRINTLEVEL > 1
1133 /* Print info */
1134 if( w->stgs->verbose ) printProgress(w->info);
1135 #endif
1136
1137 /* SAFEGUARD: Backtrack to best previously seen iterate if
1138 *
1139 * - the update was bad such that the primal residual PRES has increased by a factor of SAFEGUARD, or
1140 * - the gap became negative
1141 *
1142 * If the safeguard is activated, the solver tests if reduced precision has been reached, and reports
1143 * accordingly. If not even reduced precision is reached, ECOS returns the flag ECOS_NUMERICS.
1144 */
1145 if( w->info->iter > 0 && (w->info->pres > SAFEGUARD*pres_prev || w->info->gap < 0) ){
1146 #if PRINTLEVEL > 1
1147 if( w->stgs->verbose ) deleteLastProgressLine( w->info );
1148 if( w->stgs->verbose ) PRINTTEXT("Unreliable search direction detected, recovering best iterate (%i) and stopping.\n", (int)w->best_info->iter);
1149 #endif
1150 restoreBestIterate( w );
1151
1152 /* Determine whether we have reached at least reduced accuracy */
1153 exitcode = checkExitConditions( w, ECOS_INACC_OFFSET );
1154
1155 /* if not, exit anyways */
1156 if( exitcode == ECOS_NOT_CONVERGED_YET ){
1157 exitcode = ECOS_NUMERICS;
1158 #if PRINTLEVEL > 0
1159 if( w->stgs->verbose ) PRINTTEXT("\nNUMERICAL PROBLEMS (reached feastol=%3.1e, reltol=%3.1e, abstol=%3.1e).", MAX(w->info->dres, w->info->pres), w->info->relgap, w->info->gap);
1160 #endif
1161 break;
1162 } else {
1163 break;
1164 }
1165 }
1166 pres_prev = w->info->pres;
1167
1168
1169 /* Check termination criteria to full precision and exit if necessary */
1170 exitcode = checkExitConditions( w, 0 );
1171 #if CTRLC > 0
1172 interrupted = check_ctrlc();
1173 #endif
1174 if( exitcode == ECOS_NOT_CONVERGED_YET ){
1175
1176 /*
1177 * Full precision has not been reached yet. Check for two more cases of exit:
1178 * (i) min step size, in which case we assume we won't make progress any more, and
1179 * (ii) maximum number of iterations reached
1180 * If these two are not fulfilled, another iteration will be made.
1181 */
1182
1183 /* Did the line search cock up? (zero step length) */
1184 if( w->info->iter > 0 && w->info->step == STEPMIN*GAMMA ){
1185 #if PRINTLEVEL > 0
1186 if( w->stgs->verbose ) deleteLastProgressLine( w->info );
1187 if( w->stgs->verbose ) PRINTTEXT("No further progress possible, recovering best iterate (%i) and stopping.", (int)w->best_info->iter );
1188 #endif
1189 restoreBestIterate( w );
1190
1191 /* Determine whether we have reached reduced precision */
1192 exitcode = checkExitConditions( w, ECOS_INACC_OFFSET );
1193 if( exitcode == ECOS_NOT_CONVERGED_YET ){
1194 exitcode = ECOS_NUMERICS;
1195 #if PRINTLEVEL > 0
1196 if( w->stgs->verbose ) PRINTTEXT("\nNUMERICAL PROBLEMS (reached feastol=%3.1e, reltol=%3.1e, abstol=%3.1e).", MAX(w->info->dres, w->info->pres), w->info->relgap, w->info->gap);
1197 #endif
1198 }
1199 break;
1200 }
1201 /* MAXIT reached? */
1202 else if( interrupted || w->info->iter == w->stgs->maxit ){
1203
1204 #if PRINTLEVEL > 0
1205 const char *what = interrupted ? "SIGINT intercepted" : "Maximum number of iterations reached";
1206 #endif
1207 /* Determine whether current iterate is better than what we had so far */
1208 if( compareStatistics( w->info, w->best_info) ){
1209 #if PRINTLEVEL > 0
1210 if( w->stgs->verbose )
1211 PRINTTEXT("%s, stopping.\n",what);
1212 #endif
1213 } else
1214 {
1215 #if PRINTLEVEL > 0
1216 if( w->stgs->verbose )
1217 PRINTTEXT("%s, recovering best iterate (%i) and stopping.\n", what, (int)w->best_info->iter);
1218 #endif
1219 restoreBestIterate( w );
1220 }
1221
1222 /* Determine whether we have reached reduced precision */
1223 exitcode = checkExitConditions( w, ECOS_INACC_OFFSET );
1224 if( exitcode == ECOS_NOT_CONVERGED_YET ){
1225 exitcode = interrupted ? ECOS_SIGINT : ECOS_MAXIT;
1226 #if PRINTLEVEL > 0
1227 if( w->stgs->verbose ) {
1228 const char* what = interrupted ? "INTERRUPTED" : "RAN OUT OF ITERATIONS";
1229 PRINTTEXT("\n%s (reached feastol=%3.1e, reltol=%3.1e, abstol=%3.1e).", what, MAX(w->info->dres, w->info->pres), w->info->relgap, w->info->gap);
1230 }
1231 #endif
1232 }
1233 break;
1234 }
1235 /* stuck on NAN? */
1236 else if( isnan(w->info->pcost) ){
1237 #if PRINTLEVEL > 0
1238 const char *what = "Reached NAN dead end";
1239 #endif
1240 /* Determine whether current iterate is better than what we had so far */
1241 if( compareStatistics( w->info, w->best_info) ){
1242 #if PRINTLEVEL > 0
1243 if( w->stgs->verbose )
1244 PRINTTEXT("%s, stopping.\n",what);
1245 #endif
1246 } else {
1247 #if PRINTLEVEL > 0
1248 if( w->stgs->verbose )
1249 PRINTTEXT("%s, recovering best iterate (%i) and stopping.\n", what, (int)w->best_info->iter);
1250 #endif
1251 restoreBestIterate( w );
1252 }
1253
1254 /* Determine whether we have reached reduced precision */
1255 exitcode = checkExitConditions( w, ECOS_INACC_OFFSET );
1256 if( exitcode == ECOS_NOT_CONVERGED_YET ){
1257 exitcode = ECOS_NUMERICS;
1258 #if PRINTLEVEL > 0
1259 if( w->stgs->verbose ) {
1260 const char* what = interrupted ? "INTERRUPTED" : "RAN OUT OF ITERATIONS";
1261 PRINTTEXT("\n%s (reached feastol=%3.1e, reltol=%3.1e, abstol=%3.1e).", what, MAX(w->info->dres, w->info->pres), w->info->relgap, w->info->gap);
1262 }
1263 #endif
1264 }
1265 break;
1266 }
1267 } else {
1268
1269 /* Full precision has been reached, stop solver */
1270 break;
1271 }
1272
1273
1274
1275 /* SAFEGUARD:
1276 * Check whether current iterate is worth keeping as the best solution so far,
1277 * before doing another iteration
1278 */
1279 if (w->info->iter == 0) {
1280 /* we're at the first iterate, so there's nothing to compare yet */
1281 saveIterateAsBest( w );
1282 } else if( compareStatistics( w->info, w->best_info) ){
1283 /* PRINTTEXT("Better solution found, saving as best so far \n"); */
1284 saveIterateAsBest( w );
1285 }
1286
1287
1288 /* Compute scalings */
1289 #ifdef EXPCONE
1290 if( updateScalings(w->C, w->s, w->z, w->lambda, w->info->mu) == OUTSIDE_CONE ){
1291 #else
1292 if( updateScalings(w->C, w->s, w->z, w->lambda) == OUTSIDE_CONE ){
1293 #endif
1294 /* SAFEGUARD: we have to recover here */
1295 #if PRINTLEVEL > 0
1296 if( w->stgs->verbose ) deleteLastProgressLine( w->info );
1297 if( w->stgs->verbose ) PRINTTEXT("Slacks/multipliers leaving the cone, recovering best iterate (%i) and stopping.\n", (int)w->best_info->iter);
1298 #endif
1299 restoreBestIterate( w );
1300
1301 /* Determine whether we have reached at least reduced accuracy */
1302 exitcode = checkExitConditions( w, ECOS_INACC_OFFSET );
1303 if( exitcode == ECOS_NOT_CONVERGED_YET ){
1304 #if PRINTLEVEL > 0
1305 if( w->stgs->verbose ) PRINTTEXT("\nNUMERICAL PROBLEMS (reached feastol=%3.1e, reltol=%3.1e, abstol=%3.1e).", MAX(w->info->dres, w->info->pres), w->info->relgap, w->info->gap);
1306 #endif
1307 return ECOS_OUTCONE;
1308
1309 } else {
1310 break;
1311 }
1312 }
1313
1314 /* Update KKT matrix with scalings */
1315 kkt_update(w->KKT->PKPt, w->KKT->PK, w->C);
1316
1317 #if DEBUG > 0
1318 /* DEBUG: Store matrix to be factored */
1319 sprintf(fn, "PKPt_updated_%02i.txt", (int)w->info->iter);
1320 dumpSparseMatrix(w->KKT->PKPt, fn);
1321 #endif
1322 /* factor KKT matrix */
1323 #if PROFILING > 1
1324 tic(&tfactor);
1325 KKT_FACTOR_RETURN_CODE = kkt_factor(w->KKT, w->stgs->eps, w->stgs->delta, &w->info->tfactor_t1, &w->info->tfactor_t2);
1326 w->info->tfactor += toc(&tfactor);
1327 #else
1328 KKT_FACTOR_RETURN_CODE = kkt_factor(w->KKT, w->stgs->eps, w->stgs->delta);
1329 #endif
1330
1331 #if DEBUG > 0
1332 /* DEBUG: store factor */
1333 sprintf(fn, "PKPt_factor_%02i.txt", (int)w->info->iter);
1334 dumpSparseMatrix(w->KKT->L, fn);
1335 #endif
1336
1337 /* check if factorization was successful, exit otherwise */
1338 if( KKT_FACTOR_RETURN_CODE != KKT_OK ){
1339 #if PRINTLEVEL > 0
1340 if( w->stgs->verbose ) PRINTTEXT("\nProblem in factoring KKT system, aborting.");
1341 #endif
1342 return ECOS_FATAL;
1343 }
1344
1345 /* Solve for RHS1, which is used later also in combined direction */
1346 #if PROFILING > 1
1347 tic(&tkktsolve);
1348 #endif
1349 w->info->nitref1 = kkt_solve(w->KKT, w->A, w->G, w->KKT->RHS1, w->KKT->dx1, w->KKT->dy1, w->KKT->dz1, w->n, w->p, w->m, w->C, 0, w->stgs->nitref);
1350 #if PROFILING > 1
1351 w->info->tkktsolve += toc(&tkktsolve);
1352 #endif
1353
1354 #if DEBUG > 0 && PRINTLEVEL > 2
1355 /* Print result of linear system solve */
1356 printDenseMatrix(w->KKT->dx1, 1, 5, "dx1(1:5)");
1357 printDenseMatrix(w->KKT->dy1, 1, 5, "dy1(1:5)");
1358 printDenseMatrix(w->KKT->dz1, 1, 5, "dz1(1:5)");
1359 #endif
1360
1361 /* AFFINE SEARCH DIRECTION (predictor, need dsaff and dzaff only) */
1362 RHS_affine(w);
1363 #if PROFILING > 1
1364 tic(&tkktsolve);
1365 #endif
1366 w->info->nitref2 = kkt_solve(w->KKT, w->A, w->G, w->KKT->RHS2, w->KKT->dx2, w->KKT->dy2, w->KKT->dz2, w->n, w->p, w->m, w->C, 0, w->stgs->nitref);
1367 #if PROFILING > 1
1368 w->info->tkktsolve += toc(&tkktsolve);
1369 #endif
1370
1371 /* dtau_denom = kap/tau - (c'*x1 + by1 + h'*z1); */
1372 dtau_denom = w->kap/w->tau - eddot(w->n, w->c, w->KKT->dx1) - eddot(w->p, w->b, w->KKT->dy1) - eddot(w->m, w->h, w->KKT->dz1);
1373
1374 /* dtauaff = (dt + c'*x2 + by2 + h'*z2) / dtau_denom; */
1375 dtauaff = (w->rt - w->kap + eddot(w->n, w->c, w->KKT->dx2) + eddot(w->p, w->b, w->KKT->dy2) + eddot(w->m, w->h, w->KKT->dz2)) / dtau_denom;
1376
1377 /* dzaff = dz2 + dtau_aff*dz1 */
1378 /* let dz2 = dzaff we use this in the linesearch for unsymmetric cones*/
1379 /* and w_times_dzaff = Wdz_aff*/
1380 /* and dz2 = dz2+dtau_aff*dz1 will store the unscaled dz*/
1381 for( i=0; i<w->m; i++ ){ w->KKT->dz2[i] = w->KKT->dz2[i] + dtauaff*w->KKT->dz1[i]; }
1382 scale(w->KKT->dz2, w->C, w->W_times_dzaff);
1383
1384 /* W\dsaff = -W*dzaff -lambda; */
1385 for( i=0; i<w->m; i++ ){ w->dsaff_by_W[i] = -w->W_times_dzaff[i] - w->lambda[i]; }
1386
1387 /* dkapaff = -(bkap + kap*dtauaff)/tau; bkap = kap*tau*/
1388 dkapaff = -w->kap - w->kap/w->tau*dtauaff;
1389
1390 /* Line search on W\dsaff and W*dzaff */
1391 w->info->step_aff = lineSearch(w->lambda, w->dsaff_by_W, w->W_times_dzaff, w->tau, dtauaff, w->kap, dkapaff, w->C, w->KKT);
1392
1393 #ifdef EXPCONE
1394
1395 /* The present value of step_aff keeps the primal and dual symmetric slacks close to the boundary
1396 * Here w->W_times_dzaff contains dzaff for the exponential cones
1397 * However W\dsaff is not initialized for the exponential cones.
1398 * first we must calculate dsaff
1399 * From \muH(x)dxaff + dsaff = -s we set dsaff to be
1400 * dsaff = -s - \mu H(z)dzaff
1401 */
1402
1403 /* Calculate muH(s)dzaff and store in dsaff_by_W */
1404
1405 /* For the linesearch we need to compute x's+tk at every step size so
1406 * we need to calculate the search direction for the symmetric cones too
1407 */
1408
1409 /* ds = W*ds_by_W */
1410 scale(w->dsaff_by_W, w->C, w->dsaff);
1411
1412 /* Zero out */
1413 for(i=fc; i<w->m; i++)
1414 {
1415 w->dsaff[i] = 0.0;
1416 }
1417
1418 scaleToAddExpcone(w->dsaff,w->KKT->dz2, w->C->expc, w->C->nexc, fc);
1419
1420 /* Add -s ds = -s - muH(z)dz */
1421 for(i=fc; i<w->m; i++)
1422 {
1423 w->dsaff[i] = -w->dsaff[i]-w->s[i];
1424 }
1425 /* Now dsaff_by_W has dsaff and W_times_dzaff dzaff */
1426
1427
1428 #if PRINTLEVEL > 2
1429 PRINTTEXT("Backtrack affine\n");
1430 #endif
1431 /* Set the affine backtracking counter to zero */
1432 w->info->affBack = 0;
1433 if(w->C->nexc>0)
1434 {
1435 w->info->step_aff = expConeLineSearch(w,dtauaff,dkapaff,1);
1436 /* Sum all the backtracking steps */
1437 /* Sum the backtracking counts */
1438 w->info->affBack+= w->info->pob;
1439 w->info->affBack+=w->info->cb;
1440 w->info->affBack+= w->info->cob;
1441 w->info->affBack+=w->info->pb;
1442 w->info->affBack+=w->info->db;
1443
1444 #if PRINTLEVEL > 2
1445 PRINTTEXT("Affine backtracking %i %i %i %i %i \n",\
1446 (int)w->info->pob,\
1447 (int)w->info->cb,\
1448 (int)w->info->cob,\
1449 (int)w->info->pb,\
1450 (int)w->info->db);
1451 #endif
1452 /* Detect linesearch failure, recover and end. */
1453 if(w->info->step_aff == -1)
1454 {
1455 #if PRINTLEVEL > 2
1456 PRINTTEXT("Affine backtrack failed\n");
1457 #endif
1458 w->info->step_aff = 0.0;
1459 }
1460
1461 }
1462 #endif /* END EXPCONE */
1463
1464 /* Centering parameter */
1465 sigma = 1.0 - w->info->step_aff;
1466 sigma = sigma*sigma*sigma;
1467 if( sigma > SIGMAMAX ) sigma = SIGMAMAX;
1468 if( sigma < SIGMAMIN ) sigma = SIGMAMIN;
1469 w->info->sigma = sigma;
1470
1471
1472 /* COMBINED SEARCH DIRECTION */
1473 RHS_combined(w);
1474 #if PROFILING > 1
1475 tic(&tkktsolve);
1476 #endif
1477 w->info->nitref3 = kkt_solve(w->KKT, w->A, w->G, w->KKT->RHS2, w->KKT->dx2, w->KKT->dy2, w->KKT->dz2, w->n, w->p, w->m, w->C, 0, w->stgs->nitref);
1478 #if PROFILING > 1
1479 w->info->tkktsolve += toc(&tkktsolve);
1480 #endif
1481
1482 /* bkap = kap*tau + dkapaff*dtauaff - sigma*info.mu; */
1483 bkap = w->kap*w->tau + dkapaff*dtauaff - sigma*w->info->mu;
1484
1485 /* dtau = ((1-sigma)*rt - bkap/tau + c'*x2 + by2 + h'*z2) / dtau_denom; */
1486 dtau = ((1-sigma)*w->rt - bkap/w->tau + eddot(w->n, w->c, w->KKT->dx2) + eddot(w->p, w->b, w->KKT->dy2) + eddot(w->m, w->h, w->KKT->dz2)) / dtau_denom;
1487
1488 /* dx = x2 + dtau*x1; dy = y2 + dtau*y1; dz = z2 + dtau*z1; */
1489 for( i=0; i < w->n; i++ ){ w->KKT->dx2[i] += dtau*w->KKT->dx1[i]; }
1490 for( i=0; i < w->p; i++ ){ w->KKT->dy2[i] += dtau*w->KKT->dy1[i]; }
1491 for( i=0; i < w->m; i++ ){ w->KKT->dz2[i] += dtau*w->KKT->dz1[i]; }
1492
1493 /* ds_by_W = -(lambda \ bs + conelp_timesW(scaling,dz,dims)); */
1494 /* note that at this point w->dsaff_by_W holds already (lambda \ ds) */
1495 scale(w->KKT->dz2, w->C, w->W_times_dzaff);
1496 for( i=0; i < w->m; i++ ){ w->dsaff_by_W[i] = -(w->dsaff_by_W[i] + w->W_times_dzaff[i]); }
1497
1498 /* dkap = -(bkap + kap*dtau)/tau; */
1499 dkap = -(bkap + w->kap*dtau)/w->tau;
1500
1501 /* Line search on combined direction */
1502 w->info->step = lineSearch(w->lambda, w->dsaff_by_W, w->W_times_dzaff, w->tau, dtau, w->kap, dkap, w->C, w->KKT) * w->stgs->gamma;
1503
1504 /* Bring ds to the final unscaled form */
1505 /* ds = W*ds_by_W */
1506 scale(w->dsaff_by_W, w->C, w->dsaff);
1507
1508 #ifdef EXPCONE
1509
1510 /* Zero out the section corresponding to the exponential vars */
1511 for(i=fc; i<w->m; i++)
1512 {
1513 w->dsaff[i] = 0.0;
1514 }
1515 /* Calculate muH(z)dzaff and store in dsaff */
1516 scaleToAddExpcone(w->dsaff,w->KKT->dz2,w->C->expc,w->C->nexc,fc);
1517
1518 /* We have modified g to contain s+sigma*mu*g(z) */
1519 k = fc;
1520 for(i=0;i<w->C->nexc;i++)
1521 {
1522 w->dsaff[k] = -w->dsaff[k] - w->C->expc[i].g[0];
1523 k++;
1524 w->dsaff[k] = -w->dsaff[k] - w->C->expc[i].g[1];
1525 k++;
1526 w->dsaff[k] = -w->dsaff[k] - w->C->expc[i].g[2];
1527 k++;
1528 }
1529
1530 /* Set the backtracking counter to zero */
1531 w->info->cmbBack = 0;
1532 if(w->C->nexc>0)
1533 {
1534 w->info->step = expConeLineSearch(w,dtau,dkap,0);
1535 /* Sum the backtracking counts */
1536 w->info->cmbBack+= w->info->cb;
1537 w->info->cmbBack+= w->info->cob;
1538 w->info->cmbBack+= w->info->pb;
1539 w->info->cmbBack+= w->info->db;
1540
1541 #if PRINTLEVEL > 2
1542 PRINTTEXT("Combined backtracking %i %i %i %i \n",\
1543 (int)w->info->cb,\
1544 (int)w->info->cob,\
1545 (int)w->info->pb,\
1546 (int)w->info->db);
1547 #endif
1548
1549 /* Detect linesearch failure, recover and end. */
1550 if(w->info->step == -1)
1551 {
1552
1553 #if PRINTLEVEL > 1
1554 if( w->stgs->verbose ){
1555 PRINTTEXT("Combined backtracking failed %i %i %i %i sigma %g\n",\
1556 (int)w->info->cb,\
1557 (int)w->info->cob,\
1558 (int)w->info->pb,\
1559 (int)w->info->db,\
1560 sigma);
1561
1562 PRINTTEXT("Combined line search failed, recovering best iterate (%i) and stopping.\n",\
1563 (int)w->best_info->iter);
1564 }
1565 #endif
1566 restoreBestIterate( w );
1567
1568 /* Determine whether we have reached at least reduced accuracy */
1569 exitcode = checkExitConditions( w, ECOS_INACC_OFFSET );
1570
1571 /* if not, exit anyways */
1572 if( exitcode == ECOS_NOT_CONVERGED_YET ){
1573 exitcode = ECOS_NUMERICS;
1574 #if PRINTLEVEL > 0
1575 if( w->stgs->verbose ) PRINTTEXT("\nNUMERICAL PROBLEMS (reached feastol=%3.1e, reltol=%3.1e, abstol=%3.1e).", MAX(w->info->dres, w->info->pres), w->info->relgap, w->info->gap);
1576 #endif
1577 break;
1578 } else {
1579 break;
1580 }
1581 }
1582
1583 }
1584
1585 #endif /* endif of EXPCONE */
1586
1587
1588 /* Update variables */
1589 for( i=0; i < w->n; i++ ){ w->x[i] += w->info->step * w->KKT->dx2[i]; }
1590 for( i=0; i < w->p; i++ ){ w->y[i] += w->info->step * w->KKT->dy2[i]; }
1591 for( i=0; i < w->m; i++ ){ w->z[i] += w->info->step * w->KKT->dz2[i]; }
1592 for( i=0; i < w->m; i++ ){ w->s[i] += w->info->step * w->dsaff[i]; }
1593 w->kap += w->info->step * dkap;
1594 w->tau += w->info->step * dtau;
1595 }
1596
1597 /* scale variables back */
1598 backscale(w);
1599
1600 /* stop timer */
1601 #if PROFILING > 0
1602 w->info->tsolve = toc(&tsolve);
1603 #endif
1604
1605 #if PRINTLEVEL > 0
1606 #if PROFILING > 0
1607 if( w->stgs->verbose ) PRINTTEXT("\nRuntime: %f seconds.", w->info->tsetup + w->info->tsolve);
1608 #endif
1609 if( w->stgs->verbose ) PRINTTEXT("\n\n");
1610 #endif
1611
1612 #if CTRLC > 0
1613 remove_ctrlc();
1614 #endif
1615 return exitcode;
1616 }
1617
1618
1619 /*
1620 * Updates one element of the RHS vector h of inequalities
1621 * After the call, w->h[idx] = value (but equilibrated)
1622 */
1623 void ecos_updateDataEntry_h(pwork* w, idxint idx, pfloat value)
1624 {
1625 #if EQUILIBRATE > 0
1626 w->h[idx] = value / w->Gequil[idx];
1627 #else
1628 w->h[idx] = value;
1629 #endif
1630 }
1631
1632
1633 /*
1634 * Updates one element of the OBJ vector c of inequalities
1635 * After the call, w->c[idx] = value (but equilibrated)
1636 */
1637 void ecos_updateDataEntry_c(pwork* w, idxint idx, pfloat value)
1638 {
1639 w->c[idx] = value;
1640 }
1641
1642
1643 /*
1644 * Updates numerical data for G, A, c, h, and b,
1645 * and re-equilibrates.
1646 * Then updates the corresponding KKT entries.
1647 */
1648 void ECOS_updateData(pwork *w, pfloat *Gpr, pfloat *Apr,
1649 pfloat* c, pfloat* h, pfloat* b)
1650 {
1651
1652 idxint k;
1653
1654 #if defined EQUILIBRATE && EQUILIBRATE > 0
1655 /* Only unequilibrate the old data if the new data is in a different location in memory. */
1656 /* This is required for backward compatibility since existing code might need this step. */
1657 if (
1658 ((Gpr + w->G->nnz < w->G->pr) || (w->G->pr + w->G->nnz < Gpr)) &&
1659 ((Apr + w->A->nnz < w->A->pr) || (w->A->pr + w->A->nnz < Apr)) &&
1660 ((c + w->n < w->c) || (w->c + w->n < c)) &&
1661 ((h + w->m < w->h) || (w->h + w->m < h)) &&
1662 ((b + w->p < w->b) || (w->b + w->p < b))
1663 ){
1664 unset_equilibration(w);
1665 }
1666 #endif
1667
1668 /* update pointers */
1669 if(w->G){
1670 w->G->pr = Gpr;
1671 w->h = h;
1672 }
1673 if(w->A){
1674 w->A->pr = Apr;
1675 w->b = b;
1676 }
1677 w->c = c;
1678
1679 #if defined EQUILIBRATE && EQUILIBRATE > 0
1680 /* Here, we need to equilibrate unconditionally since all new data is unequilibrated. */
1681 set_equilibration(w);
1682 #endif
1683
1684 /* update KKT matrix */
1685 if(w->A){
1686 for (k=0; k<w->A->nnz; k++){
1687 w->KKT->PKPt->pr[w->KKT->PK[w->AtoK[k]]] = Apr[k];
1688 }
1689 }
1690 if(w->G){
1691 for (k=0; k<w->G->nnz; k++){
1692 w->KKT->PKPt->pr[w->KKT->PK[w->GtoK[k]]] = Gpr[k];
1693 }
1694 }
1695 }