glpios06.c

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

         lb = mir->lb[k];
         kk = mir->vlb[k];
         if (kk != 0)
         {  xassert(lb != -DBL_MAX);
            xassert(!mir->isint[k]);
            xassert(mir->isint[kk]);
            lb *= mir->x[kk];
         }
         /* check lower bound */
         if (lb != -DBL_MAX)
         {  eps = 1e-6 * (1.0 + fabs(lb));
            xassert(mir->x[k] >= lb - eps);
         }
         /* determine upper bound */
         ub = mir->ub[k];
         kk = mir->vub[k];
         if (kk != 0)
         {  xassert(ub != +DBL_MAX);
            xassert(!mir->isint[k]);
            xassert(mir->isint[kk]);
            ub *= mir->x[kk];
         }
         /* check upper bound */
         if (ub != +DBL_MAX)
         {  eps = 1e-6 * (1.0 + fabs(ub));
            xassert(mir->x[k] <= ub + eps);
         }
      }
      return;
}
#endif

static void initial_agg_row(glp_tree *tree, struct MIR *mir, int i)
{     /* use original i-th row as initial aggregated constraint */
      glp_prob *mip = tree->mip;
      int m = mir->m;
      GLPAIJ *aij;
      xassert(1 <= i && i <= m);
      xassert(!mir->skip[i]);
      /* mark i-th row in order not to use it in the same aggregated
         constraint */
      mir->skip[i] = 2;
      mir->agg_cnt = 1;
      mir->agg_row[1] = i;
      /* use x[i] - sum a[i,j] * x[m+j] = 0, where x[i] is auxiliary
         variable of row i, x[m+j] are structural variables */
      ios_clear_vec(mir->agg_vec);
      ios_set_vj(mir->agg_vec, i, 1.0);
      for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
         ios_set_vj(mir->agg_vec, m + aij->col->j, - aij->val);
      mir->agg_rhs = 0.0;
#if _MIR_DEBUG
      ios_check_vec(mir->agg_vec);
#endif
      return;
}

#if _MIR_DEBUG
static void check_agg_row(struct MIR *mir)
{     /* check aggregated constraint */
      int m = mir->m;
      int n = mir->n;
      int j, k;
      double r, big;
      /* 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->agg_vec->nnz; j++)
      {  k = mir->agg_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         r += mir->agg_vec->val[j] * mir->x[k];
         if (big < fabs(mir->agg_vec->val[j]))
            big = fabs(mir->agg_vec->val[j]);
      }
      r -= mir->agg_rhs;
      if (big < fabs(mir->agg_rhs))
         big = fabs(mir->agg_rhs);
      /* the residual must be close to zero */
      xassert(fabs(r) <= 1e-6 * big);
      return;
}
#endif

static void subst_fixed_vars(struct MIR *mir)
{     /* substitute fixed variables into aggregated constraint */
      int m = mir->m;
      int n = mir->n;
      int j, k;
      for (j = 1; j <= mir->agg_vec->nnz; j++)
      {  k = mir->agg_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (mir->vlb[k] == 0 && mir->vub[k] == 0 &&
             mir->lb[k] == mir->ub[k])
         {  /* x[k] is fixed */
            mir->agg_rhs -= mir->agg_vec->val[j] * mir->lb[k];
            mir->agg_vec->val[j] = 0.0;
         }
      }
      /* remove terms corresponding to fixed variables */
      ios_clean_vec(mir->agg_vec, DBL_EPSILON);
#if _MIR_DEBUG
      ios_check_vec(mir->agg_vec);
#endif
      return;
}

