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