1.1 --- /dev/null	Thu Jan 01 00:00:00 1970 +0000
     1.2 +++ b/deps/glpk/src/glpipm.c	Sun Nov 06 20:59:10 2011 +0100
     1.3 @@ -0,0 +1,1143 @@
     1.4 +/* glpipm.c */
     1.5 +
     1.6 +/***********************************************************************
     1.7 +*  This code is part of GLPK (GNU Linear Programming Kit).
     1.8 +*
     1.9 +*  Copyright (C) 2000, 2001, 2002, 2003, 2004, 2005, 2006, 2007, 2008,
    1.10 +*  2009, 2010, 2011 Andrew Makhorin, Department for Applied Informatics,
    1.11 +*  Moscow Aviation Institute, Moscow, Russia. All rights reserved.
    1.12 +*  E-mail: <mao@gnu.org>.
    1.13 +*
    1.14 +*  GLPK is free software: you can redistribute it and/or modify it
    1.15 +*  under the terms of the GNU General Public License as published by
    1.16 +*  the Free Software Foundation, either version 3 of the License, or
    1.17 +*  (at your option) any later version.
    1.18 +*
    1.19 +*  GLPK is distributed in the hope that it will be useful, but WITHOUT
    1.20 +*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
    1.21 +*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
    1.22 +*  License for more details.
    1.23 +*
    1.24 +*  You should have received a copy of the GNU General Public License
    1.25 +*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
    1.26 +***********************************************************************/
    1.27 +
    1.28 +#include "glpipm.h"
    1.29 +#include "glpmat.h"
    1.30 +
    1.31 +#define ITER_MAX 100
    1.32 +/* maximal number of iterations */
    1.33 +
    1.34 +struct csa
    1.35 +{     /* common storage area */
    1.36 +      /*--------------------------------------------------------------*/
    1.37 +      /* LP data */
    1.38 +      int m;
    1.39 +      /* number of rows (equality constraints) */
    1.40 +      int n;
    1.41 +      /* number of columns (structural variables) */
    1.42 +      int *A_ptr; /* int A_ptr[1+m+1]; */
    1.43 +      int *A_ind; /* int A_ind[A_ptr[m+1]]; */
    1.44 +      double *A_val; /* double A_val[A_ptr[m+1]]; */
    1.45 +      /* mxn-matrix A in storage-by-rows format */
    1.46 +      double *b; /* double b[1+m]; */
    1.47 +      /* m-vector b of right-hand sides */
    1.48 +      double *c; /* double c[1+n]; */
    1.49 +      /* n-vector c of objective coefficients; c[0] is constant term of
    1.50 +         the objective function */
    1.51 +      /*--------------------------------------------------------------*/
    1.52 +      /* LP solution */
    1.53 +      double *x; /* double x[1+n]; */
    1.54 +      double *y; /* double y[1+m]; */
    1.55 +      double *z; /* double z[1+n]; */
    1.56 +      /* current point in primal-dual space; the best point on exit */
    1.57 +      /*--------------------------------------------------------------*/
    1.58 +      /* control parameters */
    1.59 +      const glp_iptcp *parm;
    1.60 +      /*--------------------------------------------------------------*/
    1.61 +      /* working arrays and variables */
    1.62 +      double *D; /* double D[1+n]; */
    1.63 +      /* diagonal nxn-matrix D = X*inv(Z), where X = diag(x[j]) and
    1.64 +         Z = diag(z[j]) */
    1.65 +      int *P; /* int P[1+m+m]; */
    1.66 +      /* permutation mxm-matrix P used to minimize fill-in in Cholesky
    1.67 +         factorization */
    1.68 +      int *S_ptr; /* int S_ptr[1+m+1]; */
    1.69 +      int *S_ind; /* int S_ind[S_ptr[m+1]]; */
    1.70 +      double *S_val; /* double S_val[S_ptr[m+1]]; */
    1.71 +      double *S_diag; /* double S_diag[1+m]; */
    1.72 +      /* symmetric mxm-matrix S = P*A*D*A'*P' whose upper triangular
    1.73 +         part without diagonal elements is stored in S_ptr, S_ind, and
    1.74 +         S_val in storage-by-rows format, diagonal elements are stored
    1.75 +         in S_diag */
    1.76 +      int *U_ptr; /* int U_ptr[1+m+1]; */
    1.77 +      int *U_ind; /* int U_ind[U_ptr[m+1]]; */
    1.78 +      double *U_val; /* double U_val[U_ptr[m+1]]; */
    1.79 +      double *U_diag; /* double U_diag[1+m]; */
    1.80 +      /* upper triangular mxm-matrix U defining Cholesky factorization
    1.81 +         S = U'*U; its non-diagonal elements are stored in U_ptr, U_ind,
    1.82 +         U_val in storage-by-rows format, diagonal elements are stored
    1.83 +         in U_diag */
    1.84 +      int iter;
    1.85 +      /* iteration number (0, 1, 2, ...); iter = 0 corresponds to the
    1.86 +         initial point */
    1.87 +      double obj;
    1.88 +      /* current value of the objective function */
    1.89 +      double rpi;
    1.90 +      /* relative primal infeasibility rpi = ||A*x-b||/(1+||b||) */
    1.91 +      double rdi;
    1.92 +      /* relative dual infeasibility rdi = ||A'*y+z-c||/(1+||c||) */
    1.93 +      double gap;
    1.94 +      /* primal-dual gap = |c'*x-b'*y|/(1+|c'*x|) which is a relative
    1.95 +         difference between primal and dual objective functions */
    1.96 +      double phi;
    1.97 +      /* merit function phi = ||A*x-b||/max(1,||b||) +
    1.98 +                            + ||A'*y+z-c||/max(1,||c||) +
    1.99 +                            + |c'*x-b'*y|/max(1,||b||,||c||) */
   1.100 +      double mu;
   1.101 +      /* duality measure mu = x'*z/n (used as barrier parameter) */
   1.102 +      double rmu;
   1.103 +      /* rmu = max(||A*x-b||,||A'*y+z-c||)/mu */
   1.104 +      double rmu0;
   1.105 +      /* the initial value of rmu on iteration 0 */
   1.106 +      double *phi_min; /* double phi_min[1+ITER_MAX]; */
   1.107 +      /* phi_min[k] = min(phi[k]), where phi[k] is the value of phi on
   1.108 +         k-th iteration, 0 <= k <= iter */
   1.109 +      int best_iter;
   1.110 +      /* iteration number, on which the value of phi reached its best
   1.111 +         (minimal) value */
   1.112 +      double *best_x; /* double best_x[1+n]; */
   1.113 +      double *best_y; /* double best_y[1+m]; */
   1.114 +      double *best_z; /* double best_z[1+n]; */
   1.115 +      /* best point (in the sense of the merit function phi) which has
   1.116 +         been reached on iteration iter_best */
   1.117 +      double best_obj;
   1.118 +      /* objective value at the best point */
   1.119 +      double *dx_aff; /* double dx_aff[1+n]; */
   1.120 +      double *dy_aff; /* double dy_aff[1+m]; */
   1.121 +      double *dz_aff; /* double dz_aff[1+n]; */
   1.122 +      /* affine scaling direction */
   1.123 +      double alfa_aff_p, alfa_aff_d;
   1.124 +      /* maximal primal and dual stepsizes in affine scaling direction,
   1.125 +         on which x and z are still non-negative */
   1.126 +      double mu_aff;
   1.127 +      /* duality measure mu_aff = x_aff'*z_aff/n in the boundary point
   1.128 +         x_aff' = x+alfa_aff_p*dx_aff, z_aff' = z+alfa_aff_d*dz_aff */
   1.129 +      double sigma;
   1.130 +      /* Mehrotra's heuristic parameter (0 <= sigma <= 1) */
   1.131 +      double *dx_cc; /* double dx_cc[1+n]; */
   1.132 +      double *dy_cc; /* double dy_cc[1+m]; */
   1.133 +      double *dz_cc; /* double dz_cc[1+n]; */
   1.134 +      /* centering corrector direction */
   1.135 +      double *dx; /* double dx[1+n]; */
   1.136 +      double *dy; /* double dy[1+m]; */
   1.137 +      double *dz; /* double dz[1+n]; */
   1.138 +      /* final combined direction dx = dx_aff+dx_cc, dy = dy_aff+dy_cc,
   1.139 +         dz = dz_aff+dz_cc */
   1.140 +      double alfa_max_p;
   1.141 +      double alfa_max_d;
   1.142 +      /* maximal primal and dual stepsizes in combined direction, on
   1.143 +         which x and z are still non-negative */
   1.144 +};
   1.145 +
   1.146 +/***********************************************************************
   1.147 +*  initialize - allocate and initialize common storage area
   1.148 +*
   1.149 +*  This routine allocates and initializes the common storage area (CSA)
   1.150 +*  used by interior-point method routines. */
   1.151 +
   1.152 +static void initialize(struct csa *csa)
   1.153 +{     int m = csa->m;
   1.154 +      int n = csa->n;
   1.155 +      int i;
   1.156 +      if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.157 +         xprintf("Matrix A has %d non-zeros\n", csa->A_ptr[m+1]-1);
   1.158 +      csa->D = xcalloc(1+n, sizeof(double));
   1.159 +      /* P := I */
   1.160 +      csa->P = xcalloc(1+m+m, sizeof(int));
   1.161 +      for (i = 1; i <= m; i++) csa->P[i] = csa->P[m+i] = i;
   1.162 +      /* S := A*A', symbolically */
   1.163 +      csa->S_ptr = xcalloc(1+m+1, sizeof(int));
   1.164 +      csa->S_ind = adat_symbolic(m, n, csa->P, csa->A_ptr, csa->A_ind,
   1.165 +         csa->S_ptr);
   1.166 +      if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.167 +         xprintf("Matrix S = A*A' has %d non-zeros (upper triangle)\n",
   1.168 +            csa->S_ptr[m+1]-1 + m);
   1.169 +      /* determine P using specified ordering algorithm */
   1.170 +      if (csa->parm->ord_alg == GLP_ORD_NONE)
   1.171 +      {  if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.172 +            xprintf("Original ordering is being used\n");
   1.173 +         for (i = 1; i <= m; i++)
   1.174 +            csa->P[i] = csa->P[m+i] = i;
   1.175 +      }
   1.176 +      else if (csa->parm->ord_alg == GLP_ORD_QMD)
   1.177 +      {  if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.178 +            xprintf("Minimum degree ordering (QMD)...\n");
   1.179 +         min_degree(m, csa->S_ptr, csa->S_ind, csa->P);
   1.180 +      }
   1.181 +      else if (csa->parm->ord_alg == GLP_ORD_AMD)
   1.182 +      {  if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.183 +            xprintf("Approximate minimum degree ordering (AMD)...\n");
   1.184 +         amd_order1(m, csa->S_ptr, csa->S_ind, csa->P);
   1.185 +      }
   1.186 +      else if (csa->parm->ord_alg == GLP_ORD_SYMAMD)
   1.187 +      {  if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.188 +            xprintf("Approximate minimum degree ordering (SYMAMD)...\n")
   1.189 +               ;
   1.190 +         symamd_ord(m, csa->S_ptr, csa->S_ind, csa->P);
   1.191 +      }
   1.192 +      else
   1.193 +         xassert(csa != csa);
   1.194 +      /* S := P*A*A'*P', symbolically */
   1.195 +      xfree(csa->S_ind);
   1.196 +      csa->S_ind = adat_symbolic(m, n, csa->P, csa->A_ptr, csa->A_ind,
   1.197 +         csa->S_ptr);
   1.198 +      csa->S_val = xcalloc(csa->S_ptr[m+1], sizeof(double));
   1.199 +      csa->S_diag = xcalloc(1+m, sizeof(double));
   1.200 +      /* compute Cholesky factorization S = U'*U, symbolically */
   1.201 +      if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.202 +         xprintf("Computing Cholesky factorization S = L*L'...\n");
   1.203 +      csa->U_ptr = xcalloc(1+m+1, sizeof(int));
   1.204 +      csa->U_ind = chol_symbolic(m, csa->S_ptr, csa->S_ind, csa->U_ptr);
   1.205 +      if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.206 +         xprintf("Matrix L has %d non-zeros\n", csa->U_ptr[m+1]-1 + m);
   1.207 +      csa->U_val = xcalloc(csa->U_ptr[m+1], sizeof(double));
   1.208 +      csa->U_diag = xcalloc(1+m, sizeof(double));
   1.209 +      csa->iter = 0;
   1.210 +      csa->obj = 0.0;
   1.211 +      csa->rpi = 0.0;
   1.212 +      csa->rdi = 0.0;
   1.213 +      csa->gap = 0.0;
   1.214 +      csa->phi = 0.0;
   1.215 +      csa->mu = 0.0;
   1.216 +      csa->rmu = 0.0;
   1.217 +      csa->rmu0 = 0.0;
   1.218 +      csa->phi_min = xcalloc(1+ITER_MAX, sizeof(double));
   1.219 +      csa->best_iter = 0;
   1.220 +      csa->best_x = xcalloc(1+n, sizeof(double));
   1.221 +      csa->best_y = xcalloc(1+m, sizeof(double));
   1.222 +      csa->best_z = xcalloc(1+n, sizeof(double));
   1.223 +      csa->best_obj = 0.0;
   1.224 +      csa->dx_aff = xcalloc(1+n, sizeof(double));
   1.225 +      csa->dy_aff = xcalloc(1+m, sizeof(double));
   1.226 +      csa->dz_aff = xcalloc(1+n, sizeof(double));
   1.227 +      csa->alfa_aff_p = 0.0;
   1.228 +      csa->alfa_aff_d = 0.0;
   1.229 +      csa->mu_aff = 0.0;
   1.230 +      csa->sigma = 0.0;
   1.231 +      csa->dx_cc = xcalloc(1+n, sizeof(double));
   1.232 +      csa->dy_cc = xcalloc(1+m, sizeof(double));
   1.233 +      csa->dz_cc = xcalloc(1+n, sizeof(double));
   1.234 +      csa->dx = csa->dx_aff;
   1.235 +      csa->dy = csa->dy_aff;
   1.236 +      csa->dz = csa->dz_aff;
   1.237 +      csa->alfa_max_p = 0.0;
   1.238 +      csa->alfa_max_d = 0.0;
   1.239 +      return;
   1.240 +}
   1.241 +
   1.242 +/***********************************************************************
   1.243 +*  A_by_vec - compute y = A*x
   1.244 +*
   1.245 +*  This routine computes matrix-vector product y = A*x, where A is the
   1.246 +*  constraint matrix. */
   1.247 +
   1.248 +static void A_by_vec(struct csa *csa, double x[], double y[])
   1.249 +{     /* compute y = A*x */
   1.250 +      int m = csa->m;
   1.251 +      int *A_ptr = csa->A_ptr;
   1.252 +      int *A_ind = csa->A_ind;
   1.253 +      double *A_val = csa->A_val;
   1.254 +      int i, t, beg, end;
   1.255 +      double temp;
   1.256 +      for (i = 1; i <= m; i++)
   1.257 +      {  temp = 0.0;
   1.258 +         beg = A_ptr[i], end = A_ptr[i+1];
   1.259 +         for (t = beg; t < end; t++) temp += A_val[t] * x[A_ind[t]];
   1.260 +         y[i] = temp;
   1.261 +      }
   1.262 +      return;
   1.263 +}
   1.264 +
   1.265 +/***********************************************************************
   1.266 +*  AT_by_vec - compute y = A'*x
   1.267 +*
   1.268 +*  This routine computes matrix-vector product y = A'*x, where A' is a
   1.269 +*  matrix transposed to the constraint matrix A. */
   1.270 +
   1.271 +static void AT_by_vec(struct csa *csa, double x[], double y[])
   1.272 +{     /* compute y = A'*x, where A' is transposed to A */
   1.273 +      int m = csa->m;
   1.274 +      int n = csa->n;
   1.275 +      int *A_ptr = csa->A_ptr;
   1.276 +      int *A_ind = csa->A_ind;
   1.277 +      double *A_val = csa->A_val;
   1.278 +      int i, j, t, beg, end;
   1.279 +      double temp;
   1.280 +      for (j = 1; j <= n; j++) y[j] = 0.0;
   1.281 +      for (i = 1; i <= m; i++)
   1.282 +      {  temp = x[i];
   1.283 +         if (temp == 0.0) continue;
   1.284 +         beg = A_ptr[i], end = A_ptr[i+1];
   1.285 +         for (t = beg; t < end; t++) y[A_ind[t]] += A_val[t] * temp;
   1.286 +      }
   1.287 +      return;
   1.288 +}
   1.289 +
   1.290 +/***********************************************************************
   1.291 +*  decomp_NE - numeric factorization of matrix S = P*A*D*A'*P'
   1.292 +*
   1.293 +*  This routine implements numeric phase of Cholesky factorization of
   1.294 +*  the matrix S = P*A*D*A'*P', which is a permuted matrix of the normal
   1.295 +*  equation system. Matrix D is assumed to be already computed. */
   1.296 +
   1.297 +static void decomp_NE(struct csa *csa)
   1.298 +{     adat_numeric(csa->m, csa->n, csa->P, csa->A_ptr, csa->A_ind,
   1.299 +         csa->A_val, csa->D, csa->S_ptr, csa->S_ind, csa->S_val,
   1.300 +         csa->S_diag);
   1.301 +      chol_numeric(csa->m, csa->S_ptr, csa->S_ind, csa->S_val,
   1.302 +         csa->S_diag, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag);
   1.303 +      return;
   1.304 +}
   1.305 +
   1.306 +/***********************************************************************
   1.307 +*  solve_NE - solve normal equation system
   1.308 +*
   1.309 +*  This routine solves the normal equation system:
   1.310 +*
   1.311 +*     A*D*A'*y = h.
   1.312 +*
   1.313 +*  It is assumed that the matrix A*D*A' has been previously factorized
   1.314 +*  by the routine decomp_NE.
   1.315 +*
   1.316 +*  On entry the array y contains the vector of right-hand sides h. On
   1.317 +*  exit this array contains the computed vector of unknowns y.
   1.318 +*
   1.319 +*  Once the vector y has been computed the routine checks for numeric
   1.320 +*  stability. If the residual vector:
   1.321 +*
   1.322 +*     r = A*D*A'*y - h
   1.323 +*
   1.324 +*  is relatively small, the routine returns zero, otherwise non-zero is
   1.325 +*  returned. */
   1.326 +
   1.327 +static int solve_NE(struct csa *csa, double y[])
   1.328 +{     int m = csa->m;
   1.329 +      int n = csa->n;
   1.330 +      int *P = csa->P;
   1.331 +      int i, j, ret = 0;
   1.332 +      double *h, *r, *w;
   1.333 +      /* save vector of right-hand sides h */
   1.334 +      h = xcalloc(1+m, sizeof(double));
   1.335 +      for (i = 1; i <= m; i++) h[i] = y[i];
   1.336 +      /* solve normal equation system (A*D*A')*y = h */
   1.337 +      /* since S = P*A*D*A'*P' = U'*U, then A*D*A' = P'*U'*U*P, so we
   1.338 +         have inv(A*D*A') = P'*inv(U)*inv(U')*P */
   1.339 +      /* w := P*h */
   1.340 +      w = xcalloc(1+m, sizeof(double));
   1.341 +      for (i = 1; i <= m; i++) w[i] = y[P[i]];
   1.342 +      /* w := inv(U')*w */
   1.343 +      ut_solve(m, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag, w);
   1.344 +      /* w := inv(U)*w */
   1.345 +      u_solve(m, csa->U_ptr, csa->U_ind, csa->U_val, csa->U_diag, w);
   1.346 +      /* y := P'*w */
   1.347 +      for (i = 1; i <= m; i++) y[i] = w[P[m+i]];
   1.348 +      xfree(w);
   1.349 +      /* compute residual vector r = A*D*A'*y - h */
   1.350 +      r = xcalloc(1+m, sizeof(double));
   1.351 +      /* w := A'*y */
   1.352 +      w = xcalloc(1+n, sizeof(double));
   1.353 +      AT_by_vec(csa, y, w);
   1.354 +      /* w := D*w */
   1.355 +      for (j = 1; j <= n; j++) w[j] *= csa->D[j];
   1.356 +      /* r := A*w */
   1.357 +      A_by_vec(csa, w, r);
   1.358 +      xfree(w);
   1.359 +      /* r := r - h */
   1.360 +      for (i = 1; i <= m; i++) r[i] -= h[i];
   1.361 +      /* check for numeric stability */
   1.362 +      for (i = 1; i <= m; i++)
   1.363 +      {  if (fabs(r[i]) / (1.0 + fabs(h[i])) > 1e-4)
   1.364 +         {  ret = 1;
   1.365 +            break;
   1.366 +         }
   1.367 +      }
   1.368 +      xfree(h);
   1.369 +      xfree(r);
   1.370 +      return ret;
   1.371 +}
   1.372 +
   1.373 +/***********************************************************************
   1.374 +*  solve_NS - solve Newtonian system
   1.375 +*
   1.376 +*  This routine solves the Newtonian system:
   1.377 +*
   1.378 +*     A*dx               = p
   1.379 +*
   1.380 +*           A'*dy +   dz = q
   1.381 +*
   1.382 +*     Z*dx        + X*dz = r
   1.383 +*
   1.384 +*  where X = diag(x[j]), Z = diag(z[j]), by reducing it to the normal
   1.385 +*  equation system:
   1.386 +*
   1.387 +*     (A*inv(Z)*X*A')*dy = A*inv(Z)*(X*q-r)+p
   1.388 +*
   1.389 +*  (it is assumed that the matrix A*inv(Z)*X*A' has been factorized by
   1.390 +*  the routine decomp_NE).
   1.391 +*
   1.392 +*  Once vector dy has been computed the routine computes vectors dx and
   1.393 +*  dz as follows:
   1.394 +*
   1.395 +*     dx = inv(Z)*(X*(A'*dy-q)+r)
   1.396 +*
   1.397 +*     dz = inv(X)*(r-Z*dx)
   1.398 +*
   1.399 +*  The routine solve_NS returns the same code which was reported by the
   1.400 +*  routine solve_NE (see above). */
   1.401 +
   1.402 +static int solve_NS(struct csa *csa, double p[], double q[], double r[],
   1.403 +      double dx[], double dy[], double dz[])
   1.404 +{     int m = csa->m;
   1.405 +      int n = csa->n;
   1.406 +      double *x = csa->x;
   1.407 +      double *z = csa->z;
   1.408 +      int i, j, ret;
   1.409 +      double *w = dx;
   1.410 +      /* compute the vector of right-hand sides A*inv(Z)*(X*q-r)+p for
   1.411 +         the normal equation system */
   1.412 +      for (j = 1; j <= n; j++)
   1.413 +         w[j] = (x[j] * q[j] - r[j]) / z[j];
   1.414 +      A_by_vec(csa, w, dy);
   1.415 +      for (i = 1; i <= m; i++) dy[i] += p[i];
   1.416 +      /* solve the normal equation system to compute vector dy */
   1.417 +      ret = solve_NE(csa, dy);
   1.418 +      /* compute vectors dx and dz */
   1.419 +      AT_by_vec(csa, dy, dx);
   1.420 +      for (j = 1; j <= n; j++)
   1.421 +      {  dx[j] = (x[j] * (dx[j] - q[j]) + r[j]) / z[j];
   1.422 +         dz[j] = (r[j] - z[j] * dx[j]) / x[j];
   1.423 +      }
   1.424 +      return ret;
   1.425 +}
   1.426 +
   1.427 +/***********************************************************************
   1.428 +*  initial_point - choose initial point using Mehrotra's heuristic
   1.429 +*
   1.430 +*  This routine chooses a starting point using a heuristic proposed in
   1.431 +*  the paper:
   1.432 +*
   1.433 +*  S. Mehrotra. On the implementation of a primal-dual interior point
   1.434 +*  method. SIAM J. on Optim., 2(4), pp. 575-601, 1992.
   1.435 +*
   1.436 +*  The starting point x in the primal space is chosen as a solution of
   1.437 +*  the following least squares problem:
   1.438 +*
   1.439 +*     minimize    ||x||
   1.440 +*
   1.441 +*     subject to  A*x = b
   1.442 +*
   1.443 +*  which can be computed explicitly as follows:
   1.444 +*
   1.445 +*     x = A'*inv(A*A')*b
   1.446 +*
   1.447 +*  Similarly, the starting point (y, z) in the dual space is chosen as
   1.448 +*  a solution of the following least squares problem:
   1.449 +*
   1.450 +*     minimize    ||z||
   1.451 +*
   1.452 +*     subject to  A'*y + z = c
   1.453 +*
   1.454 +*  which can be computed explicitly as follows:
   1.455 +*
   1.456 +*     y = inv(A*A')*A*c
   1.457 +*
   1.458 +*     z = c - A'*y
   1.459 +*
   1.460 +*  However, some components of the vectors x and z may be non-positive
   1.461 +*  or close to zero, so the routine uses a Mehrotra's heuristic to find
   1.462 +*  a more appropriate starting point. */
   1.463 +
   1.464 +static void initial_point(struct csa *csa)
   1.465 +{     int m = csa->m;
   1.466 +      int n = csa->n;
   1.467 +      double *b = csa->b;
   1.468 +      double *c = csa->c;
   1.469 +      double *x = csa->x;
   1.470 +      double *y = csa->y;
   1.471 +      double *z = csa->z;
   1.472 +      double *D = csa->D;
   1.473 +      int i, j;
   1.474 +      double dp, dd, ex, ez, xz;
   1.475 +      /* factorize A*A' */
   1.476 +      for (j = 1; j <= n; j++) D[j] = 1.0;
   1.477 +      decomp_NE(csa);
   1.478 +      /* x~ = A'*inv(A*A')*b */
   1.479 +      for (i = 1; i <= m; i++) y[i] = b[i];
   1.480 +      solve_NE(csa, y);
   1.481 +      AT_by_vec(csa, y, x);
   1.482 +      /* y~ = inv(A*A')*A*c */
   1.483 +      A_by_vec(csa, c, y);
   1.484 +      solve_NE(csa, y);
   1.485 +      /* z~ = c - A'*y~ */
   1.486 +      AT_by_vec(csa, y,z);
   1.487 +      for (j = 1; j <= n; j++) z[j] = c[j] - z[j];
   1.488 +      /* use Mehrotra's heuristic in order to choose more appropriate
   1.489 +         starting point with positive components of vectors x and z */
   1.490 +      dp = dd = 0.0;
   1.491 +      for (j = 1; j <= n; j++)
   1.492 +      {  if (dp < -1.5 * x[j]) dp = -1.5 * x[j];
   1.493 +         if (dd < -1.5 * z[j]) dd = -1.5 * z[j];
   1.494 +      }
   1.495 +      /* note that b = 0 involves x = 0, and c = 0 involves y = 0 and
   1.496 +         z = 0, so we need to be careful */
   1.497 +      if (dp == 0.0) dp = 1.5;
   1.498 +      if (dd == 0.0) dd = 1.5;
   1.499 +      ex = ez = xz = 0.0;
   1.500 +      for (j = 1; j <= n; j++)
   1.501 +      {  ex += (x[j] + dp);
   1.502 +         ez += (z[j] + dd);
   1.503 +         xz += (x[j] + dp) * (z[j] + dd);
   1.504 +      }
   1.505 +      dp += 0.5 * (xz / ez);
   1.506 +      dd += 0.5 * (xz / ex);
   1.507 +      for (j = 1; j <= n; j++)
   1.508 +      {  x[j] += dp;
   1.509 +         z[j] += dd;
   1.510 +         xassert(x[j] > 0.0 && z[j] > 0.0);
   1.511 +      }
   1.512 +      return;
   1.513 +}
   1.514 +
   1.515 +/***********************************************************************
   1.516 +*  basic_info - perform basic computations at the current point
   1.517 +*
   1.518 +*  This routine computes the following quantities at the current point:
   1.519 +*
   1.520 +*  1) value of the objective function:
   1.521 +*
   1.522 +*     F = c'*x + c[0]
   1.523 +*
   1.524 +*  2) relative primal infeasibility:
   1.525 +*
   1.526 +*     rpi = ||A*x-b|| / (1+||b||)
   1.527 +*
   1.528 +*  3) relative dual infeasibility:
   1.529 +*
   1.530 +*     rdi = ||A'*y+z-c|| / (1+||c||)
   1.531 +*
   1.532 +*  4) primal-dual gap (relative difference between the primal and the
   1.533 +*     dual objective function values):
   1.534 +*
   1.535 +*     gap = |c'*x-b'*y| / (1+|c'*x|)
   1.536 +*
   1.537 +*  5) merit function:
   1.538 +*
   1.539 +*     phi = ||A*x-b|| / max(1,||b||) + ||A'*y+z-c|| / max(1,||c||) +
   1.540 +*
   1.541 +*         + |c'*x-b'*y| / max(1,||b||,||c||)
   1.542 +*
   1.543 +*  6) duality measure:
   1.544 +*
   1.545 +*     mu = x'*z / n
   1.546 +*
   1.547 +*  7) the ratio of infeasibility to mu:
   1.548 +*
   1.549 +*     rmu = max(||A*x-b||,||A'*y+z-c||) / mu
   1.550 +*
   1.551 +*  where ||*|| denotes euclidian norm, *' denotes transposition. */
   1.552 +
   1.553 +static void basic_info(struct csa *csa)
   1.554 +{     int m = csa->m;
   1.555 +      int n = csa->n;
   1.556 +      double *b = csa->b;
   1.557 +      double *c = csa->c;
   1.558 +      double *x = csa->x;
   1.559 +      double *y = csa->y;
   1.560 +      double *z = csa->z;
   1.561 +      int i, j;
   1.562 +      double norm1, bnorm, norm2, cnorm, cx, by, *work, temp;
   1.563 +      /* compute value of the objective function */
   1.564 +      temp = c[0];
   1.565 +      for (j = 1; j <= n; j++) temp += c[j] * x[j];
   1.566 +      csa->obj = temp;
   1.567 +      /* norm1 = ||A*x-b|| */
   1.568 +      work = xcalloc(1+m, sizeof(double));
   1.569 +      A_by_vec(csa, x, work);
   1.570 +      norm1 = 0.0;
   1.571 +      for (i = 1; i <= m; i++)
   1.572 +         norm1 += (work[i] - b[i]) * (work[i] - b[i]);
   1.573 +      norm1 = sqrt(norm1);
   1.574 +      xfree(work);
   1.575 +      /* bnorm = ||b|| */
   1.576 +      bnorm = 0.0;
   1.577 +      for (i = 1; i <= m; i++) bnorm += b[i] * b[i];
   1.578 +      bnorm = sqrt(bnorm);
   1.579 +      /* compute relative primal infeasibility */
   1.580 +      csa->rpi = norm1 / (1.0 + bnorm);
   1.581 +      /* norm2 = ||A'*y+z-c|| */
   1.582 +      work = xcalloc(1+n, sizeof(double));
   1.583 +      AT_by_vec(csa, y, work);
   1.584 +      norm2 = 0.0;
   1.585 +      for (j = 1; j <= n; j++)
   1.586 +         norm2 += (work[j] + z[j] - c[j]) * (work[j] + z[j] - c[j]);
   1.587 +      norm2 = sqrt(norm2);
   1.588 +      xfree(work);
   1.589 +      /* cnorm = ||c|| */
   1.590 +      cnorm = 0.0;
   1.591 +      for (j = 1; j <= n; j++) cnorm += c[j] * c[j];
   1.592 +      cnorm = sqrt(cnorm);
   1.593 +      /* compute relative dual infeasibility */
   1.594 +      csa->rdi = norm2 / (1.0 + cnorm);
   1.595 +      /* by = b'*y */
   1.596 +      by = 0.0;
   1.597 +      for (i = 1; i <= m; i++) by += b[i] * y[i];
   1.598 +      /* cx = c'*x */
   1.599 +      cx = 0.0;
   1.600 +      for (j = 1; j <= n; j++) cx += c[j] * x[j];
   1.601 +      /* compute primal-dual gap */
   1.602 +      csa->gap = fabs(cx - by) / (1.0 + fabs(cx));
   1.603 +      /* compute merit function */
   1.604 +      csa->phi = 0.0;
   1.605 +      csa->phi += norm1 / (bnorm > 1.0 ? bnorm : 1.0);
   1.606 +      csa->phi += norm2 / (cnorm > 1.0 ? cnorm : 1.0);
   1.607 +      temp = 1.0;
   1.608 +      if (temp < bnorm) temp = bnorm;
   1.609 +      if (temp < cnorm) temp = cnorm;
   1.610 +      csa->phi += fabs(cx - by) / temp;
   1.611 +      /* compute duality measure */
   1.612 +      temp = 0.0;
   1.613 +      for (j = 1; j <= n; j++) temp += x[j] * z[j];
   1.614 +      csa->mu = temp / (double)n;
   1.615 +      /* compute the ratio of infeasibility to mu */
   1.616 +      csa->rmu = (norm1 > norm2 ? norm1 : norm2) / csa->mu;
   1.617 +      return;
   1.618 +}
   1.619 +
   1.620 +/***********************************************************************
   1.621 +*  make_step - compute next point using Mehrotra's technique
   1.622 +*
   1.623 +*  This routine computes the next point using the predictor-corrector
   1.624 +*  technique proposed in the paper:
   1.625 +*
   1.626 +*  S. Mehrotra. On the implementation of a primal-dual interior point
   1.627 +*  method. SIAM J. on Optim., 2(4), pp. 575-601, 1992.
   1.628 +*
   1.629 +*  At first, the routine computes so called affine scaling (predictor)
   1.630 +*  direction (dx_aff,dy_aff,dz_aff) which is a solution of the system:
   1.631 +*
   1.632 +*     A*dx_aff                       = b - A*x
   1.633 +*
   1.634 +*               A'*dy_aff +   dz_aff = c - A'*y - z
   1.635 +*
   1.636 +*     Z*dx_aff            + X*dz_aff = - X*Z*e
   1.637 +*
   1.638 +*  where (x,y,z) is the current point, X = diag(x[j]), Z = diag(z[j]),
   1.639 +*  e = (1,...,1)'.
   1.640 +*
   1.641 +*  Then, the routine computes the centering parameter sigma, using the
   1.642 +*  following Mehrotra's heuristic:
   1.643 +*
   1.644 +*     alfa_aff_p = inf{0 <= alfa <= 1 | x+alfa*dx_aff >= 0}
   1.645 +*
   1.646 +*     alfa_aff_d = inf{0 <= alfa <= 1 | z+alfa*dz_aff >= 0}
   1.647 +*
   1.648 +*     mu_aff = (x+alfa_aff_p*dx_aff)'*(z+alfa_aff_d*dz_aff)/n
   1.649 +*
   1.650 +*     sigma = (mu_aff/mu)^3
   1.651 +*
   1.652 +*  where alfa_aff_p is the maximal stepsize along the affine scaling
   1.653 +*  direction in the primal space, alfa_aff_d is the maximal stepsize
   1.654 +*  along the same direction in the dual space.
   1.655 +*
   1.656 +*  After determining sigma the routine computes so called centering
   1.657 +*  (corrector) direction (dx_cc,dy_cc,dz_cc) which is the solution of
   1.658 +*  the system:
   1.659 +*
   1.660 +*     A*dx_cc                     = 0
   1.661 +*
   1.662 +*              A'*dy_cc +   dz_cc = 0
   1.663 +*
   1.664 +*     Z*dx_cc           + X*dz_cc = sigma*mu*e - X*Z*e
   1.665 +*
   1.666 +*  Finally, the routine computes the combined direction
   1.667 +*
   1.668 +*     (dx,dy,dz) = (dx_aff,dy_aff,dz_aff) + (dx_cc,dy_cc,dz_cc)
   1.669 +*
   1.670 +*  and determines maximal primal and dual stepsizes along the combined
   1.671 +*  direction:
   1.672 +*
   1.673 +*     alfa_max_p = inf{0 <= alfa <= 1 | x+alfa*dx >= 0}
   1.674 +*
   1.675 +*     alfa_max_d = inf{0 <= alfa <= 1 | z+alfa*dz >= 0}
   1.676 +*
   1.677 +*  In order to prevent the next point to be too close to the boundary
   1.678 +*  of the positive ortant, the routine decreases maximal stepsizes:
   1.679 +*
   1.680 +*     alfa_p = gamma_p * alfa_max_p
   1.681 +*
   1.682 +*     alfa_d = gamma_d * alfa_max_d
   1.683 +*
   1.684 +*  where gamma_p and gamma_d are scaling factors, and computes the next
   1.685 +*  point:
   1.686 +*
   1.687 +*     x_new = x + alfa_p * dx
   1.688 +*
   1.689 +*     y_new = y + alfa_d * dy
   1.690 +*
   1.691 +*     z_new = z + alfa_d * dz
   1.692 +*
   1.693 +*  which becomes the current point on the next iteration. */
   1.694 +
   1.695 +static int make_step(struct csa *csa)
   1.696 +{     int m = csa->m;
   1.697 +      int n = csa->n;
   1.698 +      double *b = csa->b;
   1.699 +      double *c = csa->c;
   1.700 +      double *x = csa->x;
   1.701 +      double *y = csa->y;
   1.702 +      double *z = csa->z;
   1.703 +      double *dx_aff = csa->dx_aff;
   1.704 +      double *dy_aff = csa->dy_aff;
   1.705 +      double *dz_aff = csa->dz_aff;
   1.706 +      double *dx_cc = csa->dx_cc;
   1.707 +      double *dy_cc = csa->dy_cc;
   1.708 +      double *dz_cc = csa->dz_cc;
   1.709 +      double *dx = csa->dx;
   1.710 +      double *dy = csa->dy;
   1.711 +      double *dz = csa->dz;
   1.712 +      int i, j, ret = 0;
   1.713 +      double temp, gamma_p, gamma_d, *p, *q, *r;
   1.714 +      /* allocate working arrays */
   1.715 +      p = xcalloc(1+m, sizeof(double));
   1.716 +      q = xcalloc(1+n, sizeof(double));
   1.717 +      r = xcalloc(1+n, sizeof(double));
   1.718 +      /* p = b - A*x */
   1.719 +      A_by_vec(csa, x, p);
   1.720 +      for (i = 1; i <= m; i++) p[i] = b[i] - p[i];
   1.721 +      /* q = c - A'*y - z */
   1.722 +      AT_by_vec(csa, y,q);
   1.723 +      for (j = 1; j <= n; j++) q[j] = c[j] - q[j] - z[j];
   1.724 +      /* r = - X * Z * e */
   1.725 +      for (j = 1; j <= n; j++) r[j] = - x[j] * z[j];
   1.726 +      /* solve the first Newtonian system */
   1.727 +      if (solve_NS(csa, p, q, r, dx_aff, dy_aff, dz_aff))
   1.728 +      {  ret = 1;
   1.729 +         goto done;
   1.730 +      }
   1.731 +      /* alfa_aff_p = inf{0 <= alfa <= 1 | x + alfa*dx_aff >= 0} */
   1.732 +      /* alfa_aff_d = inf{0 <= alfa <= 1 | z + alfa*dz_aff >= 0} */
   1.733 +      csa->alfa_aff_p = csa->alfa_aff_d = 1.0;
   1.734 +      for (j = 1; j <= n; j++)
   1.735 +      {  if (dx_aff[j] < 0.0)
   1.736 +         {  temp = - x[j] / dx_aff[j];
   1.737 +            if (csa->alfa_aff_p > temp) csa->alfa_aff_p = temp;
   1.738 +         }
   1.739 +         if (dz_aff[j] < 0.0)
   1.740 +         {  temp = - z[j] / dz_aff[j];
   1.741 +            if (csa->alfa_aff_d > temp) csa->alfa_aff_d = temp;
   1.742 +         }
   1.743 +      }
   1.744 +      /* mu_aff = (x+alfa_aff_p*dx_aff)' * (z+alfa_aff_d*dz_aff) / n */
   1.745 +      temp = 0.0;
   1.746 +      for (j = 1; j <= n; j++)
   1.747 +         temp += (x[j] + csa->alfa_aff_p * dx_aff[j]) *
   1.748 +                 (z[j] + csa->alfa_aff_d * dz_aff[j]);
   1.749 +      csa->mu_aff = temp / (double)n;
   1.750 +      /* sigma = (mu_aff/mu)^3 */
   1.751 +      temp = csa->mu_aff / csa->mu;
   1.752 +      csa->sigma = temp * temp * temp;
   1.753 +      /* p = 0 */
   1.754 +      for (i = 1; i <= m; i++) p[i] = 0.0;
   1.755 +      /* q = 0 */
   1.756 +      for (j = 1; j <= n; j++) q[j] = 0.0;
   1.757 +      /* r = sigma * mu * e - X * Z * e */
   1.758 +      for (j = 1; j <= n; j++)
   1.759 +         r[j] = csa->sigma * csa->mu - dx_aff[j] * dz_aff[j];
   1.760 +      /* solve the second Newtonian system with the same coefficients
   1.761 +         but with altered right-hand sides */
   1.762 +      if (solve_NS(csa, p, q, r, dx_cc, dy_cc, dz_cc))
   1.763 +      {  ret = 1;
   1.764 +         goto done;
   1.765 +      }
   1.766 +      /* (dx,dy,dz) = (dx_aff,dy_aff,dz_aff) + (dx_cc,dy_cc,dz_cc) */
   1.767 +      for (j = 1; j <= n; j++) dx[j] = dx_aff[j] + dx_cc[j];
   1.768 +      for (i = 1; i <= m; i++) dy[i] = dy_aff[i] + dy_cc[i];
   1.769 +      for (j = 1; j <= n; j++) dz[j] = dz_aff[j] + dz_cc[j];
   1.770 +      /* alfa_max_p = inf{0 <= alfa <= 1 | x + alfa*dx >= 0} */
   1.771 +      /* alfa_max_d = inf{0 <= alfa <= 1 | z + alfa*dz >= 0} */
   1.772 +      csa->alfa_max_p = csa->alfa_max_d = 1.0;
   1.773 +      for (j = 1; j <= n; j++)
   1.774 +      {  if (dx[j] < 0.0)
   1.775 +         {  temp = - x[j] / dx[j];
   1.776 +            if (csa->alfa_max_p > temp) csa->alfa_max_p = temp;
   1.777 +         }
   1.778 +         if (dz[j] < 0.0)
   1.779 +         {  temp = - z[j] / dz[j];
   1.780 +            if (csa->alfa_max_d > temp) csa->alfa_max_d = temp;
   1.781 +         }
   1.782 +      }
   1.783 +      /* determine scale factors (not implemented yet) */
   1.784 +      gamma_p = 0.90;
   1.785 +      gamma_d = 0.90;
   1.786 +      /* compute the next point */
   1.787 +      for (j = 1; j <= n; j++)
   1.788 +      {  x[j] += gamma_p * csa->alfa_max_p * dx[j];
   1.789 +         xassert(x[j] > 0.0);
   1.790 +      }
   1.791 +      for (i = 1; i <= m; i++)
   1.792 +         y[i] += gamma_d * csa->alfa_max_d * dy[i];
   1.793 +      for (j = 1; j <= n; j++)
   1.794 +      {  z[j] += gamma_d * csa->alfa_max_d * dz[j];
   1.795 +         xassert(z[j] > 0.0);
   1.796 +      }
   1.797 +done: /* free working arrays */
   1.798 +      xfree(p);
   1.799 +      xfree(q);
   1.800 +      xfree(r);
   1.801 +      return ret;
   1.802 +}
   1.803 +
   1.804 +/***********************************************************************
   1.805 +*  terminate - deallocate common storage area
   1.806 +*
   1.807 +*  This routine frees all memory allocated to the common storage area
   1.808 +*  used by interior-point method routines. */
   1.809 +
   1.810 +static void terminate(struct csa *csa)
   1.811 +{     xfree(csa->D);
   1.812 +      xfree(csa->P);
   1.813 +      xfree(csa->S_ptr);
   1.814 +      xfree(csa->S_ind);
   1.815 +      xfree(csa->S_val);
   1.816 +      xfree(csa->S_diag);
   1.817 +      xfree(csa->U_ptr);
   1.818 +      xfree(csa->U_ind);
   1.819 +      xfree(csa->U_val);
   1.820 +      xfree(csa->U_diag);
   1.821 +      xfree(csa->phi_min);
   1.822 +      xfree(csa->best_x);
   1.823 +      xfree(csa->best_y);
   1.824 +      xfree(csa->best_z);
   1.825 +      xfree(csa->dx_aff);
   1.826 +      xfree(csa->dy_aff);
   1.827 +      xfree(csa->dz_aff);
   1.828 +      xfree(csa->dx_cc);
   1.829 +      xfree(csa->dy_cc);
   1.830 +      xfree(csa->dz_cc);
   1.831 +      return;
   1.832 +}
   1.833 +
   1.834 +/***********************************************************************
   1.835 +*  ipm_main - main interior-point method routine
   1.836 +*
   1.837 +*  This is a main routine of the primal-dual interior-point method.
   1.838 +*
   1.839 +*  The routine ipm_main returns one of the following codes:
   1.840 +*
   1.841 +*  0 - optimal solution found;
   1.842 +*  1 - problem has no feasible (primal or dual) solution;
   1.843 +*  2 - no convergence;
   1.844 +*  3 - iteration limit exceeded;
   1.845 +*  4 - numeric instability on solving Newtonian system.
   1.846 +*
   1.847 +*  In case of non-zero return code the routine returns the best point,
   1.848 +*  which has been reached during optimization. */
   1.849 +
   1.850 +static int ipm_main(struct csa *csa)
   1.851 +{     int m = csa->m;
   1.852 +      int n = csa->n;
   1.853 +      int i, j, status;
   1.854 +      double temp;
   1.855 +      /* choose initial point using Mehrotra's heuristic */
   1.856 +      if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.857 +         xprintf("Guessing initial point...\n");
   1.858 +      initial_point(csa);
   1.859 +      /* main loop starts here */
   1.860 +      if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.861 +         xprintf("Optimization begins...\n");
   1.862 +      for (;;)
   1.863 +      {  /* perform basic computations at the current point */
   1.864 +         basic_info(csa);
   1.865 +         /* save initial value of rmu */
   1.866 +         if (csa->iter == 0) csa->rmu0 = csa->rmu;
   1.867 +         /* accumulate values of min(phi[k]) and save the best point */
   1.868 +         xassert(csa->iter <= ITER_MAX);
   1.869 +         if (csa->iter == 0 || csa->phi_min[csa->iter-1] > csa->phi)
   1.870 +         {  csa->phi_min[csa->iter] = csa->phi;
   1.871 +            csa->best_iter = csa->iter;
   1.872 +            for (j = 1; j <= n; j++) csa->best_x[j] = csa->x[j];
   1.873 +            for (i = 1; i <= m; i++) csa->best_y[i] = csa->y[i];
   1.874 +            for (j = 1; j <= n; j++) csa->best_z[j] = csa->z[j];
   1.875 +            csa->best_obj = csa->obj;
   1.876 +         }
   1.877 +         else
   1.878 +            csa->phi_min[csa->iter] = csa->phi_min[csa->iter-1];
   1.879 +         /* display information at the current point */
   1.880 +         if (csa->parm->msg_lev >= GLP_MSG_ON)
   1.881 +            xprintf("%3d: obj = %17.9e; rpi = %8.1e; rdi = %8.1e; gap ="
   1.882 +               " %8.1e\n", csa->iter, csa->obj, csa->rpi, csa->rdi,
   1.883 +               csa->gap);
   1.884 +         /* check if the current point is optimal */
   1.885 +         if (csa->rpi < 1e-8 && csa->rdi < 1e-8 && csa->gap < 1e-8)
   1.886 +         {  if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.887 +               xprintf("OPTIMAL SOLUTION FOUND\n");
   1.888 +            status = 0;
   1.889 +            break;
   1.890 +         }
   1.891 +         /* check if the problem has no feasible solution */
   1.892 +         temp = 1e5 * csa->phi_min[csa->iter];
   1.893 +         if (temp < 1e-8) temp = 1e-8;
   1.894 +         if (csa->phi >= temp)
   1.895 +         {  if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.896 +               xprintf("PROBLEM HAS NO FEASIBLE PRIMAL/DUAL SOLUTION\n")
   1.897 +                  ;
   1.898 +            status = 1;
   1.899 +            break;
   1.900 +         }
   1.901 +         /* check for very slow convergence or divergence */
   1.902 +         if (((csa->rpi >= 1e-8 || csa->rdi >= 1e-8) && csa->rmu /
   1.903 +               csa->rmu0 >= 1e6) ||
   1.904 +               (csa->iter >= 30 && csa->phi_min[csa->iter] >= 0.5 *
   1.905 +               csa->phi_min[csa->iter - 30]))
   1.906 +         {  if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.907 +               xprintf("NO CONVERGENCE; SEARCH TERMINATED\n");
   1.908 +            status = 2;
   1.909 +            break;
   1.910 +         }
   1.911 +         /* check for maximal number of iterations */
   1.912 +         if (csa->iter == ITER_MAX)
   1.913 +         {  if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.914 +               xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n");
   1.915 +            status = 3;
   1.916 +            break;
   1.917 +         }
   1.918 +         /* start the next iteration */
   1.919 +         csa->iter++;
   1.920 +         /* factorize normal equation system */
   1.921 +         for (j = 1; j <= n; j++) csa->D[j] = csa->x[j] / csa->z[j];
   1.922 +         decomp_NE(csa);
   1.923 +         /* compute the next point using Mehrotra's predictor-corrector
   1.924 +            technique */
   1.925 +         if (make_step(csa))
   1.926 +         {  if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.927 +               xprintf("NUMERIC INSTABILITY; SEARCH TERMINATED\n");
   1.928 +            status = 4;
   1.929 +            break;
   1.930 +         }
   1.931 +      }
   1.932 +      /* restore the best point */
   1.933 +      if (status != 0)
   1.934 +      {  for (j = 1; j <= n; j++) csa->x[j] = csa->best_x[j];
   1.935 +         for (i = 1; i <= m; i++) csa->y[i] = csa->best_y[i];
   1.936 +         for (j = 1; j <= n; j++) csa->z[j] = csa->best_z[j];
   1.937 +         if (csa->parm->msg_lev >= GLP_MSG_ALL)
   1.938 +            xprintf("Best point %17.9e was reached on iteration %d\n",
   1.939 +               csa->best_obj, csa->best_iter);
   1.940 +      }
   1.941 +      /* return to the calling program */
   1.942 +      return status;
   1.943 +}
   1.944 +
   1.945 +/***********************************************************************
   1.946 +*  NAME
   1.947 +*
   1.948 +*  ipm_solve - core LP solver based on the interior-point method
   1.949 +*
   1.950 +*  SYNOPSIS
   1.951 +*
   1.952 +*  #include "glpipm.h"
   1.953 +*  int ipm_solve(glp_prob *P, const glp_iptcp *parm);
   1.954 +*
   1.955 +*  DESCRIPTION
   1.956 +*
   1.957 +*  The routine ipm_solve is a core LP solver based on the primal-dual
   1.958 +*  interior-point method.
   1.959 +*
   1.960 +*  The routine assumes the following standard formulation of LP problem
   1.961 +*  to be solved:
   1.962 +*
   1.963 +*     minimize
   1.964 +*
   1.965 +*        F = c[0] + c[1]*x[1] + c[2]*x[2] + ... + c[n]*x[n]
   1.966 +*
   1.967 +*     subject to linear constraints
   1.968 +*
   1.969 +*        a[1,1]*x[1] + a[1,2]*x[2] + ... + a[1,n]*x[n] = b[1]
   1.970 +*
   1.971 +*        a[2,1]*x[1] + a[2,2]*x[2] + ... + a[2,n]*x[n] = b[2]
   1.972 +*
   1.973 +*              . . . . . .
   1.974 +*
   1.975 +*        a[m,1]*x[1] + a[m,2]*x[2] + ... + a[m,n]*x[n] = b[m]
   1.976 +*
   1.977 +*     and non-negative variables
   1.978 +*
   1.979 +*        x[1] >= 0, x[2] >= 0, ..., x[n] >= 0
   1.980 +*
   1.981 +*  where:
   1.982 +*  F                    is the objective function;
   1.983 +*  x[1], ..., x[n]      are (structural) variables;
   1.984 +*  c[0]                 is a constant term of the objective function;
   1.985 +*  c[1], ..., c[n]      are objective coefficients;
   1.986 +*  a[1,1], ..., a[m,n]  are constraint coefficients;
   1.987 +*  b[1], ..., b[n]      are right-hand sides.
   1.988 +*
   1.989 +*  The solution is three vectors x, y, and z, which are stored by the
   1.990 +*  routine in the arrays x, y, and z, respectively. These vectors
   1.991 +*  correspond to the best primal-dual point found during optimization.
   1.992 +*  They are approximate solution of the following system (which is the
   1.993 +*  Karush-Kuhn-Tucker optimality conditions):
   1.994 +*
   1.995 +*     A*x      = b      (primal feasibility condition)
   1.996 +*
   1.997 +*     A'*y + z = c      (dual feasibility condition)
   1.998 +*
   1.999 +*     x'*z     = 0      (primal-dual complementarity condition)
  1.1000 +*
  1.1001 +*     x >= 0, z >= 0    (non-negativity condition)
  1.1002 +*
  1.1003 +*  where:
  1.1004 +*  x[1], ..., x[n]      are primal (structural) variables;
  1.1005 +*  y[1], ..., y[m]      are dual variables (Lagrange multipliers) for
  1.1006 +*                       equality constraints;
  1.1007 +*  z[1], ..., z[n]      are dual variables (Lagrange multipliers) for
  1.1008 +*                       non-negativity constraints.
  1.1009 +*
  1.1010 +*  RETURNS
  1.1011 +*
  1.1012 +*  0  LP has been successfully solved.
  1.1013 +*
  1.1014 +*  GLP_ENOCVG
  1.1015 +*     No convergence.
  1.1016 +*
  1.1017 +*  GLP_EITLIM
  1.1018 +*     Iteration limit exceeded.
  1.1019 +*
  1.1020 +*  GLP_EINSTAB
  1.1021 +*     Numeric instability on solving Newtonian system.
  1.1022 +*
  1.1023 +*  In case of non-zero return code the routine returns the best point,
  1.1024 +*  which has been reached during optimization. */
  1.1025 +
  1.1026 +int ipm_solve(glp_prob *P, const glp_iptcp *parm)
  1.1027 +{     struct csa _dsa, *csa = &_dsa;
  1.1028 +      int m = P->m;
  1.1029 +      int n = P->n;
  1.1030 +      int nnz = P->nnz;
  1.1031 +      GLPROW *row;
  1.1032 +      GLPCOL *col;
  1.1033 +      GLPAIJ *aij;
  1.1034 +      int i, j, loc, ret, *A_ind, *A_ptr;
  1.1035 +      double dir, *A_val, *b, *c, *x, *y, *z;
  1.1036 +      xassert(m > 0);
  1.1037 +      xassert(n > 0);
  1.1038 +      /* allocate working arrays */
  1.1039 +      A_ptr = xcalloc(1+m+1, sizeof(int));
  1.1040 +      A_ind = xcalloc(1+nnz, sizeof(int));
  1.1041 +      A_val = xcalloc(1+nnz, sizeof(double));
  1.1042 +      b = xcalloc(1+m, sizeof(double));
  1.1043 +      c = xcalloc(1+n, sizeof(double));
  1.1044 +      x = xcalloc(1+n, sizeof(double));
  1.1045 +      y = xcalloc(1+m, sizeof(double));
  1.1046 +      z = xcalloc(1+n, sizeof(double));
  1.1047 +      /* prepare rows and constraint coefficients */
  1.1048 +      loc = 1;
  1.1049 +      for (i = 1; i <= m; i++)
  1.1050 +      {  row = P->row[i];
  1.1051 +         xassert(row->type == GLP_FX);
  1.1052 +         b[i] = row->lb * row->rii;
  1.1053 +         A_ptr[i] = loc;
  1.1054 +         for (aij = row->ptr; aij != NULL; aij = aij->r_next)
  1.1055 +         {  A_ind[loc] = aij->col->j;
  1.1056 +            A_val[loc] = row->rii * aij->val * aij->col->sjj;
  1.1057 +            loc++;
  1.1058 +         }
  1.1059 +      }
  1.1060 +      A_ptr[m+1] = loc;
  1.1061 +      xassert(loc-1 == nnz);
  1.1062 +      /* prepare columns and objective coefficients */
  1.1063 +      if (P->dir == GLP_MIN)
  1.1064 +         dir = +1.0;
  1.1065 +      else if (P->dir == GLP_MAX)
  1.1066 +         dir = -1.0;
  1.1067 +      else
  1.1068 +         xassert(P != P);
  1.1069 +      c[0] = dir * P->c0;
  1.1070 +      for (j = 1; j <= n; j++)
  1.1071 +      {  col = P->col[j];
  1.1072 +         xassert(col->type == GLP_LO && col->lb == 0.0);
  1.1073 +         c[j] = dir * col->coef * col->sjj;
  1.1074 +      }
  1.1075 +      /* allocate and initialize the common storage area */
  1.1076 +      csa->m = m;
  1.1077 +      csa->n = n;
  1.1078 +      csa->A_ptr = A_ptr;
  1.1079 +      csa->A_ind = A_ind;
  1.1080 +      csa->A_val = A_val;
  1.1081 +      csa->b = b;
  1.1082 +      csa->c = c;
  1.1083 +      csa->x = x;
  1.1084 +      csa->y = y;
  1.1085 +      csa->z = z;
  1.1086 +      csa->parm = parm;
  1.1087 +      initialize(csa);
  1.1088 +      /* solve LP with the interior-point method */
  1.1089 +      ret = ipm_main(csa);
  1.1090 +      /* deallocate the common storage area */
  1.1091 +      terminate(csa);
  1.1092 +      /* determine solution status */
  1.1093 +      if (ret == 0)
  1.1094 +      {  /* optimal solution found */
  1.1095 +         P->ipt_stat = GLP_OPT;
  1.1096 +         ret = 0;
  1.1097 +      }
  1.1098 +      else if (ret == 1)
  1.1099 +      {  /* problem has no feasible (primal or dual) solution */
  1.1100 +         P->ipt_stat = GLP_NOFEAS;
  1.1101 +         ret = 0;
  1.1102 +      }
  1.1103 +      else if (ret == 2)
  1.1104 +      {  /* no convergence */
  1.1105 +         P->ipt_stat = GLP_INFEAS;
  1.1106 +         ret = GLP_ENOCVG;
  1.1107 +      }
  1.1108 +      else if (ret == 3)
  1.1109 +      {  /* iteration limit exceeded */
  1.1110 +         P->ipt_stat = GLP_INFEAS;
  1.1111 +         ret = GLP_EITLIM;
  1.1112 +      }
  1.1113 +      else if (ret == 4)
  1.1114 +      {  /* numeric instability on solving Newtonian system */
  1.1115 +         P->ipt_stat = GLP_INFEAS;
  1.1116 +         ret = GLP_EINSTAB;
  1.1117 +      }
  1.1118 +      else
  1.1119 +         xassert(ret != ret);
  1.1120 +      /* store row solution components */
  1.1121 +      for (i = 1; i <= m; i++)
  1.1122 +      {  row = P->row[i];
  1.1123 +         row->pval = row->lb;
  1.1124 +         row->dval = dir * y[i] * row->rii;
  1.1125 +      }
  1.1126 +      /* store column solution components */
  1.1127 +      P->ipt_obj = P->c0;
  1.1128 +      for (j = 1; j <= n; j++)
  1.1129 +      {  col = P->col[j];
  1.1130 +         col->pval = x[j] * col->sjj;
  1.1131 +         col->dval = dir * z[j] / col->sjj;
  1.1132 +         P->ipt_obj += col->coef * col->pval;
  1.1133 +      }
  1.1134 +      /* free working arrays */
  1.1135 +      xfree(A_ptr);
  1.1136 +      xfree(A_ind);
  1.1137 +      xfree(A_val);
  1.1138 +      xfree(b);
  1.1139 +      xfree(c);
  1.1140 +      xfree(x);
  1.1141 +      xfree(y);
  1.1142 +      xfree(z);
  1.1143 +      return ret;
  1.1144 +}
  1.1145 +
  1.1146 +/* eof */