lemon-project-template-glpk

annotate deps/glpk/src/glpipm.c @ 9:33de93886c88

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