static void bound_subst_heur(struct MIR *mir)
{     /* bound substitution heuristic */
      int m = mir->m;
      int n = mir->n;
      int j, k, kk;
      double d1, d2;
      for (j = 1; j <= mir->agg_vec->nnz; j++)
      {  k = mir->agg_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (mir->isint[k]) continue; /* skip integer variable */
         /* compute distance from x[k] to its lower bound */
         kk = mir->vlb[k];
         if (kk == 0)
         {  if (mir->lb[k] == -DBL_MAX)
               d1 = DBL_MAX;
            else
               d1 = mir->x[k] - mir->lb[k];
         }
         else
         {  xassert(1 <= kk && kk <= m+n);
            xassert(mir->isint[kk]);
            xassert(mir->lb[k] != -DBL_MAX);
            d1 = mir->x[k] - mir->lb[k] * mir->x[kk];
         }
         /* compute distance from x[k] to its upper bound */
         kk = mir->vub[k];
         if (kk == 0)
         {  if (mir->vub[k] == +DBL_MAX)
               d2 = DBL_MAX;
            else
               d2 = mir->ub[k] - mir->x[k];
         }
         else
         {  xassert(1 <= kk && kk <= m+n);
            xassert(mir->isint[kk]);
            xassert(mir->ub[k] != +DBL_MAX);
            d2 = mir->ub[k] * mir->x[kk] - mir->x[k];
         }
         /* x[k] cannot be free */
         xassert(d1 != DBL_MAX || d2 != DBL_MAX);
         /* choose the bound which is closer to x[k] */
         xassert(mir->subst[k] == '?');
         if (d1 <= d2)
            mir->subst[k] = 'L';
         else
            mir->subst[k] = 'U';
      }
      return;
}

static void build_mod_row(struct MIR *mir)
{     /* substitute bounds and build modified constraint */
      int m = mir->m;
      int n = mir->n;
      int j, jj, k, kk;
      /* initially modified constraint is aggregated constraint */
      ios_copy_vec(mir->mod_vec, mir->agg_vec);
      mir->mod_rhs = mir->agg_rhs;
#if _MIR_DEBUG
      ios_check_vec(mir->mod_vec);
#endif
      /* substitute bounds for continuous variables; note that due to
         substitution of variable bounds additional terms may appear in
         modified constraint */
      for (j = mir->mod_vec->nnz; j >= 1; j--)
      {  k = mir->mod_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (mir->isint[k]) continue; /* skip integer variable */
         if (mir->subst[k] == 'L')
         {  /* x[k] = (lower bound) + x'[k] */
            xassert(mir->lb[k] != -DBL_MAX);
            kk = mir->vlb[k];
            if (kk == 0)
            {  /* x[k] = lb[k] + x'[k] */
               mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
            }
            else
            {  /* x[k] = lb[k] * x[kk] + x'[k] */
               xassert(mir->isint[kk]);
               jj = mir->mod_vec->pos[kk];
               if (jj == 0)
               {  ios_set_vj(mir->mod_vec, kk, 1.0);
                  jj = mir->mod_vec->pos[kk];
                  mir->mod_vec->val[jj] = 0.0;
               }
               mir->mod_vec->val[jj] +=
                  mir->mod_vec->val[j] * mir->lb[k];
            }
         }
         else if (mir->subst[k] == 'U')
         {  /* x[k] = (upper bound) - x'[k] */
            xassert(mir->ub[k] != +DBL_MAX);
            kk = mir->vub[k];
            if (kk == 0)
            {  /* x[k] = ub[k] - x'[k] */
               mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
            }
            else
            {  /* x[k] = ub[k] * x[kk] - x'[k] */
               xassert(mir->isint[kk]);
               jj = mir->mod_vec->pos[kk];
               if (jj == 0)
               {  ios_set_vj(mir->mod_vec, kk, 1.0);
                  jj = mir->mod_vec->pos[kk];
                  mir->mod_vec->val[jj] = 0.0;
               }
               mir->mod_vec->val[jj] +=
                  mir->mod_vec->val[j] * mir->ub[k];
            }
            mir->mod_vec->val[j] = - mir->mod_vec->val[j];
         }
         else
            xassert(k != k);
      }
#if _MIR_DEBUG
      ios_check_vec(mir->mod_vec);
#endif
      /* substitute bounds for integer variables */
      for (j = 1; j <= mir->mod_vec->nnz; j++)
      {  k = mir->mod_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (!mir->isint[k]) continue; /* skip continuous variable */
         xassert(mir->subst[k] == '?');
         xassert(mir->vlb[k] == 0 && mir->vub[k] == 0);
         xassert(mir->lb[k] != -DBL_MAX && mir->ub[k] != +DBL_MAX);
         if (fabs(mir->lb[k]) <= fabs(mir->ub[k]))
         {  /* x[k] = lb[k] + x'[k] */
            mir->subst[k] = 'L';
            mir->mod_rhs -= mir->mod_vec->val[j] * mir->lb[k];
         }
         else
         {  /* x[k] = ub[k] - x'[k] */
            mir->subst[k] = 'U';
            mir->mod_rhs -= mir->mod_vec->val[j] * mir->ub[k];
            mir->mod_vec->val[j] = - mir->mod_vec->val[j];
         }
      }
#if _MIR_DEBUG
      ios_check_vec(mir->mod_vec);
#endif
      return;
}

