glpios01.c

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

      }
      m = tree->orig_m;
      /* restore original attributes of rows and columns */
      xassert(m == tree->orig_m);
      xassert(n == tree->n);
      for (i = 1; i <= m; i++)
      {  glp_set_row_bnds(mip, i, tree->orig_type[i],
            tree->orig_lb[i], tree->orig_ub[i]);
         glp_set_row_stat(mip, i, tree->orig_stat[i]);
         mip->row[i]->prim = tree->orig_prim[i];
         mip->row[i]->dual = tree->orig_dual[i];
      }
      for (j = 1; j <= n; j++)
      {  glp_set_col_bnds(mip, j, tree->orig_type[m+j],
            tree->orig_lb[m+j], tree->orig_ub[m+j]);
         glp_set_col_stat(mip, j, tree->orig_stat[m+j]);
         mip->col[j]->prim = tree->orig_prim[m+j];
         mip->col[j]->dual = tree->orig_dual[m+j];
      }
      mip->pbs_stat = mip->dbs_stat = GLP_FEAS;
      mip->obj_val = tree->orig_obj;
      /* delete the branch-and-bound tree */
      xassert(tree->local != NULL);
      ios_delete_pool(tree, tree->local);
      dmp_delete_pool(tree->pool);
      xfree(tree->orig_type);
      xfree(tree->orig_lb);
      xfree(tree->orig_ub);
      xfree(tree->orig_stat);
      xfree(tree->orig_prim);
      xfree(tree->orig_dual);
      xfree(tree->slot);
      if (tree->root_type != NULL) xfree(tree->root_type);
      if (tree->root_lb != NULL) xfree(tree->root_lb);
      if (tree->root_ub != NULL) xfree(tree->root_ub);
      if (tree->root_stat != NULL) xfree(tree->root_stat);
      xfree(tree->non_int);
      xfree(tree->n_ref);
      xfree(tree->c_ref);
      xfree(tree->j_ref);
      if (tree->pcost != NULL) ios_pcost_free(tree);
      xfree(tree->iwrk);
      xfree(tree->dwrk);
      scg_delete_graph(tree->g);
      if (tree->pred_type != NULL) xfree(tree->pred_type);
      if (tree->pred_lb != NULL) xfree(tree->pred_lb);
      if (tree->pred_ub != NULL) xfree(tree->pred_ub);
      if (tree->pred_stat != NULL) xfree(tree->pred_stat);
#if 0
      xassert(tree->cut_gen == NULL);
#endif
      xassert(tree->mir_gen == NULL);
      xassert(tree->clq_gen == NULL);
      xfree(tree);
      mip->tree = NULL;
      return;
}

/***********************************************************************
*  NAME
*
*  ios_eval_degrad - estimate obj. degrad. for down- and up-branches
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  void ios_eval_degrad(glp_tree *tree, int j, double *dn, double *up);
*
*  DESCRIPTION
*
*  Given optimal basis to LP relaxation of the current subproblem the
*  routine ios_eval_degrad performs the dual ratio test to compute the
*  objective values in the adjacent basis for down- and up-branches,
*  which are stored in locations *dn and *up, assuming that x[j] is a
*  variable chosen to branch upon. */

