glpipm.c

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

      cx = 0.0;
      for (j = 1; j <= n; j++) cx += c[j] * x[j];
      /* compute primal-dual gap */
      dsa->gap = fabs(cx - by) / (1.0 + fabs(cx));
      /* compute merit function */
      dsa->phi = 0.0;
      dsa->phi += norm1 / (bnorm > 1.0 ? bnorm : 1.0);
      dsa->phi += norm2 / (cnorm > 1.0 ? cnorm : 1.0);
      temp = 1.0;
      if (temp < bnorm) temp = bnorm;
      if (temp < cnorm) temp = cnorm;
      dsa->phi += fabs(cx - by) / temp;
      /* compute duality measure */
      temp = 0.0;
      for (j = 1; j <= n; j++) temp += x[j] * z[j];
      dsa->mu = temp / (double)n;
      /* compute the ratio of infeasibility to mu */
      dsa->rmu = (norm1 > norm2 ? norm1 : norm2) / dsa->mu;
      return;
}

/*----------------------------------------------------------------------
-- make_step - compute next point using Mehrotra's technique.
--
-- This routine computes the next point using the predictor-corrector
-- technique 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.
--
-- At first the routine computes so called affine scaling (predictor)
-- direction (dx_aff,dy_aff,dz_aff) which is a solution of the system:
--
--    A*dx_aff                       = b - A*x
--              A'*dy_aff +   dz_aff = c - A'*y - z
--    Z*dx_aff            + X*dz_aff = - X*Z*e
--
-- where (x,y,z) is the current point, X = diag(x[j]), Z = diag(z[j]),
-- e = (1,...,1)'.
--
-- Then the routine computes the centering parameter sigma, using the
-- following Mehrotra's heuristic:
--
--    alfa_aff_p = inf{0 <= alfa <= 1 | x+alfa*dx_aff >= 0}
--
--    alfa_aff_d = inf{0 <= alfa <= 1 | z+alfa*dz_aff >= 0}
--
--    mu_aff = (x+alfa_aff_p*dx_aff)'*(z+alfa_aff_d*dz_aff)/n
--
--    sigma = (mu_aff/mu)^3
--
-- where alfa_aff_p is the maximal stepsize along the affine scaling
-- direction in the primal space, alfa_aff_d is the maximal stepsize
-- along the same direction in the dual space.
--
-- After determining sigma the routine computes so called centering
-- (corrector) direction (dx_cc,dy_cc,dz_cc) which is the solution of
-- the system:
--
--    A*dx_cc                     = 0
--             A'*dy_cc +   dz_cc = 0
--    Z*dx_cc           + X*dz_cc = sigma*mu*e - X*Z*e
--
-- Finally, the routine computes the combined direction
--
--    (dx,dy,dz) = (dx_aff,dy_aff,dz_aff) + (dx_cc,dy_cc,dz_cc)
--
-- and determines maximal primal and dual stepsizes along the combined
-- direction:
--
--    alfa_max_p = inf{0 <= alfa <= 1 | x+alfa*dx >= 0}
--
--    alfa_max_d = inf{0 <= alfa <= 1 | z+alfa*dz >= 0}
--
-- In order to prevent the next point to be too close to the boundary
-- of the positive ortant, the routine decreases maximal stepsizes:
--
--    alfa_p = gamma_p * alfa_max_p
--
--    alfa_d = gamma_d * alfa_max_d
--
-- where gamma_p and gamma_d are scaling factors, and computes the next
-- point:
--
--    x_new = x + alfa_p * dx
--
--    y_new = y + alfa_d * dy
--
--    z_new = z + alfa_d * dz
--
-- which becomes the current point on the next iteration. */

