glpipm.c

著名的大规模线性规划求解器源码GLPK.C语言版本,可以修剪.内有详细帮助文档.

-- a matrix transposed to the constraint matrix A. */

static void AT_by_vec(struct dsa *dsa, double x[], double y[])
{     /* compute y = A'*x, where A' is transposed to A */
      int m = dsa->m;
      int n = dsa->n;
      int *A_ptr = dsa->A_ptr;
      int *A_ind = dsa->A_ind;
      double *A_val = dsa->A_val;
      int i, j, t, beg, end;
      double temp;
      for (j = 1; j <= n; j++) y[j] = 0.0;
      for (i = 1; i <= m; i++)
      {  temp = x[i];
         if (temp == 0.0) continue;
         beg = A_ptr[i], end = A_ptr[i+1];
         for (t = beg; t < end; t++) y[A_ind[t]] += A_val[t] * temp;
      }
      return;
}

/*----------------------------------------------------------------------
-- decomp_NE - numeric factorization of matrix S = P*A*D*A'*P'.
--
-- This routine implements numeric phase of Cholesky factorization of
-- the matrix S = P*A*D*A'*P', which is a permuted matrix of the normal
-- equation system. Matrix D is assumed to be already computed */

static void decomp_NE(struct dsa *dsa)
{     adat_numeric(dsa->m, dsa->n, dsa->P, dsa->A_ptr, dsa->A_ind,
         dsa->A_val, dsa->D, dsa->S_ptr, dsa->S_ind, dsa->S_val,
         dsa->S_diag);
      chol_numeric(dsa->m, dsa->S_ptr, dsa->S_ind, dsa->S_val,
         dsa->S_diag, dsa->U_ptr, dsa->U_ind, dsa->U_val, dsa->U_diag);
      return;
}

/*----------------------------------------------------------------------
-- solve_NE - solve normal equation system.
--
-- This routine solves the normal equation system
--
--    A*D*A'*y = h.
--
-- It is assumed that the matrix A*D*A' has been previously factorized
-- by the routine decomp_NE.
--
-- On entry the array y contains the vector of right-hand sides h. On
-- exit this array contains the computed vector of unknowns y.
--
-- Once the vector y has been computed the routine checks for numeric
-- stability. If the residual vector
--
--    r = A*D*A'*y - h
--
-- is relatively small, the routine returns zero, otherwise non-zero is
-- returned. */

static int solve_NE(struct dsa *dsa, double y[])
{     int m = dsa->m;
      int n = dsa->n;
      int *P = dsa->P;
      int i, j, ret = 0;
      double *h, *r, *w;
      /* save vector of right-hand sides h */
      h = xcalloc(1+m, sizeof(double));
      for (i = 1; i <= m; i++) h[i] = y[i];
      /* solve normal equation system (A*D*A')*y = h */
      /* since S = P*A*D*A'*P' = U'*U, then A*D*A' = P'*U'*U*P, so we
         have inv(A*D*A') = P'*inv(U)*inv(U')*P */
      /* w := P*h */
      w = xcalloc(1+m, sizeof(double));
      for (i = 1; i <= m; i++) w[i] = y[P[i]];
      /* w := inv(U')*w */
      ut_solve(m, dsa->U_ptr, dsa->U_ind, dsa->U_val, dsa->U_diag, w);
      /* w := inv(U)*w */
      u_solve(m, dsa->U_ptr, dsa->U_ind, dsa->U_val, dsa->U_diag, w);
      /* y := P'*w */
      for (i = 1; i <= m; i++) y[i] = w[P[m+i]];
      xfree(w);
      /* compute residual vector r = A*D*A'*y - h */
      r = xcalloc(1+m, sizeof(double));
      /* w := A'*y */
      w = xcalloc(1+n, sizeof(double));
      AT_by_vec(dsa, y, w);
      /* w := D*w */
      for (j = 1; j <= n; j++) w[j] *= dsa->D[j];
      /* r := A*w */
      A_by_vec(dsa, w, r);
      xfree(w);
      /* r := r - h */
      for (i = 1; i <= m; i++) r[i] -= h[i];
      /* check for numeric stability */
      for (i = 1; i <= m; i++)
      {  if (fabs(r[i]) / (1.0 + fabs(h[i])) > 1e-4)
         {  ret = 1;
            break;
         }
      }
      xfree(h);
      xfree(r);
      return ret;
}

