glpios06.c

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

{     int j;
      double *aa, bb;
      aa = alpha, bb = b;
      for (j = 1; j <= n; j++)
      {  aa[j] = a[j] / delta;
         if (cset[j])
            aa[j] = - aa[j], bb -= a[j] * u[j];
      }
      bb /= delta;
      if (mir_ineq(n, aa, bb, alpha, beta, gamma)) return 1;
      for (j = 1; j <= n; j++)
      {  if (cset[j])
            alpha[j] = - alpha[j], *beta += alpha[j] * u[j];
      }
      *gamma /= delta;
      return 0;
}

/***********************************************************************
*  cmir_sep - c-MIR separation heuristic
*
*  Given the mixed knapsack set
*
*      MK              |N|
*     X   = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s,
*                      +      +  j in N
*
*             x[j] <= u[j]}
*
*                           *   *
*  and a fractional point (x , s ), this routine tries to construct
*  c-MIR inequality
*
*     sum   alpha[j] * x[j] <= beta + gamma * s,
*    j in N
*                      MK
*  which is valid for X   and has (desirably maximal) violation at the
*  fractional point given. This is attained by choosing an appropriate
*  set C of variables to be complemented and a divisor delta > 0, which
*  together define corresponding c-MIR inequality.
*
*  If a violated c-MIR inequality has been successfully constructed,
*  the routine returns its violation:
*
*                       *                      *
*     sum   alpha[j] * x [j] - beta - gamma * s ,
*    j in N
*
*  which is positive. In case of failure the routine returns zero. */

struct vset { int j; double v; };

static int cmir_cmp(const void *p1, const void *p2)
{     const struct vset *v1 = p1, *v2 = p2;
      if (v1->v < v2->v) return -1;
      if (v1->v > v2->v) return +1;
      return 0;
}

static double cmir_sep(const int n, const double a[], const double b,
      const double u[], const double x[], const double s,
      double alpha[], double *beta, double *gamma)
{     int fail, j, k, nv, v;
      double delta, eps, d_try[1+3], r, r_best;
      char *cset;
      struct vset *vset;
      /* allocate working arrays */
      cset = xcalloc(1+n, sizeof(char));
      vset = xcalloc(1+n, sizeof(struct vset));
      /* choose initial C */
      for (j = 1; j <= n; j++)
         cset[j] = (char)(x[j] >= 0.5 * u[j]);
      /* choose initial delta */
      r_best = delta = 0.0;
      for (j = 1; j <= n; j++)
      {  xassert(a[j] != 0.0);
         /* if x[j] is close to its bounds, skip it */
         eps = 1e-9 * (1.0 + fabs(u[j]));
         if (x[j] < eps || x[j] > u[j] - eps) continue;
         /* try delta = |a[j]| to construct c-MIR inequality */
         fail = cmir_ineq(n, a, b, u, cset, fabs(a[j]), alpha, beta,
            gamma);
         if (fail) continue;
         /* compute violation */
         r = - (*beta) - (*gamma) * s;
         for (k = 1; k <= n; k++) r += alpha[k] * x[k];
         if (r_best < r) r_best = r, delta = fabs(a[j]);
      }
      if (r_best < 0.001) r_best = 0.0;
      if (r_best == 0.0) goto done;
      xassert(delta > 0.0);
      /* try to increase violation by dividing delta by 2, 4, and 8,
         respectively */
      d_try[1] = delta / 2.0;
      d_try[2] = delta / 4.0;
      d_try[3] = delta / 8.0;
      for (j = 1; j <= 3; j++)
      {  /* construct c-MIR inequality */
         fail = cmir_ineq(n, a, b, u, cset, d_try[j], alpha, beta,
            gamma);
         if (fail) continue;
         /* compute violation */
         r = - (*beta) - (*gamma) * s;
         for (k = 1; k <= n; k++) r += alpha[k] * x[k];
         if (r_best < r) r_best = r, delta = d_try[j];
      }
      /* build subset of variables lying strictly between their bounds
         and order it by nondecreasing values of |x[j] - u[j]/2| */
      nv = 0;
      for (j = 1; j <= n; j++)
      {  /* if x[j] is close to its bounds, skip it */
         eps = 1e-9 * (1.0 + fabs(u[j]));
         if (x[j] < eps || x[j] > u[j] - eps) continue;
         /* add x[j] to the subset */
         nv++;
         vset[nv].j = j;
         vset[nv].v = fabs(x[j] - 0.5 * u[j]);
      }
      qsort(&vset[1], nv, sizeof(struct vset), cmir_cmp);
      /* try to increase violation by successively complementing each
         variable in the subset */
      for (v = 1; v <= nv; v++)
      {  j = vset[v].j;
         /* replace x[j] by its complement or vice versa */
         cset[j] = (char)!cset[j];
         /* construct c-MIR inequality */
         fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
         /* restore the variable */
         cset[j] = (char)!cset[j];
         /* do not replace the variable in case of failure */
         if (fail) continue;
         /* compute violation */
         r = - (*beta) - (*gamma) * s;
         for (k = 1; k <= n; k++) r += alpha[k] * x[k];
         if (r_best < r) r_best = r, cset[j] = (char)!cset[j];
      }
      /* construct the best c-MIR inequality chosen */
      fail = cmir_ineq(n, a, b, u, cset, delta, alpha, beta, gamma);
      xassert(!fail);
done: /* free working arrays */
      xfree(cset);
      xfree(vset);
      /* return to the calling routine */
      return r_best;
}