void ios_eval_degrad(glp_tree *tree, int j, double *dn, double *up)
{     glp_prob *mip = tree->mip;
      int m = mip->m, n = mip->n;
      int len, kase, k, t, stat;
      double alfa, beta, gamma, delta, dz;
      int *ind = tree->iwrk;
      double *val = tree->dwrk;
      /* current basis must be optimal */
      xassert(glp_get_status(mip) == GLP_OPT);
      /* basis factorization must exist */
      xassert(glp_bf_exists(mip));
      /* obtain (fractional) value of x[j] in optimal basic solution
         to LP relaxation of the current subproblem */
      xassert(1 <= j && j <= n);
      beta = mip->col[j]->prim;
      /* since the value of x[j] is fractional, it is basic; compute
         corresponding row of the simplex table */
      len = lpx_eval_tab_row(mip, m+j, ind, val);
      /* kase < 0 means down-branch; kase > 0 means up-branch */
      for (kase = -1; kase <= +1; kase += 2)
      {  /* for down-branch we introduce new upper bound floor(beta)
            for x[j]; similarly, for up-branch we introduce new lower
            bound ceil(beta) for x[j]; in the current basis this new
            upper/lower bound is violated, so in the adjacent basis
            x[j] will leave the basis and go to its new upper/lower
            bound; we need to know which non-basic variable x[k] should
            enter the basis to keep dual feasibility */
         k = lpx_dual_ratio_test(mip, len, ind, val, kase, 1e-7);
         /* if no variable has been chosen, current basis being primal
            infeasible due to the new upper/lower bound of x[j] is dual
            unbounded, therefore, LP relaxation to corresponding branch
            has no primal feasible solution */
         if (k == 0)
         {  if (mip->dir == GLP_MIN)
            {  if (kase < 0)
                  *dn = +DBL_MAX;
               else
                  *up = +DBL_MAX;
            }
            else if (mip->dir == GLP_MAX)
            {  if (kase < 0)
                  *dn = -DBL_MAX;
               else
                  *up = -DBL_MAX;
            }
            else
               xassert(mip != mip);
            continue;
         }
         xassert(1 <= k && k <= m+n);
         /* row of the simplex table corresponding to specified basic
            variable x[j] is the following:
               x[j] = ... + alfa * x[k] + ... ;
            we need to know influence coefficient, alfa, at non-basic
            variable x[k] chosen with the dual ratio test */
         for (t = 1; t <= len; t++)
            if (ind[t] == k) break;
         xassert(1 <= t && t <= len);
         alfa = val[t];
         /* determine status and reduced cost of variable x[k] */
         if (k <= m)
         {  stat = mip->row[k]->stat;
            gamma = mip->row[k]->dual;
         }
         else
         {  stat = mip->col[k-m]->stat;
            gamma = mip->col[k-m]->dual;
         }
         /* x[k] cannot be basic or fixed non-basic */
         xassert(stat == GLP_NL || stat == GLP_NU || stat == GLP_NF);
         /* if the current basis is dual degenerative, some reduced
            costs, which are close to zero, may have wrong sign due to
            round-off errors, so correct the sign of gamma */
         if (mip->dir == GLP_MIN)
         {  if (stat == GLP_NL && gamma < 0.0 ||
                stat == GLP_NU && gamma > 0.0 ||
                stat == GLP_NF) gamma = 0.0;
         }
         else if (mip->dir == GLP_MAX)
         {  if (stat == GLP_NL && gamma > 0.0 ||
                stat == GLP_NU && gamma < 0.0 ||
                stat == GLP_NF) gamma = 0.0;
         }
         else
            xassert(mip != mip);
         /* determine the change of x[j] in the adjacent basis:
            delta x[j] = new x[j] - old x[j] */
         delta = (kase < 0 ? floor(beta) : ceil(beta)) - beta;
         /* compute the change of x[k] in the adjacent basis:
            delta x[k] = new x[k] - old x[k] = delta x[j] / alfa */
         delta /= alfa;
         /* compute the change of the objective in the adjacent basis:
            delta z = new z - old z = gamma * delta x[k] */
         dz = gamma * delta;
         if (mip->dir == GLP_MIN)
            xassert(dz >= 0.0);
         else if (mip->dir == GLP_MAX)
            xassert(dz <= 0.0);
         else
            xassert(mip != mip);
         /* compute the new objective value in the adjacent basis:
            new z = old z + delta z */
         if (kase < 0)
            *dn = mip->obj_val + dz;
         else
            *up = mip->obj_val + dz;
      }
      /*xprintf("obj = %g; dn = %g; up = %g\n",
         mip->obj_val, *dn, *up);*/
      return;
}

/***********************************************************************
*  NAME
*
*  ios_round_bound - improve local bound by rounding
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  double ios_round_bound(glp_tree *tree, double bound);
*
*  RETURNS
*
*  For the given local bound for any integer feasible solution to the
*  current subproblem the routine ios_round_bound returns an improved
*  local bound for the same integer feasible solution.
*
*  BACKGROUND
*
*  Let the current subproblem has the following objective function:
*
*     z =   sum  c[j] * x[j] + s >= b,                               (1)
*         j in J
*
*  where J = {j: c[j] is non-zero and integer, x[j] is integer}, s is
*  the sum of terms corresponding to fixed variables, b is an initial
*  local bound (minimization).
*
*  From (1) it follows that:
*
*     d *  sum  (c[j] / d) * x[j] + s >= b,                          (2)
*        j in J
*
*  or, equivalently,
*
*     sum  (c[j] / d) * x[j] >= (b - s) / d = h,                     (3)
*   j in J
*
*  where d = gcd(c[j]). Since the left-hand side of (3) is integer,
*  h = (b - s) / d can be rounded up to the nearest integer:
*
*     h' = ceil(h) = (b' - s) / d,                                   (4)
*
*  that gives an rounded, improved local bound:
*
*     b' = d * h' + s.                                               (5)
*
*  In case of maximization '>=' in (1) should be replaced by '<=' that
*  leads to the following formula:
*
*     h' = floor(h) = (b' - s) / d,                                  (6)
*
*  which should used in the same way as (4).
*
*  NOTE: If b is a valid local bound for a child of the current
*        subproblem, b' is also valid for that child subproblem. */