static int make_step(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 *dx_aff = dsa->dx_aff;
      double *dy_aff = dsa->dy_aff;
      double *dz_aff = dsa->dz_aff;
      double *dx_cc = dsa->dx_cc;
      double *dy_cc = dsa->dy_cc;
      double *dz_cc = dsa->dz_cc;
      double *dx = dsa->dx;
      double *dy = dsa->dy;
      double *dz = dsa->dz;
      int i, j, ret = 0;
      double temp, gamma_p, gamma_d, *p, *q, *r;
      /* allocate working arrays */
      p = xcalloc(1+m, sizeof(double));
      q = xcalloc(1+n, sizeof(double));
      r = xcalloc(1+n, sizeof(double));
      /* p = b - A*x */
      A_by_vec(dsa, x, p);
      for (i = 1; i <= m; i++) p[i] = b[i] - p[i];
      /* q = c - A'*y - z */
      AT_by_vec(dsa, y,q);
      for (j = 1; j <= n; j++) q[j] = c[j] - q[j] - z[j];
      /* r = - X * Z * e */
      for (j = 1; j <= n; j++) r[j] = - x[j] * z[j];
      /* solve the first Newtonian system */
      if (solve_NS(dsa, p, q, r, dx_aff, dy_aff, dz_aff))
      {  ret = 1;
         goto done;
      }
      /* alfa_aff_p = inf{0 <= alfa <= 1 | x + alfa*dx_aff >= 0} */
      /* alfa_aff_d = inf{0 <= alfa <= 1 | z + alfa*dz_aff >= 0} */
      dsa->alfa_aff_p = dsa->alfa_aff_d = 1.0;
      for (j = 1; j <= n; j++)
      {  if (dx_aff[j] < 0.0)
         {  temp = - x[j] / dx_aff[j];
            if (dsa->alfa_aff_p > temp) dsa->alfa_aff_p = temp;
         }
         if (dz_aff[j] < 0.0)
         {  temp = - z[j] / dz_aff[j];
            if (dsa->alfa_aff_d > temp) dsa->alfa_aff_d = temp;
         }
      }
      /* mu_aff = (x+alfa_aff_p*dx_aff)' * (z+alfa_aff_d*dz_aff) / n */
      temp = 0.0;
      for (j = 1; j <= n; j++)
         temp += (x[j] + dsa->alfa_aff_p * dx_aff[j]) *
                 (z[j] + dsa->alfa_aff_d * dz_aff[j]);
      dsa->mu_aff = temp / (double)n;
      /* sigma = (mu_aff/mu)^3 */
      temp = dsa->mu_aff / dsa->mu;
      dsa->sigma = temp * temp * temp;
      /* p = 0 */
      for (i = 1; i <= m; i++) p[i] = 0.0;
      /* q = 0 */
      for (j = 1; j <= n; j++) q[j] = 0.0;
      /* r = sigma * mu * e - X * Z * e */
      for (j = 1; j <= n; j++)
         r[j] = dsa->sigma * dsa->mu - dx_aff[j] * dz_aff[j];
      /* solve the second Newtonian system with the same coefficients
         but with altered right-hand sides */
      if (solve_NS(dsa, p, q, r, dx_cc, dy_cc, dz_cc))
      {  ret = 1;
         goto done;
      }
      /* (dx,dy,dz) = (dx_aff,dy_aff,dz_aff) + (dx_cc,dy_cc,dz_cc) */
      for (j = 1; j <= n; j++) dx[j] = dx_aff[j] + dx_cc[j];
      for (i = 1; i <= m; i++) dy[i] = dy_aff[i] + dy_cc[i];
      for (j = 1; j <= n; j++) dz[j] = dz_aff[j] + dz_cc[j];
      /* alfa_max_p = inf{0 <= alfa <= 1 | x + alfa*dx >= 0} */
      /* alfa_max_d = inf{0 <= alfa <= 1 | z + alfa*dz >= 0} */
      dsa->alfa_max_p = dsa->alfa_max_d = 1.0;
      for (j = 1; j <= n; j++)
      {  if (dx[j] < 0.0)
         {  temp = - x[j] / dx[j];
            if (dsa->alfa_max_p > temp) dsa->alfa_max_p = temp;
         }
         if (dz[j] < 0.0)
         {  temp = - z[j] / dz[j];
            if (dsa->alfa_max_d > temp) dsa->alfa_max_d = temp;
         }
      }
      /* determine scale factors (not implemented yet) */
      gamma_p = 0.90;
      gamma_d = 0.90;
      /* compute the next point */
      for (j = 1; j <= n; j++)
      {  x[j] += gamma_p * dsa->alfa_max_p * dx[j];
         xassert(x[j] > 0.0);
      }
      for (i = 1; i <= m; i++)
         y[i] += gamma_d * dsa->alfa_max_d * dy[i];
      for (j = 1; j <= n; j++)
      {  z[j] += gamma_d * dsa->alfa_max_d * dz[j];
         xassert(z[j] > 0.0);
      }
done: /* free working arrays */
      xfree(p);
      xfree(q);
      xfree(r);
      return ret;
}

/*----------------------------------------------------------------------
-- terminate - deallocate working area.
--
-- This routine frees all memory allocated to the working area used by
-- all interior-point method routines. */

static void terminate(struct dsa *dsa)
{     xfree(dsa->D);
      xfree(dsa->P);
      xfree(dsa->S_ptr);
      xfree(dsa->S_ind);
      xfree(dsa->S_val);
      xfree(dsa->S_diag);
      xfree(dsa->U_ptr);
      xfree(dsa->U_ind);
      xfree(dsa->U_val);
      xfree(dsa->U_diag);
      xfree(dsa->phi_min);
      xfree(dsa->best_x);
      xfree(dsa->best_y);
      xfree(dsa->best_z);
      xfree(dsa->dx_aff);
      xfree(dsa->dy_aff);
      xfree(dsa->dz_aff);
      xfree(dsa->dx_cc);
      xfree(dsa->dy_cc);
      xfree(dsa->dz_cc);
      return;
}