static double generate(struct MIR *mir)
{     /* try to generate violated c-MIR cut for modified constraint */
      int m = mir->m;
      int n = mir->n;
      int j, k, kk, nint;
      double s, *u, *x, *alpha, r_best = 0.0, b, beta, gamma;
      ios_copy_vec(mir->cut_vec, mir->mod_vec);
      mir->cut_rhs = mir->mod_rhs;
      /* remove small terms, which can appear due to substitution of
         variable bounds */
      ios_clean_vec(mir->cut_vec, DBL_EPSILON);
#if _MIR_DEBUG
      ios_check_vec(mir->cut_vec);
#endif
      /* remove positive continuous terms to obtain MK relaxation */
      for (j = 1; j <= mir->cut_vec->nnz; j++)
      {  k = mir->cut_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (!mir->isint[k] && mir->cut_vec->val[j] > 0.0)
            mir->cut_vec->val[j] = 0.0;
      }
      ios_clean_vec(mir->cut_vec, 0.0);
#if _MIR_DEBUG
      ios_check_vec(mir->cut_vec);
#endif
      /* move integer terms to the beginning of the sparse vector and
         determine the number of integer variables */
      nint = 0;
      for (j = 1; j <= mir->cut_vec->nnz; j++)
      {  k = mir->cut_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (mir->isint[k])
         {  double temp;
            nint++;
            /* interchange elements [nint] and [j] */
            kk = mir->cut_vec->ind[nint];
            mir->cut_vec->pos[k] = nint;
            mir->cut_vec->pos[kk] = j;
            mir->cut_vec->ind[nint] = k;
            mir->cut_vec->ind[j] = kk;
            temp = mir->cut_vec->val[nint];
            mir->cut_vec->val[nint] = mir->cut_vec->val[j];
            mir->cut_vec->val[j] = temp;
         }
      }
#if _MIR_DEBUG
      ios_check_vec(mir->cut_vec);
#endif
      /* if there is no integer variable, nothing to generate */
      if (nint == 0) goto done;
      /* allocate working arrays */
      u = xcalloc(1+nint, sizeof(double));
      x = xcalloc(1+nint, sizeof(double));
      alpha = xcalloc(1+nint, sizeof(double));
      /* determine u and x */
      for (j = 1; j <= nint; j++)
      {  k = mir->cut_vec->ind[j];
         xassert(m+1 <= k && k <= m+n);
         xassert(mir->isint[k]);
         u[j] = mir->ub[k] - mir->lb[k];
         xassert(u[j] >= 1.0);
         if (mir->subst[k] == 'L')
            x[j] = mir->x[k] - mir->lb[k];
         else if (mir->subst[k] == 'U')
            x[j] = mir->ub[k] - mir->x[k];
         else
            xassert(k != k);
         xassert(x[j] >= -0.001);
         if (x[j] < 0.0) x[j] = 0.0;
      }
      /* compute s = - sum of continuous terms */
      s = 0.0;
      for (j = nint+1; j <= mir->cut_vec->nnz; j++)
      {  double x;
         k = mir->cut_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         /* must be continuous */
         xassert(!mir->isint[k]);
         if (mir->subst[k] == 'L')
         {  xassert(mir->lb[k] != -DBL_MAX);
            kk = mir->vlb[k];
            if (kk == 0)
               x = mir->x[k] - mir->lb[k];
            else
               x = mir->x[k] - mir->lb[k] * mir->x[kk];
         }
         else if (mir->subst[k] == 'U')
         {  xassert(mir->ub[k] != +DBL_MAX);
            kk = mir->vub[k];
            if (kk == 0)
               x = mir->ub[k] - mir->x[k];
            else
               x = mir->ub[k] * mir->x[kk] - mir->x[k];
         }
         else
            xassert(k != k);
         xassert(x >= -0.001);
         if (x < 0.0) x = 0.0;
         s -= mir->cut_vec->val[j] * x;
      }
      xassert(s >= 0.0);
      /* apply heuristic to obtain most violated c-MIR inequality */
      b = mir->cut_rhs;
      r_best = cmir_sep(nint, mir->cut_vec->val, b, u, x, s, alpha,
         &beta, &gamma);
      if (r_best == 0.0) goto skip;
      xassert(r_best > 0.0);
      /* convert to raw cut */
      /* sum alpha[j] * x[j] <= beta + gamma * s */
      for (j = 1; j <= nint; j++)
         mir->cut_vec->val[j] = alpha[j];
      for (j = nint+1; j <= mir->cut_vec->nnz; j++)
      {  k = mir->cut_vec->ind[j];
         if (k <= m+n) mir->cut_vec->val[j] *= gamma;
      }
      mir->cut_rhs = beta;
#if _MIR_DEBUG
      ios_check_vec(mir->cut_vec);
#endif
skip: /* free working arrays */
      xfree(u);
      xfree(x);
      xfree(alpha);
done: return r_best;
}

