COIN-OR::LEMON - Graph Library

source: glpk-cmake/src/glpipm.c @ 1:c445c931472f

Last change on this file since 1:c445c931472f was 1:c445c931472f, checked in by Alpar Juttner <alpar@…>, 10 years ago

Import glpk-4.45

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