lemon-project-template-glpk
view deps/glpk/src/glpipm.c @ 9:33de93886c88
Import GLPK 4.47
| author | Alpar Juttner <alpar@cs.elte.hu> |
|---|---|
| date | Sun, 06 Nov 2011 20:59:10 +0100 |
| parents | |
| children |
line source
1 /* glpipm.c */
3 /***********************************************************************
4 * This code is part of GLPK (GNU Linear Programming Kit).
5 *
6 * Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
7 * 2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
8 * Moscow Aviation Institute, Moscow, Russia. All rights reserved.
9 * E-mail: <mao@gnu.org>.
10 *
11 * GLPK is free software: you can redistribute it and/or modify it
12 * under the terms of the GNU General Public License as published by
13 * the Free Software Foundation, either version 3 of the License, or
14 * (at your option) any later version.
15 *
16 * GLPK is distributed in the hope that it will be useful, but WITHOUT
17 * ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
18 * or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
19 * License for more details.
20 *
21 * You should have received a copy of the GNU General Public License
22 * along with GLPK. If not, see <http://www.gnu.org/licenses/>.
23 ***********************************************************************/
25 #include "glpipm.h"
26 #include "glpmat.h"
28 #define ITER_MAX 100
29 /* maximal number of iterations */
31 struct csa
32 { /* common storage area */
33 /*--------------------------------------------------------------*/
34 /* LP data */
35 int m;
36 /* number of rows (equality constraints) */
37 int n;
38 /* number of columns (structural variables) */
39 int *A_ptr; /* int A_ptr[1+m+1]; */
40 int *A_ind; /* int A_ind[A_ptr[m+1]]; */
41 double *A_val; /* double A_val[A_ptr[m+1]]; */
42 /* mxn-matrix A in storage-by-rows format */
43 double *b; /* double b[1+m]; */
44 /* m-vector b of right-hand sides */
45 double *c; /* double c[1+n]; */
46 /* n-vector c of objective coefficients; c[0] is constant term of
47 the objective function */
48 /*--------------------------------------------------------------*/
49 /* LP solution */
50 double *x; /* double x[1+n]; */
51 double *y; /* double y[1+m]; */
52 double *z; /* double z[1+n]; */
53 /* current point in primal-dual space; the best point on exit */
54 /*--------------------------------------------------------------*/
55 /* control parameters */
56 const glp_iptcp *parm;
57 /*--------------------------------------------------------------*/
58 /* working arrays and variables */
59 double *D; /* double D[1+n]; */
60 /* diagonal nxn-matrix D = X*inv(Z), where X = diag(x[j]) and
61 Z = diag(z[j]) */
62 int *P; /* int P[1+m+m]; */
63 /* permutation mxm-matrix P used to minimize fill-in in Cholesky
64 factorization */
65 int *S_ptr; /* int S_ptr[1+m+1]; */
66 int *S_ind; /* int S_ind[S_ptr[m+1]]; */
67 double *S_val; /* double S_val[S_ptr[m+1]]; */
68 double *S_diag; /* double S_diag[1+m]; */
69 /* symmetric mxm-matrix S = P*A*D*A'*P' whose upper triangular
70 part without diagonal elements is stored in S_ptr, S_ind, and
71 S_val in storage-by-rows format, diagonal elements are stored
72 in S_diag */
73 int *U_ptr; /* int U_ptr[1+m+1]; */
74 int *U_ind; /* int U_ind[U_ptr[m+1]]; */
75 double *U_val; /* double U_val[U_ptr[m+1]]; */
76 double *U_diag; /* double U_diag[1+m]; */
77 /* upper triangular mxm-matrix U defining Cholesky factorization
78 S = U'*U; its non-diagonal elements are stored in U_ptr, U_ind,
79 U_val in storage-by-rows format, diagonal elements are stored
80 in U_diag */
81 int iter;
82 /* iteration number (0, 1, 2, ...); iter = 0 corresponds to the
83 initial point */
84 double obj;
85 /* current value of the objective function */
86 double rpi;
87 /* relative primal infeasibility rpi = ||A*x-b||/(1+||b||) */
88 double rdi;
89 /* relative dual infeasibility rdi = ||A'*y+z-c||/(1+||c||) */
90 double gap;
91 /* primal-dual gap = |c'*x-b'*y|/(1+|c'*x|) which is a relative
92 difference between primal and dual objective functions */
93 double phi;
94 /* merit function phi = ||A*x-b||/max(1,||b||) +
95 + ||A'*y+z-c||/max(1,||c||) +
96 + |c'*x-b'*y|/max(1,||b||,||c||) */
97 double mu;
98 /* duality measure mu = x'*z/n (used as barrier parameter) */
99 double rmu;
100 /* rmu = max(||A*x-b||,||A'*y+z-c||)/mu */
101 double rmu0;
102 /* the initial value of rmu on iteration 0 */
103 double *phi_min; /* double phi_min[1+ITER_MAX]; */
104 /* phi_min[k] = min(phi[k]), where phi[k] is the value of phi on
105 k-th iteration, 0 <= k <= iter */
106 int best_iter;
107 /* iteration number, on which the value of phi reached its best
108 (minimal) value */
109 double *best_x; /* double best_x[1+n]; */
110 double *best_y; /* double best_y[1+m]; */
111 double *best_z; /* double best_z[1+n]; */
112 /* best point (in the sense of the merit function phi) which has
113 been reached on iteration iter_best */
114 double best_obj;
115 /* objective value at the best point */
116 double *dx_aff; /* double dx_aff[1+n]; */
117 double *dy_aff; /* double dy_aff[1+m]; */
118 double *dz_aff; /* double dz_aff[1+n]; */
119 /* affine scaling direction */
120 double alfa_aff_p, alfa_aff_d;
121 /* maximal primal and dual stepsizes in affine scaling direction,
122 on which x and z are still non-negative */
123 double mu_aff;
124 /* duality measure mu_aff = x_aff'*z_aff/n in the boundary point
125 x_aff' = x+alfa_aff_p*dx_aff, z_aff' = z+alfa_aff_d*dz_aff */
126 double sigma;
127 /* Mehrotra's heuristic parameter (0 <= sigma <= 1) */
128 double *dx_cc; /* double dx_cc[1+n]; */
129 double *dy_cc; /* double dy_cc[1+m]; */
130 double *dz_cc; /* double dz_cc[1+n]; */
131 /* centering corrector direction */
132 double *dx; /* double dx[1+n]; */
133 double *dy; /* double dy[1+m]; */
134 double *dz; /* double dz[1+n]; */
135 /* final combined direction dx = dx_aff+dx_cc, dy = dy_aff+dy_cc,
136 dz = dz_aff+dz_cc */
137 double alfa_max_p;
138 double alfa_max_d;
139 /* maximal primal and dual stepsizes in combined direction, on
140 which x and z are still non-negative */
141 };
143 /***********************************************************************
144 * initialize - allocate and initialize common storage area
145 *
146 * This routine allocates and initializes the common storage area (CSA)
147 * used by interior-point method routines. */
149 static void initialize(struct csa *csa)
150 { int m = csa->m;
151 int n = csa->n;
152 int i;
153 if (csa->parm->msg_lev >= GLP_MSG_ALL)
154 xprintf("Matrix A has %d non-zeros\n", csa->A_ptr[m+1]-1);
155 csa->D = xcalloc(1+n, sizeof(double));
156 /* P := I */
157 csa->P = xcalloc(1+m+m, sizeof(int));
158 for (i = 1; i <= m; i++) csa->P[i] = csa->P[m+i] = i;
159 /* S := A*A', symbolically */
160 csa->S_ptr = xcalloc(1+m+1, sizeof(int));
161 csa->S_ind = adat_symbolic(m, n, csa->P, csa->A_ptr, csa->A_ind,
162 csa->S_ptr);
163 if (csa->parm->msg_lev >= GLP_MSG_ALL)
164 xprintf("Matrix S = A*A' has %d non-zeros (upper triangle)\n",
165 csa->S_ptr[m+1]-1 + m);
166 /* determine P using specified ordering algorithm */
167 if (csa->parm->ord_alg == GLP_ORD_NONE)
168 { if (csa->parm->msg_lev >= GLP_MSG_ALL)
169 xprintf("Original ordering is being used\n");
170 for (i = 1; i <= m; i++)
171 csa->P[i] = csa->P[m+i] = i;
172 }
173 else if (csa->parm->ord_alg == GLP_ORD_QMD)
174 { if (csa->parm->msg_lev >= GLP_MSG_ALL)
175 xprintf("Minimum degree ordering (QMD)...\n");
176 min_degree(m, csa->S_ptr, csa->S_ind, csa->P);
177 }
178 else if (csa->parm->ord_alg == GLP_ORD_AMD)
179 { if (csa->parm->msg_lev >= GLP_MSG_ALL)
180 xprintf("Approximate minimum degree ordering (AMD)...\n");
181 amd_order1(m, csa->S_ptr, csa->S_ind, csa->P);
182 }
183 else if (csa->parm->ord_alg == GLP_ORD_SYMAMD)
184 { if (csa->parm->msg_lev >= GLP_MSG_ALL)
185 xprintf("Approximate minimum degree ordering (SYMAMD)...\n")
186 ;
187 symamd_ord(m, csa->S_ptr, csa->S_ind, csa->P);
188 }
189 else
190 xassert(csa != csa);
191 /* S := P*A*A'*P', symbolically */
192 xfree(csa->S_ind);
193 csa->S_ind = adat_symbolic(m, n, csa->P, csa->A_ptr, csa->A_ind,
194 csa->S_ptr);
195 csa->S_val = xcalloc(csa->S_ptr[m+1], sizeof(double));
196 csa->S_diag = xcalloc(1+m, sizeof(double));
197 /* compute Cholesky factorization S = U'*U, symbolically */
198 if (csa->parm->msg_lev >= GLP_MSG_ALL)
199 xprintf("Computing Cholesky factorization S = L*L'...\n");
200 csa->U_ptr = xcalloc(1+m+1, sizeof(int));
201 csa->U_ind = chol_symbolic(m, csa->S_ptr, csa->S_ind, csa->U_ptr);
202 if (csa->parm->msg_lev >= GLP_MSG_ALL)
203 xprintf("Matrix L has %d non-zeros\n", csa->U_ptr[m+1]-1 + m);
204 csa->U_val = xcalloc(csa->U_ptr[m+1], sizeof(double));
205 csa->U_diag = xcalloc(1+m, sizeof(double));
206 csa->iter = 0;
207 csa->obj = 0.0;
208 csa->rpi = 0.0;
209 csa->rdi = 0.0;
210 csa->gap = 0.0;
211 csa->phi = 0.0;
212 csa->mu = 0.0;
213 csa->rmu = 0.0;
214 csa->rmu0 = 0.0;
215 csa->phi_min = xcalloc(1+ITER_MAX, sizeof(double));
216 csa->best_iter = 0;
217 csa->best_x = xcalloc(1+n, sizeof(double));
218 csa->best_y = xcalloc(1+m, sizeof(double));
219 csa->best_z = xcalloc(1+n, sizeof(double));
220 csa->best_obj = 0.0;
221 csa->dx_aff = xcalloc(1+n, sizeof(double));
222 csa->dy_aff = xcalloc(1+m, sizeof(double));
223 csa->dz_aff = xcalloc(1+n, sizeof(double));
224 csa->alfa_aff_p = 0.0;
225 csa->alfa_aff_d = 0.0;
226 csa->mu_aff = 0.0;
227 csa->sigma = 0.0;
228 csa->dx_cc = xcalloc(1+n, sizeof(double));
229 csa->dy_cc = xcalloc(1+m, sizeof(double));
230 csa->dz_cc = xcalloc(1+n, sizeof(double));
231 csa->dx = csa->dx_aff;
232 csa->dy = csa->dy_aff;
233 csa->dz = csa->dz_aff;
234 csa->alfa_max_p = 0.0;
235 csa->alfa_max_d = 0.0;
236 return;
237 }
239 /***********************************************************************
240 * A_by_vec - compute y = A*x
241 *
242 * This routine computes matrix-vector product y = A*x, where A is the
243 * constraint matrix. */
245 static void A_by_vec(struct csa *csa, double x[], double y[])
246 { /* compute y = A*x */
247 int m = csa->m;
248 int *A_ptr = csa->A_ptr;
249 int *A_ind = csa->A_ind;
250 double *A_val = csa->A_val;
251 int i, t, beg, end;
252 double temp;
253 for (i = 1; i <= m; i++)
254 { temp = 0.0;
255 beg = A_ptr[i], end = A_ptr[i+1];
256 for (t = beg; t < end; t++) temp += A_val[t] * x[A_ind[t]];
257 y[i] = temp;
258 }
259 return;
260 }
262 /***********************************************************************
263 * AT_by_vec - compute y = A'*x
264 *
265 * This routine computes matrix-vector product y = A'*x, where A' is a
266 * matrix transposed to the constraint matrix A. */
268 static void AT_by_vec(struct csa *csa, double x[], double y[])
269 { /* compute y = A'*x, where A' is transposed to A */
270 int m = csa->m;
271 int n = csa->n;
272 int *A_ptr = csa->A_ptr;
273 int *A_ind = csa->A_ind;
274 double *A_val = csa->A_val;
275 int i, j, t, beg, end;
276 double temp;
277 for (j = 1; j <= n; j++) y[j] = 0.0;
278 for (i = 1; i <= m; i++)
279 { temp = x[i];
280 if (temp == 0.0) continue;
281 beg = A_ptr[i], end = A_ptr[i+1];
282 for (t = beg; t < end; t++) y[A_ind[t]] += A_val[t] * temp;
283 }
284 return;
285 }
287 /***********************************************************************
288 * decomp_NE - numeric factorization of matrix S = P*A*D*A'*P'
289 *
290 * This routine implements numeric phase of Cholesky factorization of
291 * the matrix S = P*A*D*A'*P', which is a permuted matrix of the normal
292 * equation system. Matrix D is assumed to be already computed. */
294 static void decomp_NE(struct csa *csa)
295 { adat_numeric(csa->m, csa->n, csa->P, csa->A_ptr, csa->A_ind,
296 csa->A_val, csa->D, csa->S_ptr, csa->S_ind, csa->S_val,
297 csa->S_diag);
298 chol_numeric(csa->m, csa->S_ptr, csa->S_ind, csa->S_val,
299 csa->S_diag, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag);
300 return;
301 }
303 /***********************************************************************
304 * solve_NE - solve normal equation system
305 *
306 * This routine solves the normal equation system:
307 *
308 * A*D*A'*y = h.
309 *
310 * It is assumed that the matrix A*D*A' has been previously factorized
311 * by the routine decomp_NE.
312 *
313 * On entry the array y contains the vector of right-hand sides h. On
314 * exit this array contains the computed vector of unknowns y.
315 *
316 * Once the vector y has been computed the routine checks for numeric
317 * stability. If the residual vector:
318 *
319 * r = A*D*A'*y - h
320 *
321 * is relatively small, the routine returns zero, otherwise non-zero is
322 * returned. */
324 static int solve_NE(struct csa *csa, double y[])
325 { int m = csa->m;
326 int n = csa->n;
327 int *P = csa->P;
328 int i, j, ret = 0;
329 double *h, *r, *w;
330 /* save vector of right-hand sides h */
331 h = xcalloc(1+m, sizeof(double));
332 for (i = 1; i <= m; i++) h[i] = y[i];
333 /* solve normal equation system (A*D*A')*y = h */
334 /* since S = P*A*D*A'*P' = U'*U, then A*D*A' = P'*U'*U*P, so we
335 have inv(A*D*A') = P'*inv(U)*inv(U')*P */
336 /* w := P*h */
337 w = xcalloc(1+m, sizeof(double));
338 for (i = 1; i <= m; i++) w[i] = y[P[i]];
339 /* w := inv(U')*w */
340 ut_solve(m, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag, w);
341 /* w := inv(U)*w */
342 u_solve(m, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag, w);
343 /* y := P'*w */
344 for (i = 1; i <= m; i++) y[i] = w[P[m+i]];
345 xfree(w);
346 /* compute residual vector r = A*D*A'*y - h */
347 r = xcalloc(1+m, sizeof(double));
348 /* w := A'*y */
349 w = xcalloc(1+n, sizeof(double));
350 AT_by_vec(csa, y, w);
351 /* w := D*w */
352 for (j = 1; j <= n; j++) w[j] *= csa->D[j];
353 /* r := A*w */
354 A_by_vec(csa, w, r);
355 xfree(w);
356 /* r := r - h */
357 for (i = 1; i <= m; i++) r[i] -= h[i];
358 /* check for numeric stability */
359 for (i = 1; i <= m; i++)
360 { if (fabs(r[i]) / (1.0 + fabs(h[i])) > 1e-4)
361 { ret = 1;
362 break;
363 }
364 }
365 xfree(h);
366 xfree(r);
367 return ret;
368 }
370 /***********************************************************************
371 * solve_NS - solve Newtonian system
372 *
373 * This routine solves the Newtonian system:
374 *
375 * A*dx = p
376 *
377 * A'*dy + dz = q
378 *
379 * Z*dx + X*dz = r
380 *
381 * where X = diag(x[j]), Z = diag(z[j]), by reducing it to the normal
382 * equation system:
383 *
384 * (A*inv(Z)*X*A')*dy = A*inv(Z)*(X*q-r)+p
385 *
386 * (it is assumed that the matrix A*inv(Z)*X*A' has been factorized by
387 * the routine decomp_NE).
388 *
389 * Once vector dy has been computed the routine computes vectors dx and
390 * dz as follows:
391 *
392 * dx = inv(Z)*(X*(A'*dy-q)+r)
393 *
394 * dz = inv(X)*(r-Z*dx)
395 *
396 * The routine solve_NS returns the same code which was reported by the
397 * routine solve_NE (see above). */
399 static int solve_NS(struct csa *csa, double p[], double q[], double r[],
400 double dx[], double dy[], double dz[])
401 { int m = csa->m;
402 int n = csa->n;
403 double *x = csa->x;
404 double *z = csa->z;
405 int i, j, ret;
406 double *w = dx;
407 /* compute the vector of right-hand sides A*inv(Z)*(X*q-r)+p for
408 the normal equation system */
409 for (j = 1; j <= n; j++)
410 w[j] = (x[j] * q[j] - r[j]) / z[j];
411 A_by_vec(csa, w, dy);
412 for (i = 1; i <= m; i++) dy[i] += p[i];
413 /* solve the normal equation system to compute vector dy */
414 ret = solve_NE(csa, dy);
415 /* compute vectors dx and dz */
416 AT_by_vec(csa, dy, dx);
417 for (j = 1; j <= n; j++)
418 { dx[j] = (x[j] * (dx[j] - q[j]) + r[j]) / z[j];
419 dz[j] = (r[j] - z[j] * dx[j]) / x[j];
420 }
421 return ret;
422 }
424 /***********************************************************************
425 * initial_point - choose initial point using Mehrotra's heuristic
426 *
427 * This routine chooses a starting point using a heuristic proposed in
428 * the paper:
429 *
430 * S. Mehrotra. On the implementation of a primal-dual interior point
431 * method. SIAM J. on Optim., 2(4), pp. 575-601, 1992.
432 *
433 * The starting point x in the primal space is chosen as a solution of
434 * the following least squares problem:
435 *
436 * minimize ||x||
437 *
438 * subject to A*x = b
439 *
440 * which can be computed explicitly as follows:
441 *
442 * x = A'*inv(A*A')*b
443 *
444 * Similarly, the starting point (y, z) in the dual space is chosen as
445 * a solution of the following least squares problem:
446 *
447 * minimize ||z||
448 *
449 * subject to A'*y + z = c
450 *
451 * which can be computed explicitly as follows:
452 *
453 * y = inv(A*A')*A*c
454 *
455 * z = c - A'*y
456 *
457 * However, some components of the vectors x and z may be non-positive
458 * or close to zero, so the routine uses a Mehrotra's heuristic to find
459 * a more appropriate starting point. */
461 static void initial_point(struct csa *csa)
462 { int m = csa->m;
463 int n = csa->n;
464 double *b = csa->b;
465 double *c = csa->c;
466 double *x = csa->x;
467 double *y = csa->y;
468 double *z = csa->z;
469 double *D = csa->D;
470 int i, j;
471 double dp, dd, ex, ez, xz;
472 /* factorize A*A' */
473 for (j = 1; j <= n; j++) D[j] = 1.0;
474 decomp_NE(csa);
475 /* x~ = A'*inv(A*A')*b */
476 for (i = 1; i <= m; i++) y[i] = b[i];
477 solve_NE(csa, y);
478 AT_by_vec(csa, y, x);
479 /* y~ = inv(A*A')*A*c */
480 A_by_vec(csa, c, y);
481 solve_NE(csa, y);
482 /* z~ = c - A'*y~ */
483 AT_by_vec(csa, y,z);
484 for (j = 1; j <= n; j++) z[j] = c[j] - z[j];
485 /* use Mehrotra's heuristic in order to choose more appropriate
486 starting point with positive components of vectors x and z */
487 dp = dd = 0.0;
488 for (j = 1; j <= n; j++)
489 { if (dp < -1.5 * x[j]) dp = -1.5 * x[j];
490 if (dd < -1.5 * z[j]) dd = -1.5 * z[j];
491 }
492 /* note that b = 0 involves x = 0, and c = 0 involves y = 0 and
493 z = 0, so we need to be careful */
494 if (dp == 0.0) dp = 1.5;
495 if (dd == 0.0) dd = 1.5;
496 ex = ez = xz = 0.0;
497 for (j = 1; j <= n; j++)
498 { ex += (x[j] + dp);
499 ez += (z[j] + dd);
500 xz += (x[j] + dp) * (z[j] + dd);
501 }
502 dp += 0.5 * (xz / ez);
503 dd += 0.5 * (xz / ex);
504 for (j = 1; j <= n; j++)
505 { x[j] += dp;
506 z[j] += dd;
507 xassert(x[j] > 0.0 && z[j] > 0.0);
508 }
509 return;
510 }
512 /***********************************************************************
513 * basic_info - perform basic computations at the current point
514 *
515 * This routine computes the following quantities at the current point:
516 *
517 * 1) value of the objective function:
518 *
519 * F = c'*x + c[0]
520 *
521 * 2) relative primal infeasibility:
522 *
523 * rpi = ||A*x-b|| / (1+||b||)
524 *
525 * 3) relative dual infeasibility:
526 *
527 * rdi = ||A'*y+z-c|| / (1+||c||)
528 *
529 * 4) primal-dual gap (relative difference between the primal and the
530 * dual objective function values):
531 *
532 * gap = |c'*x-b'*y| / (1+|c'*x|)
533 *
534 * 5) merit function:
535 *
536 * phi = ||A*x-b|| / max(1,||b||) + ||A'*y+z-c|| / max(1,||c||) +
537 *
538 * + |c'*x-b'*y| / max(1,||b||,||c||)
539 *
540 * 6) duality measure:
541 *
542 * mu = x'*z / n
543 *
544 * 7) the ratio of infeasibility to mu:
545 *
546 * rmu = max(||A*x-b||,||A'*y+z-c||) / mu
547 *
548 * where ||*|| denotes euclidian norm, *' denotes transposition. */
550 static void basic_info(struct csa *csa)
551 { int m = csa->m;
552 int n = csa->n;
553 double *b = csa->b;
554 double *c = csa->c;
555 double *x = csa->x;
556 double *y = csa->y;
557 double *z = csa->z;
558 int i, j;
559 double norm1, bnorm, norm2, cnorm, cx, by, *work, temp;
560 /* compute value of the objective function */
561 temp = c[0];
562 for (j = 1; j <= n; j++) temp += c[j] * x[j];
563 csa->obj = temp;
564 /* norm1 = ||A*x-b|| */
565 work = xcalloc(1+m, sizeof(double));
566 A_by_vec(csa, x, work);
567 norm1 = 0.0;
568 for (i = 1; i <= m; i++)
569 norm1 += (work[i] - b[i]) * (work[i] - b[i]);
570 norm1 = sqrt(norm1);
571 xfree(work);
572 /* bnorm = ||b|| */
573 bnorm = 0.0;
574 for (i = 1; i <= m; i++) bnorm += b[i] * b[i];
575 bnorm = sqrt(bnorm);
576 /* compute relative primal infeasibility */
577 csa->rpi = norm1 / (1.0 + bnorm);
578 /* norm2 = ||A'*y+z-c|| */
579 work = xcalloc(1+n, sizeof(double));
580 AT_by_vec(csa, y, work);
581 norm2 = 0.0;
582 for (j = 1; j <= n; j++)
583 norm2 += (work[j] + z[j] - c[j]) * (work[j] + z[j] - c[j]);
584 norm2 = sqrt(norm2);
585 xfree(work);
586 /* cnorm = ||c|| */
587 cnorm = 0.0;
588 for (j = 1; j <= n; j++) cnorm += c[j] * c[j];
589 cnorm = sqrt(cnorm);
590 /* compute relative dual infeasibility */
591 csa->rdi = norm2 / (1.0 + cnorm);
592 /* by = b'*y */
593 by = 0.0;
594 for (i = 1; i <= m; i++) by += b[i] * y[i];
595 /* cx = c'*x */
596 cx = 0.0;
597 for (j = 1; j <= n; j++) cx += c[j] * x[j];
598 /* compute primal-dual gap */
599 csa->gap = fabs(cx - by) / (1.0 + fabs(cx));
600 /* compute merit function */
601 csa->phi = 0.0;
602 csa->phi += norm1 / (bnorm > 1.0 ? bnorm : 1.0);
603 csa->phi += norm2 / (cnorm > 1.0 ? cnorm : 1.0);
604 temp = 1.0;
605 if (temp < bnorm) temp = bnorm;
606 if (temp < cnorm) temp = cnorm;
607 csa->phi += fabs(cx - by) / temp;
608 /* compute duality measure */
609 temp = 0.0;
610 for (j = 1; j <= n; j++) temp += x[j] * z[j];
611 csa->mu = temp / (double)n;
612 /* compute the ratio of infeasibility to mu */
613 csa->rmu = (norm1 > norm2 ? norm1 : norm2) / csa->mu;
614 return;
615 }
617 /***********************************************************************
618 * make_step - compute next point using Mehrotra's technique
619 *
620 * This routine computes the next point using the predictor-corrector
621 * technique proposed in the paper:
622 *
623 * S. Mehrotra. On the implementation of a primal-dual interior point
624 * method. SIAM J. on Optim., 2(4), pp. 575-601, 1992.
625 *
626 * At first, the routine computes so called affine scaling (predictor)
627 * direction (dx_aff,dy_aff,dz_aff) which is a solution of the system:
628 *
629 * A*dx_aff = b - A*x
630 *
631 * A'*dy_aff + dz_aff = c - A'*y - z
632 *
633 * Z*dx_aff + X*dz_aff = - X*Z*e
634 *
635 * where (x,y,z) is the current point, X = diag(x[j]), Z = diag(z[j]),
636 * e = (1,...,1)'.
637 *
638 * Then, the routine computes the centering parameter sigma, using the
639 * following Mehrotra's heuristic:
640 *
641 * alfa_aff_p = inf{0 <= alfa <= 1 | x+alfa*dx_aff >= 0}
642 *
643 * alfa_aff_d = inf{0 <= alfa <= 1 | z+alfa*dz_aff >= 0}
644 *
645 * mu_aff = (x+alfa_aff_p*dx_aff)'*(z+alfa_aff_d*dz_aff)/n
646 *
647 * sigma = (mu_aff/mu)^3
648 *
649 * where alfa_aff_p is the maximal stepsize along the affine scaling
650 * direction in the primal space, alfa_aff_d is the maximal stepsize
651 * along the same direction in the dual space.
652 *
653 * After determining sigma the routine computes so called centering
654 * (corrector) direction (dx_cc,dy_cc,dz_cc) which is the solution of
655 * the system:
656 *
657 * A*dx_cc = 0
658 *
659 * A'*dy_cc + dz_cc = 0
660 *
661 * Z*dx_cc + X*dz_cc = sigma*mu*e - X*Z*e
662 *
663 * Finally, the routine computes the combined direction
664 *
665 * (dx,dy,dz) = (dx_aff,dy_aff,dz_aff) + (dx_cc,dy_cc,dz_cc)
666 *
667 * and determines maximal primal and dual stepsizes along the combined
668 * direction:
669 *
670 * alfa_max_p = inf{0 <= alfa <= 1 | x+alfa*dx >= 0}
671 *
672 * alfa_max_d = inf{0 <= alfa <= 1 | z+alfa*dz >= 0}
673 *
674 * In order to prevent the next point to be too close to the boundary
675 * of the positive ortant, the routine decreases maximal stepsizes:
676 *
677 * alfa_p = gamma_p * alfa_max_p
678 *
679 * alfa_d = gamma_d * alfa_max_d
680 *
681 * where gamma_p and gamma_d are scaling factors, and computes the next
682 * point:
683 *
684 * x_new = x + alfa_p * dx
685 *
686 * y_new = y + alfa_d * dy
687 *
688 * z_new = z + alfa_d * dz
689 *
690 * which becomes the current point on the next iteration. */
692 static int make_step(struct csa *csa)
693 { int m = csa->m;
694 int n = csa->n;
695 double *b = csa->b;
696 double *c = csa->c;
697 double *x = csa->x;
698 double *y = csa->y;
699 double *z = csa->z;
700 double *dx_aff = csa->dx_aff;
701 double *dy_aff = csa->dy_aff;
702 double *dz_aff = csa->dz_aff;
703 double *dx_cc = csa->dx_cc;
704 double *dy_cc = csa->dy_cc;
705 double *dz_cc = csa->dz_cc;
706 double *dx = csa->dx;
707 double *dy = csa->dy;
708 double *dz = csa->dz;
709 int i, j, ret = 0;
710 double temp, gamma_p, gamma_d, *p, *q, *r;
711 /* allocate working arrays */
712 p = xcalloc(1+m, sizeof(double));
713 q = xcalloc(1+n, sizeof(double));
714 r = xcalloc(1+n, sizeof(double));
715 /* p = b - A*x */
716 A_by_vec(csa, x, p);
717 for (i = 1; i <= m; i++) p[i] = b[i] - p[i];
718 /* q = c - A'*y - z */
719 AT_by_vec(csa, y,q);
720 for (j = 1; j <= n; j++) q[j] = c[j] - q[j] - z[j];
721 /* r = - X * Z * e */
722 for (j = 1; j <= n; j++) r[j] = - x[j] * z[j];
723 /* solve the first Newtonian system */
724 if (solve_NS(csa, p, q, r, dx_aff, dy_aff, dz_aff))
725 { ret = 1;
726 goto done;
727 }
728 /* alfa_aff_p = inf{0 <= alfa <= 1 | x + alfa*dx_aff >= 0} */
729 /* alfa_aff_d = inf{0 <= alfa <= 1 | z + alfa*dz_aff >= 0} */
730 csa->alfa_aff_p = csa->alfa_aff_d = 1.0;
731 for (j = 1; j <= n; j++)
732 { if (dx_aff[j] < 0.0)
733 { temp = - x[j] / dx_aff[j];
734 if (csa->alfa_aff_p > temp) csa->alfa_aff_p = temp;
735 }
736 if (dz_aff[j] < 0.0)
737 { temp = - z[j] / dz_aff[j];
738 if (csa->alfa_aff_d > temp) csa->alfa_aff_d = temp;
739 }
740 }
741 /* mu_aff = (x+alfa_aff_p*dx_aff)' * (z+alfa_aff_d*dz_aff) / n */
742 temp = 0.0;
743 for (j = 1; j <= n; j++)
744 temp += (x[j] + csa->alfa_aff_p * dx_aff[j]) *
745 (z[j] + csa->alfa_aff_d * dz_aff[j]);
746 csa->mu_aff = temp / (double)n;
747 /* sigma = (mu_aff/mu)^3 */
748 temp = csa->mu_aff / csa->mu;
749 csa->sigma = temp * temp * temp;
750 /* p = 0 */
751 for (i = 1; i <= m; i++) p[i] = 0.0;
752 /* q = 0 */
753 for (j = 1; j <= n; j++) q[j] = 0.0;
754 /* r = sigma * mu * e - X * Z * e */
755 for (j = 1; j <= n; j++)
756 r[j] = csa->sigma * csa->mu - dx_aff[j] * dz_aff[j];
757 /* solve the second Newtonian system with the same coefficients
758 but with altered right-hand sides */
759 if (solve_NS(csa, p, q, r, dx_cc, dy_cc, dz_cc))
760 { ret = 1;
761 goto done;
762 }
763 /* (dx,dy,dz) = (dx_aff,dy_aff,dz_aff) + (dx_cc,dy_cc,dz_cc) */
764 for (j = 1; j <= n; j++) dx[j] = dx_aff[j] + dx_cc[j];
765 for (i = 1; i <= m; i++) dy[i] = dy_aff[i] + dy_cc[i];
766 for (j = 1; j <= n; j++) dz[j] = dz_aff[j] + dz_cc[j];
767 /* alfa_max_p = inf{0 <= alfa <= 1 | x + alfa*dx >= 0} */
768 /* alfa_max_d = inf{0 <= alfa <= 1 | z + alfa*dz >= 0} */
769 csa->alfa_max_p = csa->alfa_max_d = 1.0;
770 for (j = 1; j <= n; j++)
771 { if (dx[j] < 0.0)
772 { temp = - x[j] / dx[j];
773 if (csa->alfa_max_p > temp) csa->alfa_max_p = temp;
774 }
775 if (dz[j] < 0.0)
776 { temp = - z[j] / dz[j];
777 if (csa->alfa_max_d > temp) csa->alfa_max_d = temp;
778 }
779 }
780 /* determine scale factors (not implemented yet) */
781 gamma_p = 0.90;
782 gamma_d = 0.90;
783 /* compute the next point */
784 for (j = 1; j <= n; j++)
785 { x[j] += gamma_p * csa->alfa_max_p * dx[j];
786 xassert(x[j] > 0.0);
787 }
788 for (i = 1; i <= m; i++)
789 y[i] += gamma_d * csa->alfa_max_d * dy[i];
790 for (j = 1; j <= n; j++)
791 { z[j] += gamma_d * csa->alfa_max_d * dz[j];
792 xassert(z[j] > 0.0);
793 }
794 done: /* free working arrays */
795 xfree(p);
796 xfree(q);
797 xfree(r);
798 return ret;
799 }
801 /***********************************************************************
802 * terminate - deallocate common storage area
803 *
804 * This routine frees all memory allocated to the common storage area
805 * used by interior-point method routines. */
807 static void terminate(struct csa *csa)
808 { xfree(csa->D);
809 xfree(csa->P);
810 xfree(csa->S_ptr);
811 xfree(csa->S_ind);
812 xfree(csa->S_val);
813 xfree(csa->S_diag);
814 xfree(csa->U_ptr);
815 xfree(csa->U_ind);
816 xfree(csa->U_val);
817 xfree(csa->U_diag);
818 xfree(csa->phi_min);
819 xfree(csa->best_x);
820 xfree(csa->best_y);
821 xfree(csa->best_z);
822 xfree(csa->dx_aff);
823 xfree(csa->dy_aff);
824 xfree(csa->dz_aff);
825 xfree(csa->dx_cc);
826 xfree(csa->dy_cc);
827 xfree(csa->dz_cc);
828 return;
829 }
831 /***********************************************************************
832 * ipm_main - main interior-point method routine
833 *
834 * This is a main routine of the primal-dual interior-point method.
835 *
836 * The routine ipm_main returns one of the following codes:
837 *
838 * 0 - optimal solution found;
839 * 1 - problem has no feasible (primal or dual) solution;
840 * 2 - no convergence;
841 * 3 - iteration limit exceeded;
842 * 4 - numeric instability on solving Newtonian system.
843 *
844 * In case of non-zero return code the routine returns the best point,
845 * which has been reached during optimization. */
847 static int ipm_main(struct csa *csa)
848 { int m = csa->m;
849 int n = csa->n;
850 int i, j, status;
851 double temp;
852 /* choose initial point using Mehrotra's heuristic */
853 if (csa->parm->msg_lev >= GLP_MSG_ALL)
854 xprintf("Guessing initial point...\n");
855 initial_point(csa);
856 /* main loop starts here */
857 if (csa->parm->msg_lev >= GLP_MSG_ALL)
858 xprintf("Optimization begins...\n");
859 for (;;)
860 { /* perform basic computations at the current point */
861 basic_info(csa);
862 /* save initial value of rmu */
863 if (csa->iter == 0) csa->rmu0 = csa->rmu;
864 /* accumulate values of min(phi[k]) and save the best point */
865 xassert(csa->iter <= ITER_MAX);
866 if (csa->iter == 0 || csa->phi_min[csa->iter-1] > csa->phi)
867 { csa->phi_min[csa->iter] = csa->phi;
868 csa->best_iter = csa->iter;
869 for (j = 1; j <= n; j++) csa->best_x[j] = csa->x[j];
870 for (i = 1; i <= m; i++) csa->best_y[i] = csa->y[i];
871 for (j = 1; j <= n; j++) csa->best_z[j] = csa->z[j];
872 csa->best_obj = csa->obj;
873 }
874 else
875 csa->phi_min[csa->iter] = csa->phi_min[csa->iter-1];
876 /* display information at the current point */
877 if (csa->parm->msg_lev >= GLP_MSG_ON)
878 xprintf("%3d: obj = %17.9e; rpi = %8.1e; rdi = %8.1e; gap ="
879 " %8.1e\n", csa->iter, csa->obj, csa->rpi, csa->rdi,
880 csa->gap);
881 /* check if the current point is optimal */
882 if (csa->rpi < 1e-8 && csa->rdi < 1e-8 && csa->gap < 1e-8)
883 { if (csa->parm->msg_lev >= GLP_MSG_ALL)
884 xprintf("OPTIMAL SOLUTION FOUND\n");
885 status = 0;
886 break;
887 }
888 /* check if the problem has no feasible solution */
889 temp = 1e5 * csa->phi_min[csa->iter];
890 if (temp < 1e-8) temp = 1e-8;
891 if (csa->phi >= temp)
892 { if (csa->parm->msg_lev >= GLP_MSG_ALL)
893 xprintf("PROBLEM HAS NO FEASIBLE PRIMAL/DUAL SOLUTION\n")
894 ;
895 status = 1;
896 break;
897 }
898 /* check for very slow convergence or divergence */
899 if (((csa->rpi >= 1e-8 || csa->rdi >= 1e-8) && csa->rmu /
900 csa->rmu0 >= 1e6) ||
901 (csa->iter >= 30 && csa->phi_min[csa->iter] >= 0.5 *
902 csa->phi_min[csa->iter - 30]))
903 { if (csa->parm->msg_lev >= GLP_MSG_ALL)
904 xprintf("NO CONVERGENCE; SEARCH TERMINATED\n");
905 status = 2;
906 break;
907 }
908 /* check for maximal number of iterations */
909 if (csa->iter == ITER_MAX)
910 { if (csa->parm->msg_lev >= GLP_MSG_ALL)
911 xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n");
912 status = 3;
913 break;
914 }
915 /* start the next iteration */
916 csa->iter++;
917 /* factorize normal equation system */
918 for (j = 1; j <= n; j++) csa->D[j] = csa->x[j] / csa->z[j];
919 decomp_NE(csa);
920 /* compute the next point using Mehrotra's predictor-corrector
921 technique */
922 if (make_step(csa))
923 { if (csa->parm->msg_lev >= GLP_MSG_ALL)
924 xprintf("NUMERIC INSTABILITY; SEARCH TERMINATED\n");
925 status = 4;
926 break;
927 }
928 }
929 /* restore the best point */
930 if (status != 0)
931 { for (j = 1; j <= n; j++) csa->x[j] = csa->best_x[j];
932 for (i = 1; i <= m; i++) csa->y[i] = csa->best_y[i];
933 for (j = 1; j <= n; j++) csa->z[j] = csa->best_z[j];
934 if (csa->parm->msg_lev >= GLP_MSG_ALL)
935 xprintf("Best point %17.9e was reached on iteration %d\n",
936 csa->best_obj, csa->best_iter);
937 }
938 /* return to the calling program */
939 return status;
940 }
942 /***********************************************************************
943 * NAME
944 *
945 * ipm_solve - core LP solver based on the interior-point method
946 *
947 * SYNOPSIS
948 *
949 * #include "glpipm.h"
950 * int ipm_solve(glp_prob *P, const glp_iptcp *parm);
951 *
952 * DESCRIPTION
953 *
954 * The routine ipm_solve is a core LP solver based on the primal-dual
955 * interior-point method.
956 *
957 * The routine assumes the following standard formulation of LP problem
958 * to be solved:
959 *
960 * minimize
961 *
962 * F = c[0] + c[1]*x[1] + c[2]*x[2] + ... + c[n]*x[n]
963 *
964 * subject to linear constraints
965 *
966 * a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1]
967 *
968 * a[2,1]*x[1] + a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2]
969 *
970 * . . . . . .
971 *
972 * a[m,1]*x[1] + a[m,2]*x[2] + ... + a[m,n]*x[n] = b[m]
973 *
974 * and non-negative variables
975 *
976 * x[1] >= 0, x[2] >= 0, ..., x[n] >= 0
977 *
978 * where:
979 * F is the objective function;
980 * x[1], ..., x[n] are (structural) variables;
981 * c[0] is a constant term of the objective function;
982 * c[1], ..., c[n] are objective coefficients;
983 * a[1,1], ..., a[m,n] are constraint coefficients;
984 * b[1], ..., b[n] are right-hand sides.
985 *
986 * The solution is three vectors x, y, and z, which are stored by the
987 * routine in the arrays x, y, and z, respectively. These vectors
988 * correspond to the best primal-dual point found during optimization.
989 * They are approximate solution of the following system (which is the
990 * Karush-Kuhn-Tucker optimality conditions):
991 *
992 * A*x = b (primal feasibility condition)
993 *
994 * A'*y + z = c (dual feasibility condition)
995 *
996 * x'*z = 0 (primal-dual complementarity condition)
997 *
998 * x >= 0, z >= 0 (non-negativity condition)
999 *
1000 * where:
1001 * x[1], ..., x[n] are primal (structural) variables;
1002 * y[1], ..., y[m] are dual variables (Lagrange multipliers) for
1003 * equality constraints;
1004 * z[1], ..., z[n] are dual variables (Lagrange multipliers) for
1005 * non-negativity constraints.
1006 *
1007 * RETURNS
1008 *
1009 * 0 LP has been successfully solved.
1010 *
1011 * GLP_ENOCVG
1012 * No convergence.
1013 *
1014 * GLP_EITLIM
1015 * Iteration limit exceeded.
1016 *
1017 * GLP_EINSTAB
1018 * Numeric instability on solving Newtonian system.
1019 *
1020 * In case of non-zero return code the routine returns the best point,
1021 * which has been reached during optimization. */
1023 int ipm_solve(glp_prob *P, const glp_iptcp *parm)
1024 { struct csa _dsa, *csa = &_dsa;
1025 int m = P->m;
1026 int n = P->n;
1027 int nnz = P->nnz;
1028 GLPROW *row;
1029 GLPCOL *col;
1030 GLPAIJ *aij;
1031 int i, j, loc, ret, *A_ind, *A_ptr;
1032 double dir, *A_val, *b, *c, *x, *y, *z;
1033 xassert(m > 0);
1034 xassert(n > 0);
1035 /* allocate working arrays */
1036 A_ptr = xcalloc(1+m+1, sizeof(int));
1037 A_ind = xcalloc(1+nnz, sizeof(int));
1038 A_val = xcalloc(1+nnz, sizeof(double));
1039 b = xcalloc(1+m, sizeof(double));
1040 c = xcalloc(1+n, sizeof(double));
1041 x = xcalloc(1+n, sizeof(double));
1042 y = xcalloc(1+m, sizeof(double));
1043 z = xcalloc(1+n, sizeof(double));
1044 /* prepare rows and constraint coefficients */
1045 loc = 1;
1046 for (i = 1; i <= m; i++)
1047 { row = P->row[i];
1048 xassert(row->type == GLP_FX);
1049 b[i] = row->lb * row->rii;
1050 A_ptr[i] = loc;
1051 for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1052 { A_ind[loc] = aij->col->j;
1053 A_val[loc] = row->rii * aij->val * aij->col->sjj;
1054 loc++;
1055 }
1056 }
1057 A_ptr[m+1] = loc;
1058 xassert(loc-1 == nnz);
1059 /* prepare columns and objective coefficients */
1060 if (P->dir == GLP_MIN)
1061 dir = +1.0;
1062 else if (P->dir == GLP_MAX)
1063 dir = -1.0;
1064 else
1065 xassert(P != P);
1066 c[0] = dir * P->c0;
1067 for (j = 1; j <= n; j++)
1068 { col = P->col[j];
1069 xassert(col->type == GLP_LO && col->lb == 0.0);
1070 c[j] = dir * col->coef * col->sjj;
1071 }
1072 /* allocate and initialize the common storage area */
1073 csa->m = m;
1074 csa->n = n;
1075 csa->A_ptr = A_ptr;
1076 csa->A_ind = A_ind;
1077 csa->A_val = A_val;
1078 csa->b = b;
1079 csa->c = c;
1080 csa->x = x;
1081 csa->y = y;
1082 csa->z = z;
1083 csa->parm = parm;
1084 initialize(csa);
1085 /* solve LP with the interior-point method */
1086 ret = ipm_main(csa);
1087 /* deallocate the common storage area */
1088 terminate(csa);
1089 /* determine solution status */
1090 if (ret == 0)
1091 { /* optimal solution found */
1092 P->ipt_stat = GLP_OPT;
1093 ret = 0;
1094 }
1095 else if (ret == 1)
1096 { /* problem has no feasible (primal or dual) solution */
1097 P->ipt_stat = GLP_NOFEAS;
1098 ret = 0;
1099 }
1100 else if (ret == 2)
1101 { /* no convergence */
1102 P->ipt_stat = GLP_INFEAS;
1103 ret = GLP_ENOCVG;
1104 }
1105 else if (ret == 3)
1106 { /* iteration limit exceeded */
1107 P->ipt_stat = GLP_INFEAS;
1108 ret = GLP_EITLIM;
1109 }
1110 else if (ret == 4)
1111 { /* numeric instability on solving Newtonian system */
1112 P->ipt_stat = GLP_INFEAS;
1113 ret = GLP_EINSTAB;
1114 }
1115 else
1116 xassert(ret != ret);
1117 /* store row solution components */
1118 for (i = 1; i <= m; i++)
1119 { row = P->row[i];
1120 row->pval = row->lb;
1121 row->dval = dir * y[i] * row->rii;
1122 }
1123 /* store column solution components */
1124 P->ipt_obj = P->c0;
1125 for (j = 1; j <= n; j++)
1126 { col = P->col[j];
1127 col->pval = x[j] * col->sjj;
1128 col->dval = dir * z[j] / col->sjj;
1129 P->ipt_obj += col->coef * col->pval;
1130 }
1131 /* free working arrays */
1132 xfree(A_ptr);
1133 xfree(A_ind);
1134 xfree(A_val);
1135 xfree(b);
1136 xfree(c);
1137 xfree(x);
1138 xfree(y);
1139 xfree(z);
1140 return ret;
1141 }
1143 /* eof */