/farmR/src/glpipm.c

https://code.google.com/p/javawfm/ · C · 983 lines · 497 code · 28 blank · 458 comment · 93 complexity · 7fb7967375c66c2e86f1f7bfd3e1779e MD5 · raw file

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