/*----------------------------------------------------------------------
-- solve_NS - solve Newtonian system.
--
-- This routine solves the Newtonian system:
--
--    A*dx               = p
--          A'*dy +   dz = q
--    Z*dx        + X*dz = r
--
-- where X = diag(x[j]), Z = diag(z[j]), by reducing it to the normal
-- equation system:
--
--    (A*inv(Z)*X*A')*dy = A*inv(Z)*(X*q-r)+p
--
-- (it is assumed that the matrix A*inv(Z)*X*A' has been factorized by
-- means of the decomp_NE routine).
--
-- Once vector dy has been computed the routine computes vectors dx and
-- dz as follows:
--
--    dx = inv(Z)*(X*(A'*dy-q)+r)
--
--    dz = inv(X)*(r-Z*dx)
--
-- The routine solve_NS returns a code reported by the routine solve_NE
-- which solves the normal equation system. */

static int solve_NS(struct dsa *dsa, double p[], double q[], double r[],
      double dx[], double dy[], double dz[])
{     int m = dsa->m;
      int n = dsa->n;
      double *x = dsa->x;
      double *z = dsa->z;
      int i, j, ret;
      double *w = dx;
      /* compute the vector of right-hand sides A*inv(Z)*(X*q-r)+p for
         the normal equation system */
      for (j = 1; j <= n; j++)
         w[j] = (x[j] * q[j] - r[j]) / z[j];
      A_by_vec(dsa, w, dy);
      for (i = 1; i <= m; i++) dy[i] += p[i];
      /* solve the normal equation system to compute vector dy */
      ret = solve_NE(dsa, dy);
      /* compute vectors dx and dz */
      AT_by_vec(dsa, dy, dx);
      for (j = 1; j <= n; j++)
      {  dx[j] = (x[j] * (dx[j] - q[j]) + r[j]) / z[j];
         dz[j] = (r[j] - z[j] * dx[j]) / x[j];
      }
      return ret;
}

/*----------------------------------------------------------------------
-- initial_point - choose initial point using Mehrotra's heuristic.
--
-- This routine chooses a starting point using a heuristic proposed in
-- the paper:
--
-- S. Mehrotra. On the implementation of a primal-dual interior point
-- method. SIAM J. on Optim., 2(4), pp. 575-601, 1992.
--
-- The starting point x in the primal space is chosen as a solution of
-- the following least squares problem:
--
--    minimize    ||x||
--
--    subject to  A*x = b
--
-- which can be computed explicitly as follows:
--
--    x = A'*inv(A*A')*b
--
-- Similarly, the starting point (y, z) in the dual space is chosen as
-- a solution of the following least squares problem:
--
--    minimize    ||z||
--
--    subject to  A'*y + z = c
--
-- which can be computed explicitly as follows:
--
--    y = inv(A*A')*A*c
--
--    z = c - A'*y
--
-- However, some components of the vectors x and z may be non-positive
-- or close to zero, so the routine uses a Mehrotra's heuristic to find
-- a more appropriate starting point. */

