PageRenderTime 95ms CodeModel.GetById 11ms app.highlight 76ms RepoModel.GetById 2ms app.codeStats 0ms

/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/***********************************************************************
  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 */