1.1 --- /dev/null Thu Jan 01 00:00:00 1970 +0000
1.2 +++ b/src/glpipm.c Mon Dec 06 13:09:21 2010 +0100
1.3 @@ -0,0 +1,1143 @@
1.4 +/* glpipm.c */
1.5 +
1.6 +/***********************************************************************
1.7 +* This code is part of GLPK (GNU Linear Programming Kit).
1.8 +*
1.9 +* Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
1.10 +* 2009, 2010 Andrew Makhorin, Department for Applied Informatics,
1.11 +* Moscow Aviation Institute, Moscow, Russia. All rights reserved.
1.12 +* E-mail: <mao@gnu.org>.
1.13 +*
1.14 +* GLPK is free software: you can redistribute it and/or modify it
1.15 +* under the terms of the GNU General Public License as published by
1.16 +* the Free Software Foundation, either version 3 of the License, or
1.17 +* (at your option) any later version.
1.18 +*
1.19 +* GLPK is distributed in the hope that it will be useful, but WITHOUT
1.20 +* ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
1.21 +* or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
1.22 +* License for more details.
1.23 +*
1.24 +* You should have received a copy of the GNU General Public License
1.25 +* along with GLPK. If not, see <http://www.gnu.org/licenses/>.
1.26 +***********************************************************************/
1.27 +
1.28 +#include "glpipm.h"
1.29 +#include "glpmat.h"
1.30 +
1.31 +#define ITER_MAX 100
1.32 +/* maximal number of iterations */
1.33 +
1.34 +struct csa
1.35 +{ /* common storage area */
1.36 + /*--------------------------------------------------------------*/
1.37 + /* LP data */
1.38 + int m;
1.39 + /* number of rows (equality constraints) */
1.40 + int n;
1.41 + /* number of columns (structural variables) */
1.42 + int *A_ptr; /* int A_ptr[1+m+1]; */
1.43 + int *A_ind; /* int A_ind[A_ptr[m+1]]; */
1.44 + double *A_val; /* double A_val[A_ptr[m+1]]; */
1.45 + /* mxn-matrix A in storage-by-rows format */
1.46 + double *b; /* double b[1+m]; */
1.47 + /* m-vector b of right-hand sides */
1.48 + double *c; /* double c[1+n]; */
1.49 + /* n-vector c of objective coefficients; c[0] is constant term of
1.50 + the objective function */
1.51 + /*--------------------------------------------------------------*/
1.52 + /* LP solution */
1.53 + double *x; /* double x[1+n]; */
1.54 + double *y; /* double y[1+m]; */
1.55 + double *z; /* double z[1+n]; */
1.56 + /* current point in primal-dual space; the best point on exit */
1.57 + /*--------------------------------------------------------------*/
1.58 + /* control parameters */
1.59 + const glp_iptcp *parm;
1.60 + /*--------------------------------------------------------------*/
1.61 + /* working arrays and variables */
1.62 + double *D; /* double D[1+n]; */
1.63 + /* diagonal nxn-matrix D = X*inv(Z), where X = diag(x[j]) and
1.64 + Z = diag(z[j]) */
1.65 + int *P; /* int P[1+m+m]; */
1.66 + /* permutation mxm-matrix P used to minimize fill-in in Cholesky
1.67 + factorization */
1.68 + int *S_ptr; /* int S_ptr[1+m+1]; */
1.69 + int *S_ind; /* int S_ind[S_ptr[m+1]]; */
1.70 + double *S_val; /* double S_val[S_ptr[m+1]]; */
1.71 + double *S_diag; /* double S_diag[1+m]; */
1.72 + /* symmetric mxm-matrix S = P*A*D*A'*P' whose upper triangular
1.73 + part without diagonal elements is stored in S_ptr, S_ind, and
1.74 + S_val in storage-by-rows format, diagonal elements are stored
1.75 + in S_diag */
1.76 + int *U_ptr; /* int U_ptr[1+m+1]; */
1.77 + int *U_ind; /* int U_ind[U_ptr[m+1]]; */
1.78 + double *U_val; /* double U_val[U_ptr[m+1]]; */
1.79 + double *U_diag; /* double U_diag[1+m]; */
1.80 + /* upper triangular mxm-matrix U defining Cholesky factorization
1.81 + S = U'*U; its non-diagonal elements are stored in U_ptr, U_ind,
1.82 + U_val in storage-by-rows format, diagonal elements are stored
1.83 + in U_diag */
1.84 + int iter;
1.85 + /* iteration number (0, 1, 2, ...); iter = 0 corresponds to the
1.86 + initial point */
1.87 + double obj;
1.88 + /* current value of the objective function */
1.89 + double rpi;
1.90 + /* relative primal infeasibility rpi = ||A*x-b||/(1+||b||) */
1.91 + double rdi;
1.92 + /* relative dual infeasibility rdi = ||A'*y+z-c||/(1+||c||) */
1.93 + double gap;
1.94 + /* primal-dual gap = |c'*x-b'*y|/(1+|c'*x|) which is a relative
1.95 + difference between primal and dual objective functions */
1.96 + double phi;
1.97 + /* merit function phi = ||A*x-b||/max(1,||b||) +
1.98 + + ||A'*y+z-c||/max(1,||c||) +
1.99 + + |c'*x-b'*y|/max(1,||b||,||c||) */
1.100 + double mu;
1.101 + /* duality measure mu = x'*z/n (used as barrier parameter) */
1.102 + double rmu;
1.103 + /* rmu = max(||A*x-b||,||A'*y+z-c||)/mu */
1.104 + double rmu0;
1.105 + /* the initial value of rmu on iteration 0 */
1.106 + double *phi_min; /* double phi_min[1+ITER_MAX]; */
1.107 + /* phi_min[k] = min(phi[k]), where phi[k] is the value of phi on
1.108 + k-th iteration, 0 <= k <= iter */
1.109 + int best_iter;
1.110 + /* iteration number, on which the value of phi reached its best
1.111 + (minimal) value */
1.112 + double *best_x; /* double best_x[1+n]; */
1.113 + double *best_y; /* double best_y[1+m]; */
1.114 + double *best_z; /* double best_z[1+n]; */
1.115 + /* best point (in the sense of the merit function phi) which has
1.116 + been reached on iteration iter_best */
1.117 + double best_obj;
1.118 + /* objective value at the best point */
1.119 + double *dx_aff; /* double dx_aff[1+n]; */
1.120 + double *dy_aff; /* double dy_aff[1+m]; */
1.121 + double *dz_aff; /* double dz_aff[1+n]; */
1.122 + /* affine scaling direction */
1.123 + double alfa_aff_p, alfa_aff_d;
1.124 + /* maximal primal and dual stepsizes in affine scaling direction,
1.125 + on which x and z are still non-negative */
1.126 + double mu_aff;
1.127 + /* duality measure mu_aff = x_aff'*z_aff/n in the boundary point
1.128 + x_aff' = x+alfa_aff_p*dx_aff, z_aff' = z+alfa_aff_d*dz_aff */
1.129 + double sigma;
1.130 + /* Mehrotra's heuristic parameter (0 <= sigma <= 1) */
1.131 + double *dx_cc; /* double dx_cc[1+n]; */
1.132 + double *dy_cc; /* double dy_cc[1+m]; */
1.133 + double *dz_cc; /* double dz_cc[1+n]; */
1.134 + /* centering corrector direction */
1.135 + double *dx; /* double dx[1+n]; */
1.136 + double *dy; /* double dy[1+m]; */
1.137 + double *dz; /* double dz[1+n]; */
1.138 + /* final combined direction dx = dx_aff+dx_cc, dy = dy_aff+dy_cc,
1.139 + dz = dz_aff+dz_cc */
1.140 + double alfa_max_p;
1.141 + double alfa_max_d;
1.142 + /* maximal primal and dual stepsizes in combined direction, on
1.143 + which x and z are still non-negative */
1.144 +};
1.145 +
1.146 +/***********************************************************************
1.147 +* initialize - allocate and initialize common storage area
1.148 +*
1.149 +* This routine allocates and initializes the common storage area (CSA)
1.150 +* used by interior-point method routines. */
1.151 +
1.152 +static void initialize(struct csa *csa)
1.153 +{ int m = csa->m;
1.154 + int n = csa->n;
1.155 + int i;
1.156 + if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.157 + xprintf("Matrix A has %d non-zeros\n", csa->A_ptr[m+1]-1);
1.158 + csa->D = xcalloc(1+n, sizeof(double));
1.159 + /* P := I */
1.160 + csa->P = xcalloc(1+m+m, sizeof(int));
1.161 + for (i = 1; i <= m; i++) csa->P[i] = csa->P[m+i] = i;
1.162 + /* S := A*A', symbolically */
1.163 + csa->S_ptr = xcalloc(1+m+1, sizeof(int));
1.164 + csa->S_ind = adat_symbolic(m, n, csa->P, csa->A_ptr, csa->A_ind,
1.165 + csa->S_ptr);
1.166 + if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.167 + xprintf("Matrix S = A*A' has %d non-zeros (upper triangle)\n",
1.168 + csa->S_ptr[m+1]-1 + m);
1.169 + /* determine P using specified ordering algorithm */
1.170 + if (csa->parm->ord_alg == GLP_ORD_NONE)
1.171 + { if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.172 + xprintf("Original ordering is being used\n");
1.173 + for (i = 1; i <= m; i++)
1.174 + csa->P[i] = csa->P[m+i] = i;
1.175 + }
1.176 + else if (csa->parm->ord_alg == GLP_ORD_QMD)
1.177 + { if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.178 + xprintf("Minimum degree ordering (QMD)...\n");
1.179 + min_degree(m, csa->S_ptr, csa->S_ind, csa->P);
1.180 + }
1.181 + else if (csa->parm->ord_alg == GLP_ORD_AMD)
1.182 + { if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.183 + xprintf("Approximate minimum degree ordering (AMD)...\n");
1.184 + amd_order1(m, csa->S_ptr, csa->S_ind, csa->P);
1.185 + }
1.186 + else if (csa->parm->ord_alg == GLP_ORD_SYMAMD)
1.187 + { if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.188 + xprintf("Approximate minimum degree ordering (SYMAMD)...\n")
1.189 + ;
1.190 + symamd_ord(m, csa->S_ptr, csa->S_ind, csa->P);
1.191 + }
1.192 + else
1.193 + xassert(csa != csa);
1.194 + /* S := P*A*A'*P', symbolically */
1.195 + xfree(csa->S_ind);
1.196 + csa->S_ind = adat_symbolic(m, n, csa->P, csa->A_ptr, csa->A_ind,
1.197 + csa->S_ptr);
1.198 + csa->S_val = xcalloc(csa->S_ptr[m+1], sizeof(double));
1.199 + csa->S_diag = xcalloc(1+m, sizeof(double));
1.200 + /* compute Cholesky factorization S = U'*U, symbolically */
1.201 + if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.202 + xprintf("Computing Cholesky factorization S = L*L'...\n");
1.203 + csa->U_ptr = xcalloc(1+m+1, sizeof(int));
1.204 + csa->U_ind = chol_symbolic(m, csa->S_ptr, csa->S_ind, csa->U_ptr);
1.205 + if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.206 + xprintf("Matrix L has %d non-zeros\n", csa->U_ptr[m+1]-1 + m);
1.207 + csa->U_val = xcalloc(csa->U_ptr[m+1], sizeof(double));
1.208 + csa->U_diag = xcalloc(1+m, sizeof(double));
1.209 + csa->iter = 0;
1.210 + csa->obj = 0.0;
1.211 + csa->rpi = 0.0;
1.212 + csa->rdi = 0.0;
1.213 + csa->gap = 0.0;
1.214 + csa->phi = 0.0;
1.215 + csa->mu = 0.0;
1.216 + csa->rmu = 0.0;
1.217 + csa->rmu0 = 0.0;
1.218 + csa->phi_min = xcalloc(1+ITER_MAX, sizeof(double));
1.219 + csa->best_iter = 0;
1.220 + csa->best_x = xcalloc(1+n, sizeof(double));
1.221 + csa->best_y = xcalloc(1+m, sizeof(double));
1.222 + csa->best_z = xcalloc(1+n, sizeof(double));
1.223 + csa->best_obj = 0.0;
1.224 + csa->dx_aff = xcalloc(1+n, sizeof(double));
1.225 + csa->dy_aff = xcalloc(1+m, sizeof(double));
1.226 + csa->dz_aff = xcalloc(1+n, sizeof(double));
1.227 + csa->alfa_aff_p = 0.0;
1.228 + csa->alfa_aff_d = 0.0;
1.229 + csa->mu_aff = 0.0;
1.230 + csa->sigma = 0.0;
1.231 + csa->dx_cc = xcalloc(1+n, sizeof(double));
1.232 + csa->dy_cc = xcalloc(1+m, sizeof(double));
1.233 + csa->dz_cc = xcalloc(1+n, sizeof(double));
1.234 + csa->dx = csa->dx_aff;
1.235 + csa->dy = csa->dy_aff;
1.236 + csa->dz = csa->dz_aff;
1.237 + csa->alfa_max_p = 0.0;
1.238 + csa->alfa_max_d = 0.0;
1.239 + return;
1.240 +}
1.241 +
1.242 +/***********************************************************************
1.243 +* A_by_vec - compute y = A*x
1.244 +*
1.245 +* This routine computes matrix-vector product y = A*x, where A is the
1.246 +* constraint matrix. */
1.247 +
1.248 +static void A_by_vec(struct csa *csa, double x[], double y[])
1.249 +{ /* compute y = A*x */
1.250 + int m = csa->m;
1.251 + int *A_ptr = csa->A_ptr;
1.252 + int *A_ind = csa->A_ind;
1.253 + double *A_val = csa->A_val;
1.254 + int i, t, beg, end;
1.255 + double temp;
1.256 + for (i = 1; i <= m; i++)
1.257 + { temp = 0.0;
1.258 + beg = A_ptr[i], end = A_ptr[i+1];
1.259 + for (t = beg; t < end; t++) temp += A_val[t] * x[A_ind[t]];
1.260 + y[i] = temp;
1.261 + }
1.262 + return;
1.263 +}
1.264 +
1.265 +/***********************************************************************
1.266 +* AT_by_vec - compute y = A'*x
1.267 +*
1.268 +* This routine computes matrix-vector product y = A'*x, where A' is a
1.269 +* matrix transposed to the constraint matrix A. */
1.270 +
1.271 +static void AT_by_vec(struct csa *csa, double x[], double y[])
1.272 +{ /* compute y = A'*x, where A' is transposed to A */
1.273 + int m = csa->m;
1.274 + int n = csa->n;
1.275 + int *A_ptr = csa->A_ptr;
1.276 + int *A_ind = csa->A_ind;
1.277 + double *A_val = csa->A_val;
1.278 + int i, j, t, beg, end;
1.279 + double temp;
1.280 + for (j = 1; j <= n; j++) y[j] = 0.0;
1.281 + for (i = 1; i <= m; i++)
1.282 + { temp = x[i];
1.283 + if (temp == 0.0) continue;
1.284 + beg = A_ptr[i], end = A_ptr[i+1];
1.285 + for (t = beg; t < end; t++) y[A_ind[t]] += A_val[t] * temp;
1.286 + }
1.287 + return;
1.288 +}
1.289 +
1.290 +/***********************************************************************
1.291 +* decomp_NE - numeric factorization of matrix S = P*A*D*A'*P'
1.292 +*
1.293 +* This routine implements numeric phase of Cholesky factorization of
1.294 +* the matrix S = P*A*D*A'*P', which is a permuted matrix of the normal
1.295 +* equation system. Matrix D is assumed to be already computed. */
1.296 +
1.297 +static void decomp_NE(struct csa *csa)
1.298 +{ adat_numeric(csa->m, csa->n, csa->P, csa->A_ptr, csa->A_ind,
1.299 + csa->A_val, csa->D, csa->S_ptr, csa->S_ind, csa->S_val,
1.300 + csa->S_diag);
1.301 + chol_numeric(csa->m, csa->S_ptr, csa->S_ind, csa->S_val,
1.302 + csa->S_diag, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag);
1.303 + return;
1.304 +}
1.305 +
1.306 +/***********************************************************************
1.307 +* solve_NE - solve normal equation system
1.308 +*
1.309 +* This routine solves the normal equation system:
1.310 +*
1.311 +* A*D*A'*y = h.
1.312 +*
1.313 +* It is assumed that the matrix A*D*A' has been previously factorized
1.314 +* by the routine decomp_NE.
1.315 +*
1.316 +* On entry the array y contains the vector of right-hand sides h. On
1.317 +* exit this array contains the computed vector of unknowns y.
1.318 +*
1.319 +* Once the vector y has been computed the routine checks for numeric
1.320 +* stability. If the residual vector:
1.321 +*
1.322 +* r = A*D*A'*y - h
1.323 +*
1.324 +* is relatively small, the routine returns zero, otherwise non-zero is
1.325 +* returned. */
1.326 +
1.327 +static int solve_NE(struct csa *csa, double y[])
1.328 +{ int m = csa->m;
1.329 + int n = csa->n;
1.330 + int *P = csa->P;
1.331 + int i, j, ret = 0;
1.332 + double *h, *r, *w;
1.333 + /* save vector of right-hand sides h */
1.334 + h = xcalloc(1+m, sizeof(double));
1.335 + for (i = 1; i <= m; i++) h[i] = y[i];
1.336 + /* solve normal equation system (A*D*A')*y = h */
1.337 + /* since S = P*A*D*A'*P' = U'*U, then A*D*A' = P'*U'*U*P, so we
1.338 + have inv(A*D*A') = P'*inv(U)*inv(U')*P */
1.339 + /* w := P*h */
1.340 + w = xcalloc(1+m, sizeof(double));
1.341 + for (i = 1; i <= m; i++) w[i] = y[P[i]];
1.342 + /* w := inv(U')*w */
1.343 + ut_solve(m, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag, w);
1.344 + /* w := inv(U)*w */
1.345 + u_solve(m, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag, w);
1.346 + /* y := P'*w */
1.347 + for (i = 1; i <= m; i++) y[i] = w[P[m+i]];
1.348 + xfree(w);
1.349 + /* compute residual vector r = A*D*A'*y - h */
1.350 + r = xcalloc(1+m, sizeof(double));
1.351 + /* w := A'*y */
1.352 + w = xcalloc(1+n, sizeof(double));
1.353 + AT_by_vec(csa, y, w);
1.354 + /* w := D*w */
1.355 + for (j = 1; j <= n; j++) w[j] *= csa->D[j];
1.356 + /* r := A*w */
1.357 + A_by_vec(csa, w, r);
1.358 + xfree(w);
1.359 + /* r := r - h */
1.360 + for (i = 1; i <= m; i++) r[i] -= h[i];
1.361 + /* check for numeric stability */
1.362 + for (i = 1; i <= m; i++)
1.363 + { if (fabs(r[i]) / (1.0 + fabs(h[i])) > 1e-4)
1.364 + { ret = 1;
1.365 + break;
1.366 + }
1.367 + }
1.368 + xfree(h);
1.369 + xfree(r);
1.370 + return ret;
1.371 +}
1.372 +
1.373 +/***********************************************************************
1.374 +* solve_NS - solve Newtonian system
1.375 +*
1.376 +* This routine solves the Newtonian system:
1.377 +*
1.378 +* A*dx = p
1.379 +*
1.380 +* A'*dy + dz = q
1.381 +*
1.382 +* Z*dx + X*dz = r
1.383 +*
1.384 +* where X = diag(x[j]), Z = diag(z[j]), by reducing it to the normal
1.385 +* equation system:
1.386 +*
1.387 +* (A*inv(Z)*X*A')*dy = A*inv(Z)*(X*q-r)+p
1.388 +*
1.389 +* (it is assumed that the matrix A*inv(Z)*X*A' has been factorized by
1.390 +* the routine decomp_NE).
1.391 +*
1.392 +* Once vector dy has been computed the routine computes vectors dx and
1.393 +* dz as follows:
1.394 +*
1.395 +* dx = inv(Z)*(X*(A'*dy-q)+r)
1.396 +*
1.397 +* dz = inv(X)*(r-Z*dx)
1.398 +*
1.399 +* The routine solve_NS returns the same code which was reported by the
1.400 +* routine solve_NE (see above). */
1.401 +
1.402 +static int solve_NS(struct csa *csa, double p[], double q[], double r[],
1.403 + double dx[], double dy[], double dz[])
1.404 +{ int m = csa->m;
1.405 + int n = csa->n;
1.406 + double *x = csa->x;
1.407 + double *z = csa->z;
1.408 + int i, j, ret;
1.409 + double *w = dx;
1.410 + /* compute the vector of right-hand sides A*inv(Z)*(X*q-r)+p for
1.411 + the normal equation system */
1.412 + for (j = 1; j <= n; j++)
1.413 + w[j] = (x[j] * q[j] - r[j]) / z[j];
1.414 + A_by_vec(csa, w, dy);
1.415 + for (i = 1; i <= m; i++) dy[i] += p[i];
1.416 + /* solve the normal equation system to compute vector dy */
1.417 + ret = solve_NE(csa, dy);
1.418 + /* compute vectors dx and dz */
1.419 + AT_by_vec(csa, dy, dx);
1.420 + for (j = 1; j <= n; j++)
1.421 + { dx[j] = (x[j] * (dx[j] - q[j]) + r[j]) / z[j];
1.422 + dz[j] = (r[j] - z[j] * dx[j]) / x[j];
1.423 + }
1.424 + return ret;
1.425 +}
1.426 +
1.427 +/***********************************************************************
1.428 +* initial_point - choose initial point using Mehrotra's heuristic
1.429 +*
1.430 +* This routine chooses a starting point using a heuristic proposed in
1.431 +* the paper:
1.432 +*
1.433 +* S. Mehrotra. On the implementation of a primal-dual interior point
1.434 +* method. SIAM J. on Optim., 2(4), pp. 575-601, 1992.
1.435 +*
1.436 +* The starting point x in the primal space is chosen as a solution of
1.437 +* the following least squares problem:
1.438 +*
1.439 +* minimize ||x||
1.440 +*
1.441 +* subject to A*x = b
1.442 +*
1.443 +* which can be computed explicitly as follows:
1.444 +*
1.445 +* x = A'*inv(A*A')*b
1.446 +*
1.447 +* Similarly, the starting point (y, z) in the dual space is chosen as
1.448 +* a solution of the following least squares problem:
1.449 +*
1.450 +* minimize ||z||
1.451 +*
1.452 +* subject to A'*y + z = c
1.453 +*
1.454 +* which can be computed explicitly as follows:
1.455 +*
1.456 +* y = inv(A*A')*A*c
1.457 +*
1.458 +* z = c - A'*y
1.459 +*
1.460 +* However, some components of the vectors x and z may be non-positive
1.461 +* or close to zero, so the routine uses a Mehrotra's heuristic to find
1.462 +* a more appropriate starting point. */
1.463 +
1.464 +static void initial_point(struct csa *csa)
1.465 +{ int m = csa->m;
1.466 + int n = csa->n;
1.467 + double *b = csa->b;
1.468 + double *c = csa->c;
1.469 + double *x = csa->x;
1.470 + double *y = csa->y;
1.471 + double *z = csa->z;
1.472 + double *D = csa->D;
1.473 + int i, j;
1.474 + double dp, dd, ex, ez, xz;
1.475 + /* factorize A*A' */
1.476 + for (j = 1; j <= n; j++) D[j] = 1.0;
1.477 + decomp_NE(csa);
1.478 + /* x~ = A'*inv(A*A')*b */
1.479 + for (i = 1; i <= m; i++) y[i] = b[i];
1.480 + solve_NE(csa, y);
1.481 + AT_by_vec(csa, y, x);
1.482 + /* y~ = inv(A*A')*A*c */
1.483 + A_by_vec(csa, c, y);
1.484 + solve_NE(csa, y);
1.485 + /* z~ = c - A'*y~ */
1.486 + AT_by_vec(csa, y,z);
1.487 + for (j = 1; j <= n; j++) z[j] = c[j] - z[j];
1.488 + /* use Mehrotra's heuristic in order to choose more appropriate
1.489 + starting point with positive components of vectors x and z */
1.490 + dp = dd = 0.0;
1.491 + for (j = 1; j <= n; j++)
1.492 + { if (dp < -1.5 * x[j]) dp = -1.5 * x[j];
1.493 + if (dd < -1.5 * z[j]) dd = -1.5 * z[j];
1.494 + }
1.495 + /* note that b = 0 involves x = 0, and c = 0 involves y = 0 and
1.496 + z = 0, so we need to be careful */
1.497 + if (dp == 0.0) dp = 1.5;
1.498 + if (dd == 0.0) dd = 1.5;
1.499 + ex = ez = xz = 0.0;
1.500 + for (j = 1; j <= n; j++)
1.501 + { ex += (x[j] + dp);
1.502 + ez += (z[j] + dd);
1.503 + xz += (x[j] + dp) * (z[j] + dd);
1.504 + }
1.505 + dp += 0.5 * (xz / ez);
1.506 + dd += 0.5 * (xz / ex);
1.507 + for (j = 1; j <= n; j++)
1.508 + { x[j] += dp;
1.509 + z[j] += dd;
1.510 + xassert(x[j] > 0.0 && z[j] > 0.0);
1.511 + }
1.512 + return;
1.513 +}
1.514 +
1.515 +/***********************************************************************
1.516 +* basic_info - perform basic computations at the current point
1.517 +*
1.518 +* This routine computes the following quantities at the current point:
1.519 +*
1.520 +* 1) value of the objective function:
1.521 +*
1.522 +* F = c'*x + c[0]
1.523 +*
1.524 +* 2) relative primal infeasibility:
1.525 +*
1.526 +* rpi = ||A*x-b|| / (1+||b||)
1.527 +*
1.528 +* 3) relative dual infeasibility:
1.529 +*
1.530 +* rdi = ||A'*y+z-c|| / (1+||c||)
1.531 +*
1.532 +* 4) primal-dual gap (relative difference between the primal and the
1.533 +* dual objective function values):
1.534 +*
1.535 +* gap = |c'*x-b'*y| / (1+|c'*x|)
1.536 +*
1.537 +* 5) merit function:
1.538 +*
1.539 +* phi = ||A*x-b|| / max(1,||b||) + ||A'*y+z-c|| / max(1,||c||) +
1.540 +*
1.541 +* + |c'*x-b'*y| / max(1,||b||,||c||)
1.542 +*
1.543 +* 6) duality measure:
1.544 +*
1.545 +* mu = x'*z / n
1.546 +*
1.547 +* 7) the ratio of infeasibility to mu:
1.548 +*
1.549 +* rmu = max(||A*x-b||,||A'*y+z-c||) / mu
1.550 +*
1.551 +* where ||*|| denotes euclidian norm, *' denotes transposition. */
1.552 +
1.553 +static void basic_info(struct csa *csa)
1.554 +{ int m = csa->m;
1.555 + int n = csa->n;
1.556 + double *b = csa->b;
1.557 + double *c = csa->c;
1.558 + double *x = csa->x;
1.559 + double *y = csa->y;
1.560 + double *z = csa->z;
1.561 + int i, j;
1.562 + double norm1, bnorm, norm2, cnorm, cx, by, *work, temp;
1.563 + /* compute value of the objective function */
1.564 + temp = c[0];
1.565 + for (j = 1; j <= n; j++) temp += c[j] * x[j];
1.566 + csa->obj = temp;
1.567 + /* norm1 = ||A*x-b|| */
1.568 + work = xcalloc(1+m, sizeof(double));
1.569 + A_by_vec(csa, x, work);
1.570 + norm1 = 0.0;
1.571 + for (i = 1; i <= m; i++)
1.572 + norm1 += (work[i] - b[i]) * (work[i] - b[i]);
1.573 + norm1 = sqrt(norm1);
1.574 + xfree(work);
1.575 + /* bnorm = ||b|| */
1.576 + bnorm = 0.0;
1.577 + for (i = 1; i <= m; i++) bnorm += b[i] * b[i];
1.578 + bnorm = sqrt(bnorm);
1.579 + /* compute relative primal infeasibility */
1.580 + csa->rpi = norm1 / (1.0 + bnorm);
1.581 + /* norm2 = ||A'*y+z-c|| */
1.582 + work = xcalloc(1+n, sizeof(double));
1.583 + AT_by_vec(csa, y, work);
1.584 + norm2 = 0.0;
1.585 + for (j = 1; j <= n; j++)
1.586 + norm2 += (work[j] + z[j] - c[j]) * (work[j] + z[j] - c[j]);
1.587 + norm2 = sqrt(norm2);
1.588 + xfree(work);
1.589 + /* cnorm = ||c|| */
1.590 + cnorm = 0.0;
1.591 + for (j = 1; j <= n; j++) cnorm += c[j] * c[j];
1.592 + cnorm = sqrt(cnorm);
1.593 + /* compute relative dual infeasibility */
1.594 + csa->rdi = norm2 / (1.0 + cnorm);
1.595 + /* by = b'*y */
1.596 + by = 0.0;
1.597 + for (i = 1; i <= m; i++) by += b[i] * y[i];
1.598 + /* cx = c'*x */
1.599 + cx = 0.0;
1.600 + for (j = 1; j <= n; j++) cx += c[j] * x[j];
1.601 + /* compute primal-dual gap */
1.602 + csa->gap = fabs(cx - by) / (1.0 + fabs(cx));
1.603 + /* compute merit function */
1.604 + csa->phi = 0.0;
1.605 + csa->phi += norm1 / (bnorm > 1.0 ? bnorm : 1.0);
1.606 + csa->phi += norm2 / (cnorm > 1.0 ? cnorm : 1.0);
1.607 + temp = 1.0;
1.608 + if (temp < bnorm) temp = bnorm;
1.609 + if (temp < cnorm) temp = cnorm;
1.610 + csa->phi += fabs(cx - by) / temp;
1.611 + /* compute duality measure */
1.612 + temp = 0.0;
1.613 + for (j = 1; j <= n; j++) temp += x[j] * z[j];
1.614 + csa->mu = temp / (double)n;
1.615 + /* compute the ratio of infeasibility to mu */
1.616 + csa->rmu = (norm1 > norm2 ? norm1 : norm2) / csa->mu;
1.617 + return;
1.618 +}
1.619 +
1.620 +/***********************************************************************
1.621 +* make_step - compute next point using Mehrotra's technique
1.622 +*
1.623 +* This routine computes the next point using the predictor-corrector
1.624 +* technique proposed in the paper:
1.625 +*
1.626 +* S. Mehrotra. On the implementation of a primal-dual interior point
1.627 +* method. SIAM J. on Optim., 2(4), pp. 575-601, 1992.
1.628 +*
1.629 +* At first, the routine computes so called affine scaling (predictor)
1.630 +* direction (dx_aff,dy_aff,dz_aff) which is a solution of the system:
1.631 +*
1.632 +* A*dx_aff = b - A*x
1.633 +*
1.634 +* A'*dy_aff + dz_aff = c - A'*y - z
1.635 +*
1.636 +* Z*dx_aff + X*dz_aff = - X*Z*e
1.637 +*
1.638 +* where (x,y,z) is the current point, X = diag(x[j]), Z = diag(z[j]),
1.639 +* e = (1,...,1)'.
1.640 +*
1.641 +* Then, the routine computes the centering parameter sigma, using the
1.642 +* following Mehrotra's heuristic:
1.643 +*
1.644 +* alfa_aff_p = inf{0 <= alfa <= 1 | x+alfa*dx_aff >= 0}
1.645 +*
1.646 +* alfa_aff_d = inf{0 <= alfa <= 1 | z+alfa*dz_aff >= 0}
1.647 +*
1.648 +* mu_aff = (x+alfa_aff_p*dx_aff)'*(z+alfa_aff_d*dz_aff)/n
1.649 +*
1.650 +* sigma = (mu_aff/mu)^3
1.651 +*
1.652 +* where alfa_aff_p is the maximal stepsize along the affine scaling
1.653 +* direction in the primal space, alfa_aff_d is the maximal stepsize
1.654 +* along the same direction in the dual space.
1.655 +*
1.656 +* After determining sigma the routine computes so called centering
1.657 +* (corrector) direction (dx_cc,dy_cc,dz_cc) which is the solution of
1.658 +* the system:
1.659 +*
1.660 +* A*dx_cc = 0
1.661 +*
1.662 +* A'*dy_cc + dz_cc = 0
1.663 +*
1.664 +* Z*dx_cc + X*dz_cc = sigma*mu*e - X*Z*e
1.665 +*
1.666 +* Finally, the routine computes the combined direction
1.667 +*
1.668 +* (dx,dy,dz) = (dx_aff,dy_aff,dz_aff) + (dx_cc,dy_cc,dz_cc)
1.669 +*
1.670 +* and determines maximal primal and dual stepsizes along the combined
1.671 +* direction:
1.672 +*
1.673 +* alfa_max_p = inf{0 <= alfa <= 1 | x+alfa*dx >= 0}
1.674 +*
1.675 +* alfa_max_d = inf{0 <= alfa <= 1 | z+alfa*dz >= 0}
1.676 +*
1.677 +* In order to prevent the next point to be too close to the boundary
1.678 +* of the positive ortant, the routine decreases maximal stepsizes:
1.679 +*
1.680 +* alfa_p = gamma_p * alfa_max_p
1.681 +*
1.682 +* alfa_d = gamma_d * alfa_max_d
1.683 +*
1.684 +* where gamma_p and gamma_d are scaling factors, and computes the next
1.685 +* point:
1.686 +*
1.687 +* x_new = x + alfa_p * dx
1.688 +*
1.689 +* y_new = y + alfa_d * dy
1.690 +*
1.691 +* z_new = z + alfa_d * dz
1.692 +*
1.693 +* which becomes the current point on the next iteration. */
1.694 +
1.695 +static int make_step(struct csa *csa)
1.696 +{ int m = csa->m;
1.697 + int n = csa->n;
1.698 + double *b = csa->b;
1.699 + double *c = csa->c;
1.700 + double *x = csa->x;
1.701 + double *y = csa->y;
1.702 + double *z = csa->z;
1.703 + double *dx_aff = csa->dx_aff;
1.704 + double *dy_aff = csa->dy_aff;
1.705 + double *dz_aff = csa->dz_aff;
1.706 + double *dx_cc = csa->dx_cc;
1.707 + double *dy_cc = csa->dy_cc;
1.708 + double *dz_cc = csa->dz_cc;
1.709 + double *dx = csa->dx;
1.710 + double *dy = csa->dy;
1.711 + double *dz = csa->dz;
1.712 + int i, j, ret = 0;
1.713 + double temp, gamma_p, gamma_d, *p, *q, *r;
1.714 + /* allocate working arrays */
1.715 + p = xcalloc(1+m, sizeof(double));
1.716 + q = xcalloc(1+n, sizeof(double));
1.717 + r = xcalloc(1+n, sizeof(double));
1.718 + /* p = b - A*x */
1.719 + A_by_vec(csa, x, p);
1.720 + for (i = 1; i <= m; i++) p[i] = b[i] - p[i];
1.721 + /* q = c - A'*y - z */
1.722 + AT_by_vec(csa, y,q);
1.723 + for (j = 1; j <= n; j++) q[j] = c[j] - q[j] - z[j];
1.724 + /* r = - X * Z * e */
1.725 + for (j = 1; j <= n; j++) r[j] = - x[j] * z[j];
1.726 + /* solve the first Newtonian system */
1.727 + if (solve_NS(csa, p, q, r, dx_aff, dy_aff, dz_aff))
1.728 + { ret = 1;
1.729 + goto done;
1.730 + }
1.731 + /* alfa_aff_p = inf{0 <= alfa <= 1 | x + alfa*dx_aff >= 0} */
1.732 + /* alfa_aff_d = inf{0 <= alfa <= 1 | z + alfa*dz_aff >= 0} */
1.733 + csa->alfa_aff_p = csa->alfa_aff_d = 1.0;
1.734 + for (j = 1; j <= n; j++)
1.735 + { if (dx_aff[j] < 0.0)
1.736 + { temp = - x[j] / dx_aff[j];
1.737 + if (csa->alfa_aff_p > temp) csa->alfa_aff_p = temp;
1.738 + }
1.739 + if (dz_aff[j] < 0.0)
1.740 + { temp = - z[j] / dz_aff[j];
1.741 + if (csa->alfa_aff_d > temp) csa->alfa_aff_d = temp;
1.742 + }
1.743 + }
1.744 + /* mu_aff = (x+alfa_aff_p*dx_aff)' * (z+alfa_aff_d*dz_aff) / n */
1.745 + temp = 0.0;
1.746 + for (j = 1; j <= n; j++)
1.747 + temp += (x[j] + csa->alfa_aff_p * dx_aff[j]) *
1.748 + (z[j] + csa->alfa_aff_d * dz_aff[j]);
1.749 + csa->mu_aff = temp / (double)n;
1.750 + /* sigma = (mu_aff/mu)^3 */
1.751 + temp = csa->mu_aff / csa->mu;
1.752 + csa->sigma = temp * temp * temp;
1.753 + /* p = 0 */
1.754 + for (i = 1; i <= m; i++) p[i] = 0.0;
1.755 + /* q = 0 */
1.756 + for (j = 1; j <= n; j++) q[j] = 0.0;
1.757 + /* r = sigma * mu * e - X * Z * e */
1.758 + for (j = 1; j <= n; j++)
1.759 + r[j] = csa->sigma * csa->mu - dx_aff[j] * dz_aff[j];
1.760 + /* solve the second Newtonian system with the same coefficients
1.761 + but with altered right-hand sides */
1.762 + if (solve_NS(csa, p, q, r, dx_cc, dy_cc, dz_cc))
1.763 + { ret = 1;
1.764 + goto done;
1.765 + }
1.766 + /* (dx,dy,dz) = (dx_aff,dy_aff,dz_aff) + (dx_cc,dy_cc,dz_cc) */
1.767 + for (j = 1; j <= n; j++) dx[j] = dx_aff[j] + dx_cc[j];
1.768 + for (i = 1; i <= m; i++) dy[i] = dy_aff[i] + dy_cc[i];
1.769 + for (j = 1; j <= n; j++) dz[j] = dz_aff[j] + dz_cc[j];
1.770 + /* alfa_max_p = inf{0 <= alfa <= 1 | x + alfa*dx >= 0} */
1.771 + /* alfa_max_d = inf{0 <= alfa <= 1 | z + alfa*dz >= 0} */
1.772 + csa->alfa_max_p = csa->alfa_max_d = 1.0;
1.773 + for (j = 1; j <= n; j++)
1.774 + { if (dx[j] < 0.0)
1.775 + { temp = - x[j] / dx[j];
1.776 + if (csa->alfa_max_p > temp) csa->alfa_max_p = temp;
1.777 + }
1.778 + if (dz[j] < 0.0)
1.779 + { temp = - z[j] / dz[j];
1.780 + if (csa->alfa_max_d > temp) csa->alfa_max_d = temp;
1.781 + }
1.782 + }
1.783 + /* determine scale factors (not implemented yet) */
1.784 + gamma_p = 0.90;
1.785 + gamma_d = 0.90;
1.786 + /* compute the next point */
1.787 + for (j = 1; j <= n; j++)
1.788 + { x[j] += gamma_p * csa->alfa_max_p * dx[j];
1.789 + xassert(x[j] > 0.0);
1.790 + }
1.791 + for (i = 1; i <= m; i++)
1.792 + y[i] += gamma_d * csa->alfa_max_d * dy[i];
1.793 + for (j = 1; j <= n; j++)
1.794 + { z[j] += gamma_d * csa->alfa_max_d * dz[j];
1.795 + xassert(z[j] > 0.0);
1.796 + }
1.797 +done: /* free working arrays */
1.798 + xfree(p);
1.799 + xfree(q);
1.800 + xfree(r);
1.801 + return ret;
1.802 +}
1.803 +
1.804 +/***********************************************************************
1.805 +* terminate - deallocate common storage area
1.806 +*
1.807 +* This routine frees all memory allocated to the common storage area
1.808 +* used by interior-point method routines. */
1.809 +
1.810 +static void terminate(struct csa *csa)
1.811 +{ xfree(csa->D);
1.812 + xfree(csa->P);
1.813 + xfree(csa->S_ptr);
1.814 + xfree(csa->S_ind);
1.815 + xfree(csa->S_val);
1.816 + xfree(csa->S_diag);
1.817 + xfree(csa->U_ptr);
1.818 + xfree(csa->U_ind);
1.819 + xfree(csa->U_val);
1.820 + xfree(csa->U_diag);
1.821 + xfree(csa->phi_min);
1.822 + xfree(csa->best_x);
1.823 + xfree(csa->best_y);
1.824 + xfree(csa->best_z);
1.825 + xfree(csa->dx_aff);
1.826 + xfree(csa->dy_aff);
1.827 + xfree(csa->dz_aff);
1.828 + xfree(csa->dx_cc);
1.829 + xfree(csa->dy_cc);
1.830 + xfree(csa->dz_cc);
1.831 + return;
1.832 +}
1.833 +
1.834 +/***********************************************************************
1.835 +* ipm_main - main interior-point method routine
1.836 +*
1.837 +* This is a main routine of the primal-dual interior-point method.
1.838 +*
1.839 +* The routine ipm_main returns one of the following codes:
1.840 +*
1.841 +* 0 - optimal solution found;
1.842 +* 1 - problem has no feasible (primal or dual) solution;
1.843 +* 2 - no convergence;
1.844 +* 3 - iteration limit exceeded;
1.845 +* 4 - numeric instability on solving Newtonian system.
1.846 +*
1.847 +* In case of non-zero return code the routine returns the best point,
1.848 +* which has been reached during optimization. */
1.849 +
1.850 +static int ipm_main(struct csa *csa)
1.851 +{ int m = csa->m;
1.852 + int n = csa->n;
1.853 + int i, j, status;
1.854 + double temp;
1.855 + /* choose initial point using Mehrotra's heuristic */
1.856 + if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.857 + xprintf("Guessing initial point...\n");
1.858 + initial_point(csa);
1.859 + /* main loop starts here */
1.860 + if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.861 + xprintf("Optimization begins...\n");
1.862 + for (;;)
1.863 + { /* perform basic computations at the current point */
1.864 + basic_info(csa);
1.865 + /* save initial value of rmu */
1.866 + if (csa->iter == 0) csa->rmu0 = csa->rmu;
1.867 + /* accumulate values of min(phi[k]) and save the best point */
1.868 + xassert(csa->iter <= ITER_MAX);
1.869 + if (csa->iter == 0 || csa->phi_min[csa->iter-1] > csa->phi)
1.870 + { csa->phi_min[csa->iter] = csa->phi;
1.871 + csa->best_iter = csa->iter;
1.872 + for (j = 1; j <= n; j++) csa->best_x[j] = csa->x[j];
1.873 + for (i = 1; i <= m; i++) csa->best_y[i] = csa->y[i];
1.874 + for (j = 1; j <= n; j++) csa->best_z[j] = csa->z[j];
1.875 + csa->best_obj = csa->obj;
1.876 + }
1.877 + else
1.878 + csa->phi_min[csa->iter] = csa->phi_min[csa->iter-1];
1.879 + /* display information at the current point */
1.880 + if (csa->parm->msg_lev >= GLP_MSG_ON)
1.881 + xprintf("%3d: obj = %17.9e; rpi = %8.1e; rdi = %8.1e; gap ="
1.882 + " %8.1e\n", csa->iter, csa->obj, csa->rpi, csa->rdi,
1.883 + csa->gap);
1.884 + /* check if the current point is optimal */
1.885 + if (csa->rpi < 1e-8 && csa->rdi < 1e-8 && csa->gap < 1e-8)
1.886 + { if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.887 + xprintf("OPTIMAL SOLUTION FOUND\n");
1.888 + status = 0;
1.889 + break;
1.890 + }
1.891 + /* check if the problem has no feasible solution */
1.892 + temp = 1e5 * csa->phi_min[csa->iter];
1.893 + if (temp < 1e-8) temp = 1e-8;
1.894 + if (csa->phi >= temp)
1.895 + { if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.896 + xprintf("PROBLEM HAS NO FEASIBLE PRIMAL/DUAL SOLUTION\n")
1.897 + ;
1.898 + status = 1;
1.899 + break;
1.900 + }
1.901 + /* check for very slow convergence or divergence */
1.902 + if (((csa->rpi >= 1e-8 || csa->rdi >= 1e-8) && csa->rmu /
1.903 + csa->rmu0 >= 1e6) ||
1.904 + (csa->iter >= 30 && csa->phi_min[csa->iter] >= 0.5 *
1.905 + csa->phi_min[csa->iter - 30]))
1.906 + { if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.907 + xprintf("NO CONVERGENCE; SEARCH TERMINATED\n");
1.908 + status = 2;
1.909 + break;
1.910 + }
1.911 + /* check for maximal number of iterations */
1.912 + if (csa->iter == ITER_MAX)
1.913 + { if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.914 + xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n");
1.915 + status = 3;
1.916 + break;
1.917 + }
1.918 + /* start the next iteration */
1.919 + csa->iter++;
1.920 + /* factorize normal equation system */
1.921 + for (j = 1; j <= n; j++) csa->D[j] = csa->x[j] / csa->z[j];
1.922 + decomp_NE(csa);
1.923 + /* compute the next point using Mehrotra's predictor-corrector
1.924 + technique */
1.925 + if (make_step(csa))
1.926 + { if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.927 + xprintf("NUMERIC INSTABILITY; SEARCH TERMINATED\n");
1.928 + status = 4;
1.929 + break;
1.930 + }
1.931 + }
1.932 + /* restore the best point */
1.933 + if (status != 0)
1.934 + { for (j = 1; j <= n; j++) csa->x[j] = csa->best_x[j];
1.935 + for (i = 1; i <= m; i++) csa->y[i] = csa->best_y[i];
1.936 + for (j = 1; j <= n; j++) csa->z[j] = csa->best_z[j];
1.937 + if (csa->parm->msg_lev >= GLP_MSG_ALL)
1.938 + xprintf("Best point %17.9e was reached on iteration %d\n",
1.939 + csa->best_obj, csa->best_iter);
1.940 + }
1.941 + /* return to the calling program */
1.942 + return status;
1.943 +}
1.944 +
1.945 +/***********************************************************************
1.946 +* NAME
1.947 +*
1.948 +* ipm_solve - core LP solver based on the interior-point method
1.949 +*
1.950 +* SYNOPSIS
1.951 +*
1.952 +* #include "glpipm.h"
1.953 +* int ipm_solve(glp_prob *P, const glp_iptcp *parm);
1.954 +*
1.955 +* DESCRIPTION
1.956 +*
1.957 +* The routine ipm_solve is a core LP solver based on the primal-dual
1.958 +* interior-point method.
1.959 +*
1.960 +* The routine assumes the following standard formulation of LP problem
1.961 +* to be solved:
1.962 +*
1.963 +* minimize
1.964 +*
1.965 +* F = c[0] + c[1]*x[1] + c[2]*x[2] + ... + c[n]*x[n]
1.966 +*
1.967 +* subject to linear constraints
1.968 +*
1.969 +* a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1]
1.970 +*
1.971 +* a[2,1]*x[1] + a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2]
1.972 +*
1.973 +* . . . . . .
1.974 +*
1.975 +* a[m,1]*x[1] + a[m,2]*x[2] + ... + a[m,n]*x[n] = b[m]
1.976 +*
1.977 +* and non-negative variables
1.978 +*
1.979 +* x[1] >= 0, x[2] >= 0, ..., x[n] >= 0
1.980 +*
1.981 +* where:
1.982 +* F is the objective function;
1.983 +* x[1], ..., x[n] are (structural) variables;
1.984 +* c[0] is a constant term of the objective function;
1.985 +* c[1], ..., c[n] are objective coefficients;
1.986 +* a[1,1], ..., a[m,n] are constraint coefficients;
1.987 +* b[1], ..., b[n] are right-hand sides.
1.988 +*
1.989 +* The solution is three vectors x, y, and z, which are stored by the
1.990 +* routine in the arrays x, y, and z, respectively. These vectors
1.991 +* correspond to the best primal-dual point found during optimization.
1.992 +* They are approximate solution of the following system (which is the
1.993 +* Karush-Kuhn-Tucker optimality conditions):
1.994 +*
1.995 +* A*x = b (primal feasibility condition)
1.996 +*
1.997 +* A'*y + z = c (dual feasibility condition)
1.998 +*
1.999 +* x'*z = 0 (primal-dual complementarity condition)
1.1000 +*
1.1001 +* x >= 0, z >= 0 (non-negativity condition)
1.1002 +*
1.1003 +* where:
1.1004 +* x[1], ..., x[n] are primal (structural) variables;
1.1005 +* y[1], ..., y[m] are dual variables (Lagrange multipliers) for
1.1006 +* equality constraints;
1.1007 +* z[1], ..., z[n] are dual variables (Lagrange multipliers) for
1.1008 +* non-negativity constraints.
1.1009 +*
1.1010 +* RETURNS
1.1011 +*
1.1012 +* 0 LP has been successfully solved.
1.1013 +*
1.1014 +* GLP_ENOCVG
1.1015 +* No convergence.
1.1016 +*
1.1017 +* GLP_EITLIM
1.1018 +* Iteration limit exceeded.
1.1019 +*
1.1020 +* GLP_EINSTAB
1.1021 +* Numeric instability on solving Newtonian system.
1.1022 +*
1.1023 +* In case of non-zero return code the routine returns the best point,
1.1024 +* which has been reached during optimization. */
1.1025 +
1.1026 +int ipm_solve(glp_prob *P, const glp_iptcp *parm)
1.1027 +{ struct csa _dsa, *csa = &_dsa;
1.1028 + int m = P->m;
1.1029 + int n = P->n;
1.1030 + int nnz = P->nnz;
1.1031 + GLPROW *row;
1.1032 + GLPCOL *col;
1.1033 + GLPAIJ *aij;
1.1034 + int i, j, loc, ret, *A_ind, *A_ptr;
1.1035 + double dir, *A_val, *b, *c, *x, *y, *z;
1.1036 + xassert(m > 0);
1.1037 + xassert(n > 0);
1.1038 + /* allocate working arrays */
1.1039 + A_ptr = xcalloc(1+m+1, sizeof(int));
1.1040 + A_ind = xcalloc(1+nnz, sizeof(int));
1.1041 + A_val = xcalloc(1+nnz, sizeof(double));
1.1042 + b = xcalloc(1+m, sizeof(double));
1.1043 + c = xcalloc(1+n, sizeof(double));
1.1044 + x = xcalloc(1+n, sizeof(double));
1.1045 + y = xcalloc(1+m, sizeof(double));
1.1046 + z = xcalloc(1+n, sizeof(double));
1.1047 + /* prepare rows and constraint coefficients */
1.1048 + loc = 1;
1.1049 + for (i = 1; i <= m; i++)
1.1050 + { row = P->row[i];
1.1051 + xassert(row->type == GLP_FX);
1.1052 + b[i] = row->lb * row->rii;
1.1053 + A_ptr[i] = loc;
1.1054 + for (aij = row->ptr; aij != NULL; aij = aij->r_next)
1.1055 + { A_ind[loc] = aij->col->j;
1.1056 + A_val[loc] = row->rii * aij->val * aij->col->sjj;
1.1057 + loc++;
1.1058 + }
1.1059 + }
1.1060 + A_ptr[m+1] = loc;
1.1061 + xassert(loc-1 == nnz);
1.1062 + /* prepare columns and objective coefficients */
1.1063 + if (P->dir == GLP_MIN)
1.1064 + dir = +1.0;
1.1065 + else if (P->dir == GLP_MAX)
1.1066 + dir = -1.0;
1.1067 + else
1.1068 + xassert(P != P);
1.1069 + c[0] = dir * P->c0;
1.1070 + for (j = 1; j <= n; j++)
1.1071 + { col = P->col[j];
1.1072 + xassert(col->type == GLP_LO && col->lb == 0.0);
1.1073 + c[j] = dir * col->coef * col->sjj;
1.1074 + }
1.1075 + /* allocate and initialize the common storage area */
1.1076 + csa->m = m;
1.1077 + csa->n = n;
1.1078 + csa->A_ptr = A_ptr;
1.1079 + csa->A_ind = A_ind;
1.1080 + csa->A_val = A_val;
1.1081 + csa->b = b;
1.1082 + csa->c = c;
1.1083 + csa->x = x;
1.1084 + csa->y = y;
1.1085 + csa->z = z;
1.1086 + csa->parm = parm;
1.1087 + initialize(csa);
1.1088 + /* solve LP with the interior-point method */
1.1089 + ret = ipm_main(csa);
1.1090 + /* deallocate the common storage area */
1.1091 + terminate(csa);
1.1092 + /* determine solution status */
1.1093 + if (ret == 0)
1.1094 + { /* optimal solution found */
1.1095 + P->ipt_stat = GLP_OPT;
1.1096 + ret = 0;
1.1097 + }
1.1098 + else if (ret == 1)
1.1099 + { /* problem has no feasible (primal or dual) solution */
1.1100 + P->ipt_stat = GLP_NOFEAS;
1.1101 + ret = 0;
1.1102 + }
1.1103 + else if (ret == 2)
1.1104 + { /* no convergence */
1.1105 + P->ipt_stat = GLP_INFEAS;
1.1106 + ret = GLP_ENOCVG;
1.1107 + }
1.1108 + else if (ret == 3)
1.1109 + { /* iteration limit exceeded */
1.1110 + P->ipt_stat = GLP_INFEAS;
1.1111 + ret = GLP_EITLIM;
1.1112 + }
1.1113 + else if (ret == 4)
1.1114 + { /* numeric instability on solving Newtonian system */
1.1115 + P->ipt_stat = GLP_INFEAS;
1.1116 + ret = GLP_EINSTAB;
1.1117 + }
1.1118 + else
1.1119 + xassert(ret != ret);
1.1120 + /* store row solution components */
1.1121 + for (i = 1; i <= m; i++)
1.1122 + { row = P->row[i];
1.1123 + row->pval = row->lb;
1.1124 + row->dval = dir * y[i] * row->rii;
1.1125 + }
1.1126 + /* store column solution components */
1.1127 + P->ipt_obj = P->c0;
1.1128 + for (j = 1; j <= n; j++)
1.1129 + { col = P->col[j];
1.1130 + col->pval = x[j] * col->sjj;
1.1131 + col->dval = dir * z[j] / col->sjj;
1.1132 + P->ipt_obj += col->coef * col->pval;
1.1133 + }
1.1134 + /* free working arrays */
1.1135 + xfree(A_ptr);
1.1136 + xfree(A_ind);
1.1137 + xfree(A_val);
1.1138 + xfree(b);
1.1139 + xfree(c);
1.1140 + xfree(x);
1.1141 + xfree(y);
1.1142 + xfree(z);
1.1143 + return ret;
1.1144 +}
1.1145 +
1.1146 +/* eof */