[1] | 1 | /* glpipm.c */ |
---|
| 2 | |
---|
| 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 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 | ***********************************************************************/ |
---|
| 24 | |
---|
| 25 | #include "glpipm.h" |
---|
| 26 | #include "glpmat.h" |
---|
| 27 | |
---|
| 28 | #define ITER_MAX 100 |
---|
| 29 | /* maximal number of iterations */ |
---|
| 30 | |
---|
| 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 | }; |
---|
| 142 | |
---|
| 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. */ |
---|
| 148 | |
---|
| 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 | } |
---|
| 238 | |
---|
| 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. */ |
---|
| 244 | |
---|
| 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 | } |
---|
| 261 | |
---|
| 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. */ |
---|
| 267 | |
---|
| 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 | } |
---|
| 286 | |
---|
| 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. */ |
---|
| 293 | |
---|
| 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 | } |
---|
| 302 | |
---|
| 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. */ |
---|
| 323 | |
---|
| 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 | } |
---|
| 369 | |
---|
| 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). */ |
---|
| 398 | |
---|
| 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 | } |
---|
| 423 | |
---|
| 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. */ |
---|
| 460 | |
---|
| 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 | } |
---|
| 511 | |
---|
| 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. */ |
---|
| 549 | |
---|
| 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 | } |
---|
| 616 | |
---|
| 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. */ |
---|
| 691 | |
---|
| 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 | } |
---|
| 800 | |
---|
| 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. */ |
---|
| 806 | |
---|
| 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 | } |
---|
| 830 | |
---|
| 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. */ |
---|
| 846 | |
---|
| 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 | } |
---|
| 941 | |
---|
| 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. */ |
---|
| 1022 | |
---|
| 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 | } |
---|
| 1142 | |
---|
| 1143 | /* eof */ |
---|