static void initial_point(struct dsa *dsa)
{     int m = dsa->m;
      int n = dsa->n;
      double *b = dsa->b;
      double *c = dsa->c;
      double *x = dsa->x;
      double *y = dsa->y;
      double *z = dsa->z;
      double *D = dsa->D;
      int i, j;
      double dp, dd, ex, ez, xz;
      /* factorize A*A' */
      for (j = 1; j <= n; j++) D[j] = 1.0;
      decomp_NE(dsa);
      /* x~ = A'*inv(A*A')*b */
      for (i = 1; i <= m; i++) y[i] = b[i];
      solve_NE(dsa, y);
      AT_by_vec(dsa, y, x);
      /* y~ = inv(A*A')*A*c */
      A_by_vec(dsa, c, y);
      solve_NE(dsa, y);
      /* z~ = c - A'*y~ */
      AT_by_vec(dsa, y,z);
      for (j = 1; j <= n; j++) z[j] = c[j] - z[j];
      /* use Mehrotra's heuristic in order to choose more appropriate
         starting point with positive components of vectors x and z */
      dp = dd = 0.0;
      for (j = 1; j <= n; j++)
      {  if (dp < -1.5 * x[j]) dp = -1.5 * x[j];
         if (dd < -1.5 * z[j]) dd = -1.5 * z[j];
      }
      /* note that b = 0 involves x = 0, and c = 0 involves y = 0 and
         z = 0, so we need to be careful */
      if (dp == 0.0) dp = 1.5;
      if (dd == 0.0) dd = 1.5;
      ex = ez = xz = 0.0;
      for (j = 1; j <= n; j++)
      {  ex += (x[j] + dp);
         ez += (z[j] + dd);
         xz += (x[j] + dp) * (z[j] + dd);
      }
      dp += 0.5 * (xz / ez);
      dd += 0.5 * (xz / ex);
      for (j = 1; j <= n; j++)
      {  x[j] += dp;
         z[j] += dd;
         xassert(x[j] > 0.0 && z[j] > 0.0);
      }
      return;
}

/*----------------------------------------------------------------------
-- basic_info - perform basic computations at the current point.
--
-- This routine computes the following quantities at the current point:
--
-- value of the objective function:
--
--    F = c'*x + c[0]
--
-- relative primal infeasibility:
--
--    rpi = ||A*x-b|| / (1+||b||)
--
-- relative dual infeasibility:
--
--    rdi = ||A'*y+z-c|| / (1+||c||)
--
-- primal-dual gap (a relative difference between the primal and the
-- dual objective function values):
--
--    gap = |c'*x-b'*y| / (1+|c'*x|)
--
-- merit function:
--
--    phi = ||A*x-b|| / max(1,||b||) + ||A'*y+z-c|| / max(1,||c||) +
--
--        + |c'*x-b'*y| / max(1,||b||,||c||)
--
-- duality measure:
--
--    mu = x'*z / n
--
-- the ratio of infeasibility to mu:
--
--    rmu = max(||A*x-b||,||A'*y+z-c||) / mu
--
-- where ||*|| denotes euclidian norm, *' denotes transposition */

static void basic_info(struct dsa *dsa)
{     int m = dsa->m;
      int n = dsa->n;
      double *b = dsa->b;
      double *c = dsa->c;
      double *x = dsa->x;
      double *y = dsa->y;
      double *z = dsa->z;
      int i, j;
      double norm1, bnorm, norm2, cnorm, cx, by, *work, temp;
      /* compute value of the objective function */
      temp = c[0];
      for (j = 1; j <= n; j++) temp += c[j] * x[j];
      dsa->obj = temp;
      /* norm1 = ||A*x-b|| */
      work = xcalloc(1+m, sizeof(double));
      A_by_vec(dsa, x, work);
      norm1 = 0.0;
      for (i = 1; i <= m; i++)
         norm1 += (work[i] - b[i]) * (work[i] - b[i]);
      norm1 = sqrt(norm1);
      xfree(work);
      /* bnorm = ||b|| */
      bnorm = 0.0;
      for (i = 1; i <= m; i++) bnorm += b[i] * b[i];
      bnorm = sqrt(bnorm);
      /* compute relative primal infeasibility */
      dsa->rpi = norm1 / (1.0 + bnorm);
      /* norm2 = ||A'*y+z-c|| */
      work = xcalloc(1+n, sizeof(double));
      AT_by_vec(dsa, y, work);
      norm2 = 0.0;
      for (j = 1; j <= n; j++)
         norm2 += (work[j] + z[j] - c[j]) * (work[j] + z[j] - c[j]);
      norm2 = sqrt(norm2);
      xfree(work);
      /* cnorm = ||c|| */
      cnorm = 0.0;
      for (j = 1; j <= n; j++) cnorm += c[j] * c[j];
      cnorm = sqrt(cnorm);
      /* compute relative dual infeasibility */
      dsa->rdi = norm2 / (1.0 + cnorm);
      /* by = b'*y */
      by = 0.0;
      for (i = 1; i <= m; i++) by += b[i] * y[i];
      /* cx = c'*x */