/*----------------------------------------------------------------------
-- ipm_main - main interior-point method routine.
--
-- This is a main routine of the primal-dual interior-point method. */

int ipm_main(int m, int n, int A_ptr[], int A_ind[], double A_val[],
      double b[], double c[], double x[], double y[], double z[])
{     struct dsa _dsa, *dsa = &_dsa;
      int i, j, status;
      double temp;
      /* allocate and initialize working area */
      xassert(m > 0);
      xassert(n > 0);
      initialize(dsa, m, n, A_ptr, A_ind, A_val, b, c, x, y, z);
      /* choose initial point using Mehrotra's heuristic */
      xprintf("lpx_interior: guessing initial point...\n");
      initial_point(dsa);
      /* main loop starts here */
      xprintf("Optimization begins...\n");
      for (;;)
      {  /* perform basic computations at the current point */
         basic_info(dsa);
         /* save initial value of rmu */
         if (dsa->iter == 0) dsa->rmu0 = dsa->rmu;
         /* accumulate values of min(phi[k]) and save the best point */
         xassert(dsa->iter <= ITER_MAX);
         if (dsa->iter == 0 || dsa->phi_min[dsa->iter-1] > dsa->phi)
         {  dsa->phi_min[dsa->iter] = dsa->phi;
            dsa->best_iter = dsa->iter;
            for (j = 1; j <= n; j++) dsa->best_x[j] = dsa->x[j];
            for (i = 1; i <= m; i++) dsa->best_y[i] = dsa->y[i];
            for (j = 1; j <= n; j++) dsa->best_z[j] = dsa->z[j];
            dsa->best_obj = dsa->obj;
         }
         else
            dsa->phi_min[dsa->iter] = dsa->phi_min[dsa->iter-1];
         /* display information at the current point */
         xprintf("%3d: obj = %17.9e; rpi = %8.1e; rdi = %8.1e; gap = %8"
            ".1e\n", dsa->iter, dsa->obj, dsa->rpi, dsa->rdi, dsa->gap);
         /* check if the current point is optimal */
         if (dsa->rpi < 1e-8 && dsa->rdi < 1e-8 && dsa->gap < 1e-8)
         {  xprintf("OPTIMAL SOLUTION FOUND\n");
            status = 0;
            break;
         }
         /* check if the problem has no feasible solutions */
         temp = 1e5 * dsa->phi_min[dsa->iter];
         if (temp < 1e-8) temp = 1e-8;
         if (dsa->phi >= temp)
         {  xprintf("PROBLEM HAS NO FEASIBLE PRIMAL/DUAL SOLUTION\n");
            status = 1;
            break;
         }
         /* check for very slow convergence or divergence */
         if (((dsa->rpi >= 1e-8 || dsa->rdi >= 1e-8) && dsa->rmu /
               dsa->rmu0 >= 1e6) ||
               (dsa->iter >= 30 && dsa->phi_min[dsa->iter] >= 0.5 *
               dsa->phi_min[dsa->iter - 30]))
         {  xprintf("NO CONVERGENCE; SEARCH TERMINATED\n");
            status = 2;
            break;
         }
         /* check for maximal number of iterations */
         if (dsa->iter == ITER_MAX)
         {  xprintf("ITERATION LIMIT EXCEEDED; SEARCH TERMINATED\n");
            status = 3;
            break;
         }
         /* start the next iteration */
         dsa->iter++;
         /* factorize normal equation system */
         for (j = 1; j <= n; j++) dsa->D[j] = dsa->x[j] / dsa->z[j];
         decomp_NE(dsa);
         /* compute the next point using Mehrotra's predictor-corrector
            technique */
         if (make_step(dsa))
         {  xprintf("NUMERIC INSTABILITY; SEARCH TERMINATED\n");
            status = 4;
            break;
         }
      }
      /* restore the best point */
      if (status != 0)
      {  for (j = 1; j <= n; j++) dsa->x[j] = dsa->best_x[j];
         for (i = 1; i <= m; i++) dsa->y[i] = dsa->best_y[i];
         for (j = 1; j <= n; j++) dsa->z[j] = dsa->best_z[j];
         xprintf("The best point %17.9e was reached on iteration %d\n",
            dsa->best_obj, dsa->best_iter);
      }
      /* deallocate working area */
      terminate(dsa);
      /* return to the calling program */
      return status;
}

/* eof */