#if _MIR_DEBUG
static void check_mod_row(struct MIR *mir)
{     /* check modified constraint */
      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->mod_vec->nnz; j++)
      {  k = mir->mod_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         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 = mir->x[k] - mir->lb[k];
            else
               x = mir->x[k] - mir->lb[k] * mir->x[kk];
         }
         else if (mir->subst[k] == 'U')
         {  /* x'[k] = (upper bound) - x[k] */
            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->mod_vec->val[j] * x;
         if (big < fabs(mir->mod_vec->val[j]))
            big = fabs(mir->mod_vec->val[j]);
      }
      r -= mir->mod_rhs;
      if (big < fabs(mir->mod_rhs))
         big = fabs(mir->mod_rhs);
      /* the residual must be close to zero */
      xassert(fabs(r) <= 1e-6 * big);
      return;
}
#endif

/***********************************************************************
*  mir_ineq - construct MIR inequality
*
*  Given the single constraint mixed integer set
*
*                    |N|
*     X = {(x,s) in Z    x R  : sum   a[j] * x[j] <= b + s},
*                    +      +  j in N
*
*  this routine constructs the mixed integer rounding (MIR) inequality
*
*     sum   alpha[j] * x[j] <= beta + gamma * s,
*    j in N
*
*  which is valid for X.
*
*  If the MIR inequality has been successfully constructed, the routine
*  returns zero. Otherwise, if b is close to nearest integer, there may
*  be numeric difficulties due to big coefficients; so in this case the
*  routine returns non-zero. */

static int mir_ineq(const int n, const double a[], const double b,
      double alpha[], double *beta, double *gamma)
{     int j;
      double f, t;
      if (fabs(b - floor(b + .5)) < 0.01)
         return 1;
      f = b - floor(b);
      for (j = 1; j <= n; j++)
      {  t = (a[j] - floor(a[j])) - f;
         if (t <= 0.0)
            alpha[j] = floor(a[j]);
         else
            alpha[j] = floor(a[j]) + t / (1.0 - f);
      }
      *beta = floor(b);
      *gamma = 1.0 / (1.0 - f);
      return 0;
}

/***********************************************************************
*  cmir_ineq - construct c-MIR inequality
*
*  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]},
*
*  a subset C of variables to be complemented, and a divisor delta > 0,
*  this routine constructs the complemented MIR (c-MIR) inequality
*
*     sum   alpha[j] * x[j] <= beta + gamma * s,
*    j in N
*                      MK
*  which is valid for X  .
*
*  If the c-MIR inequality has been successfully constructed, the
*  routine returns zero. Otherwise, if there is a risk of numerical
*  difficulties due to big coefficients (see comments to the routine
*  mir_ineq), the routine cmir_ineq returns non-zero. */

static int cmir_ineq(const int n, const double a[], const double b,
      const double u[], const char cset[], const double delta,
      double alpha[], double *beta, double *gamma)