#if _MIR_DEBUG
static void check_raw_cut(struct MIR *mir, double r_best)
{     /* check raw cut before back bound substitution */
      int m = mir->m;
      int n = mir->n;
      int j, k, kk;
      double r, big, x;
      /* compute the residual r = sum a[k] * x[k] - b and determine
         big = max(1, |a[k]|, |b|) */
      r = 0.0, big = 1.0;
      for (j = 1; j <= mir->cut_vec->nnz; j++)
      {  k = mir->cut_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (mir->subst[k] == 'L')
         {  xassert(mir->lb[k] != -DBL_MAX);
            kk = mir->vlb[k];
            if (kk == 0)
               x = mir->x[k] - mir->lb[k];
            else
               x = mir->x[k] - mir->lb[k] * mir->x[kk];
         }
         else if (mir->subst[k] == 'U')
         {  xassert(mir->ub[k] != +DBL_MAX);
            kk = mir->vub[k];
            if (kk == 0)
               x = mir->ub[k] - mir->x[k];
            else
               x = mir->ub[k] * mir->x[kk] - mir->x[k];
         }
         else
            xassert(k != k);
         r += mir->cut_vec->val[j] * x;
         if (big < fabs(mir->cut_vec->val[j]))
            big = fabs(mir->cut_vec->val[j]);
      }
      r -= mir->cut_rhs;
      if (big < fabs(mir->cut_rhs))
         big = fabs(mir->cut_rhs);
      /* the residual must be close to r_best */
      xassert(fabs(r - r_best) <= 1e-6 * big);
      return;
}
#endif

static void back_subst(struct MIR *mir)
{     /* back substitution of original bounds */
      int m = mir->m;
      int n = mir->n;
      int j, jj, k, kk;
      /* at first, restore bounds of integer variables (because on
         restoring variable bounds of continuous variables we need
         original, not shifted, bounds of integer variables) */
      for (j = 1; j <= mir->cut_vec->nnz; j++)
      {  k = mir->cut_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (!mir->isint[k]) continue; /* skip continuous */
         if (mir->subst[k] == 'L')
         {  /* x'[k] = x[k] - lb[k] */
            xassert(mir->lb[k] != -DBL_MAX);
            xassert(mir->vlb[k] == 0);
            mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
         }
         else if (mir->subst[k] == 'U')
         {  /* x'[k] = ub[k] - x[k] */
            xassert(mir->ub[k] != +DBL_MAX);
            xassert(mir->vub[k] == 0);
            mir->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
            mir->cut_vec->val[j] = - mir->cut_vec->val[j];
         }
         else
            xassert(k != k);
      }
      /* now restore bounds of continuous variables */
      for (j = 1; j <= mir->cut_vec->nnz; j++)
      {  k = mir->cut_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (mir->isint[k]) continue; /* skip integer */
         if (mir->subst[k] == 'L')
         {  /* x'[k] = x[k] - (lower bound) */
            xassert(mir->lb[k] != -DBL_MAX);
            kk = mir->vlb[k];
            if (kk == 0)
            {  /* x'[k] = x[k] - lb[k] */
               mir->cut_rhs += mir->cut_vec->val[j] * mir->lb[k];
            }
            else
            {  /* x'[k] = x[k] - lb[k] * x[kk] */
               jj = mir->cut_vec->pos[kk];
#if 0
               xassert(jj != 0);