glpios06.c

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

#else
               if (jj == 0)
               {  ios_set_vj(mir->cut_vec, kk, 1.0);
                  jj = mir->cut_vec->pos[kk];
                  xassert(jj != 0);
                  mir->cut_vec->val[jj] = 0.0;
               }
#endif
               mir->cut_vec->val[jj] -= mir->cut_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->cut_rhs -= mir->cut_vec->val[j] * mir->ub[k];
            }
            else
            {  /* x'[k] = ub[k] * x[kk] - x[k] */
               jj = mir->cut_vec->pos[kk];
               if (jj == 0)
               {  ios_set_vj(mir->cut_vec, kk, 1.0);
                  jj = mir->cut_vec->pos[kk];
                  xassert(jj != 0);
                  mir->cut_vec->val[jj] = 0.0;
               }
               mir->cut_vec->val[jj] += mir->cut_vec->val[j] *
                  mir->ub[k];
            }
            mir->cut_vec->val[j] = - mir->cut_vec->val[j];
         }
         else
            xassert(k != k);
      }
#if _MIR_DEBUG
      ios_check_vec(mir->cut_vec);
#endif
      return;
}

#if _MIR_DEBUG
static void check_cut_row(struct MIR *mir, double r_best)
{     /* check the cut after back bound substitution or elimination of
         auxiliary variables */
      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->cut_vec->nnz; j++)
      {  k = mir->cut_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         r += mir->cut_vec->val[j] * mir->x[k];
         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 subst_aux_vars(glp_tree *tree, struct MIR *mir)
{     /* final substitution to eliminate auxiliary variables */
      glp_prob *mip = tree->mip;
      int m = mir->m;
      int n = mir->n;
      GLPAIJ *aij;
      int j, k, kk, jj;
      for (j = mir->cut_vec->nnz; j >= 1; j--)
      {  k = mir->cut_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (k > m) continue; /* skip structurals */
         for (aij = mip->row[k]->ptr; aij != NULL; aij = aij->r_next)
         {  kk = m + aij->col->j; /* structural */
            jj = mir->cut_vec->pos[kk];
            if (jj == 0)
            {  ios_set_vj(mir->cut_vec, kk, 1.0);
               jj = mir->cut_vec->pos[kk];
               mir->cut_vec->val[jj] = 0.0;
            }
            mir->cut_vec->val[jj] += mir->cut_vec->val[j] * aij->val;
         }
         mir->cut_vec->val[j] = 0.0;
      }
      ios_clean_vec(mir->cut_vec, 0.0);
      return;
}

static void add_cut(glp_tree *tree, struct MIR *mir)
{     /* add constructed cut inequality to the cut pool */
      int m = mir->m;
      int n = mir->n;
      int j, k, len;
      int *ind = xcalloc(1+n, sizeof(int));
      double *val = xcalloc(1+n, sizeof(double));
      len = 0;
      for (j = mir->cut_vec->nnz; j >= 1; j--)
      {  k = mir->cut_vec->ind[j];
         xassert(m+1 <= k && k <= m+n);
         len++, ind[len] = k - m, val[len] = mir->cut_vec->val[j];
      }
#if 0
      ios_add_cut_row(tree, pool, GLP_RF_MIR, len, ind, val, GLP_UP,
         mir->cut_rhs);
#else
      glp_ios_add_row(tree, NULL, GLP_RF_MIR, 0, len, ind, val, GLP_UP,
         mir->cut_rhs);
#endif
      xfree(ind);
      xfree(val);
      return;
}

static int aggregate_row(glp_tree *tree, struct MIR *mir)
{     /* try to aggregate another row */
      glp_prob *mip = tree->mip;
      int m = mir->m;
      int n = mir->n;
      GLPAIJ *aij;
      IOSVEC *v;
      int ii, j, jj, k, kk, kappa = 0, ret = 0;
      double d1, d2, d, d_max = 0.0;
      /* choose appropriate structural variable in the aggregated row
         to be substituted */
      for (j = 1; j <= mir->agg_vec->nnz; j++)
      {  k = mir->agg_vec->ind[j];
         xassert(1 <= k && k <= m+n);
         if (k <= m) continue; /* skip auxiliary var */
         if (mir->isint[k]) continue; /* skip integer var */
         if (fabs(mir->agg_vec->val[j]) < 0.001) continue;
         /* 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);
         /* d = min(d1, d2) */
         d = (d1 <= d2 ? d1 : d2);
         xassert(d != DBL_MAX);
         /* should not be close to corresponding bound */
         if (d < 0.001) continue;
         if (d_max < d) d_max = d, kappa = k;
      }
      if (kappa == 0)
      {  /* nothing chosen */
         ret = 1;
         goto done;
      }
      /* x[kappa] has been chosen */
      xassert(m+1 <= kappa && kappa <= m+n);
      xassert(!mir->isint[kappa]);
      /* find another row, which have not been used yet, to eliminate
         x[kappa] from the aggregated row */
      for (ii = 1; ii <= m; ii++)
      {  if (mir->skip[ii]) continue;
         for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
            if (aij->col->j == kappa - m) break;
         if (aij != NULL && fabs(aij->val) >= 0.001) break;
      }
      if (ii > m)
      {  /* nothing found */
         ret = 2;
         goto done;
      }
      /* row ii has been found; include it in the aggregated list */
      mir->agg_cnt++;
      xassert(mir->agg_cnt <= MAXAGGR);
      mir->agg_row[mir->agg_cnt] = ii;
      mir->skip[ii] = 2;
      /* v := new row */
      v = ios_create_vec(m+n);
      ios_set_vj(v, ii, 1.0);
      for (aij = mip->row[ii]->ptr; aij != NULL; aij = aij->r_next)
         ios_set_vj(v, m + aij->col->j, - aij->val);
#if _MIR_DEBUG
      ios_check_vec(v);
#endif
      /* perform gaussian elimination to remove x[kappa] */
      j = mir->agg_vec->pos[kappa];
      xassert(j != 0);
      jj = v->pos[kappa];
      xassert(jj != 0);
      ios_linear_comb(mir->agg_vec,
         - mir->agg_vec->val[j] / v->val[jj], v);
      ios_delete_vec(v);
      ios_set_vj(mir->agg_vec, kappa, 0.0);
#if _MIR_DEBUG
      ios_check_vec(mir->agg_vec);
#endif
done: return ret;
}

void ios_mir_gen(glp_tree *tree, void *gen)
{     /* main routine to generate MIR cuts */
      glp_prob *mip = tree->mip;
      struct MIR *mir = gen;
      int m = mir->m;
      int n = mir->n;
      int i;
      double r_best;
      xassert(mip->m >= m);
      xassert(mip->n == n);
      /* obtain current point */
      get_current_point(tree, mir);
#if _MIR_DEBUG
      /* check current point */
      check_current_point(mir);
#endif
      /* reset bound substitution flags */
      memset(&mir->subst[1], '?', m+n);
      /* try to generate a set of violated MIR cuts */
      for (i = 1; i <= m; i++)
      {  if (mir->skip[i]) continue;
         /* use original i-th row as initial aggregated constraint */
         initial_agg_row(tree, mir, i);
loop:    ;
#if _MIR_DEBUG
         /* check aggregated row */
         check_agg_row(mir);
#endif
         /* substitute fixed variables into aggregated constraint */
         subst_fixed_vars(mir);
#if _MIR_DEBUG
         /* check aggregated row */
         check_agg_row(mir);
#endif
#if _MIR_DEBUG
         /* check bound substitution flags */
         {  int k;
            for (k = 1; k <= m+n; k++)
               xassert(mir->subst[k] == '?');
         }
#endif
         /* apply bound substitution heuristic */
         bound_subst_heur(mir);
         /* substitute bounds and build modified constraint */
         build_mod_row(mir);
#if _MIR_DEBUG
         /* check modified row */
         check_mod_row(mir);
#endif
         /* try to generate violated c-MIR cut for modified row */
         r_best = generate(mir);
         if (r_best > 0.0)
         {  /* success */
#if _MIR_DEBUG
            /* check raw cut before back bound substitution */
            check_raw_cut(mir, r_best);
#endif
            /* back substitution of original bounds */
            back_subst(mir);
#if _MIR_DEBUG
            /* check the cut after back bound substitution */
            check_cut_row(mir, r_best);
#endif
            /* final substitution to eliminate auxiliary variables */
            subst_aux_vars(tree, mir);
#if _MIR_DEBUG
            /* check the cut after elimination of auxiliaries */
            check_cut_row(mir, r_best);
#endif
            /* add constructed cut inequality to the cut pool */
            add_cut(tree, mir);
         }
         /* reset bound substitution flags */
         {  int j, k;
            for (j = 1; j <= mir->mod_vec->nnz; j++)
            {  k = mir->mod_vec->ind[j];
               xassert(1 <= k && k <= m+n);
               xassert(mir->subst[k] != '?');
               mir->subst[k] = '?';
            }
         }
         if (r_best == 0.0)
         {  /* failure */
            if (mir->agg_cnt < MAXAGGR)
            {  /* try to aggregate another row */
               if (aggregate_row(tree, mir) == 0) goto loop;
            }
         }
         /* unmark rows used in the aggregated constraint */
         {  int k, ii;
            for (k = 1; k <= mir->agg_cnt; k++)
            {  ii = mir->agg_row[k];
               xassert(1 <= ii && ii <= m);
               xassert(mir->skip[ii] == 2);
               mir->skip[ii] = 0;
            }
         }
      }
      return;
}

/***********************************************************************
*  NAME
*
*  ios_mir_term - terminate MIR cut generator
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  void ios_mir_term(void *gen);
*
*  DESCRIPTION
*
*  The routine ios_mir_term deletes the MIR cut generator working area
*  freeing all the memory allocated to it. */

void ios_mir_term(void *gen)
{     struct MIR *mir = gen;
      xfree(mir->skip);
      xfree(mir->isint);
      xfree(mir->lb);
      xfree(mir->vlb);
      xfree(mir->ub);
      xfree(mir->vub);
      xfree(mir->x);
      xfree(mir->agg_row);
      ios_delete_vec(mir->agg_vec);
      xfree(mir->subst);
      ios_delete_vec(mir->mod_vec);
      ios_delete_vec(mir->cut_vec);
      xfree(mir);
      return;
}

/* eof */