double ios_round_bound(glp_tree *tree, double bound)
{     glp_prob *mip = tree->mip;
      int n = mip->n;
      int d, j, nn, *c = tree->iwrk;
      double s, h;
      /* determine c[j] and compute s */
      nn = 0, s = mip->c0, d = 0;
      for (j = 1; j <= n; j++)
      {  GLPCOL *col = mip->col[j];
         if (col->coef == 0.0) continue;
         if (col->type == GLP_FX)
         {  /* fixed variable */
            s += col->coef * col->prim;
         }
         else
         {  /* non-fixed variable */
            if (col->kind != GLP_IV) goto skip;
            if (col->coef != floor(col->coef)) goto skip;
            if (fabs(col->coef) <= (double)INT_MAX)
               c[++nn] = (int)fabs(col->coef);
            else
               d = 1;
         }
      }
      /* compute d = gcd(c[1],...c[nn]) */
      if (d == 0)
      {  if (nn == 0) goto skip;
         d = gcdn(nn, c);
      }
      xassert(d > 0);
      /* compute new local bound */
      if (mip->dir == GLP_MIN)
      {  if (bound != +DBL_MAX)
         {  h = (bound - s) / (double)d;
            if (h >= floor(h) + 0.001)
            {  /* round up */
               h = ceil(h);
               /*xprintf("d = %d; old = %g; ", d, bound);*/
               bound = (double)d * h + s;
               /*xprintf("new = %g\n", bound);*/
            }
         }
      }
      else if (mip->dir == GLP_MAX)
      {  if (bound != -DBL_MAX)
         {  h = (bound - s) / (double)d;
            if (h <= ceil(h) - 0.001)
            {  /* round down */
               h = floor(h);
               bound = (double)d * h + s;
            }
         }
      }
      else
         xassert(mip != mip);
skip: return bound;
}

/***********************************************************************
*  NAME
*
*  ios_is_hopeful - check if subproblem is hopeful
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  int ios_is_hopeful(glp_tree *tree, double bound);
*
*  DESCRIPTION
*
*  Given the local bound of a subproblem the routine ios_is_hopeful
*  checks if the subproblem can have an integer optimal solution which
*  is better than the best one currently known.
*
*  RETURNS
*
*  If the subproblem can have a better integer optimal solution, the
*  routine returns non-zero; otherwise, if the corresponding branch can
*  be pruned, the routine returns zero. */

int ios_is_hopeful(glp_tree *tree, double bound)
{     glp_prob *mip = tree->mip;
      int ret = 1;
      double eps;
      if (mip->mip_stat == GLP_FEAS)
      {  eps = tree->parm->tol_obj * (1.0 + fabs(mip->mip_obj));
         switch (mip->dir)
         {  case GLP_MIN:
               if (bound >= mip->mip_obj - eps) ret = 0;
               break;
            case GLP_MAX:
               if (bound <= mip->mip_obj + eps) ret = 0;
               break;
            default:
               xassert(mip != mip);
         }
      }
      else
      {  switch (mip->dir)
         {  case GLP_MIN:
               if (bound == +DBL_MAX) ret = 0;
               break;
            case GLP_MAX:
               if (bound == -DBL_MAX) ret = 0;
               break;
            default:
               xassert(mip != mip);
         }
      }
      return ret;
}

/***********************************************************************
*  NAME
*
*  ios_best_node - find active node with best local bound
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  int ios_best_node(glp_tree *tree);
*
*  DESCRIPTION
*
*  The routine ios_best_node finds an active node whose local bound is
*  best among other active nodes.
*
*  It is understood that the integer optimal solution of the original
*  mip problem cannot be better than the best bound, so the best bound
*  is an lower (minimization) or upper (maximization) global bound for
*  the original problem.
*
*  RETURNS
*
*  The routine ios_best_node returns the subproblem reference number
*  for the best node. However, if the tree is empty, it returns zero. */

int ios_best_node(glp_tree *tree)
{     IOSNPD *node, *best = NULL;
      switch (tree->mip->dir)
      {  case GLP_MIN:
            /* minimization */
            for (node = tree->head; node != NULL; node = node->next)
               if (best == NULL || best->bound > node->bound)
                  best = node;
            break;
         case GLP_MAX:
            /* maximization */
            for (node = tree->head; node != NULL; node = node->next)
               if (best == NULL || best->bound < node->bound)