glpios06.c

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

/* glpios06.c (MIR cut generator) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2000,01,02,03,04,05,06,07,08,2009 Andrew Makhorin,
*  Department for Applied Informatics, Moscow Aviation Institute,
*  Moscow, Russia. All rights reserved. E-mail: <[e-mail]>.
*
*  GLPK is free software: you can redistribute it and/or modify it
*  under the terms of the GNU General Public License as published by
*  the Free Software Foundation, either version 3 of the License, or
*  (at your option) any later version.
*
*  GLPK is distributed in the hope that it will be useful, but WITHOUT
*  ANY WARRANTY; without even the implied warranty of MERCHANTABILITY
*  or FITNESS FOR A PARTICULAR PURPOSE. See the GNU General Public
*  License for more details.
*
*  You should have received a copy of the GNU General Public License
*  along with GLPK. If not, see <http://www.gnu.org/licenses/>.
***********************************************************************/

#include "glpios.h"

#define _MIR_DEBUG 0

#define MAXAGGR 5
/* maximal number of rows which can be aggregated */

struct MIR
{     /* MIR cut generator working area */
      /*--------------------------------------------------------------*/
      /* global information valid for the root subproblem */
      int m;
      /* number of rows (in the root subproblem) */
      int n;
      /* number of columns */
      char *skip; /* char skip[1+m]; */
      /* skip[i], 1 <= i <= m, is a flag that means that row i should
         not be used because (1) it is not suitable, or (2) because it
         has been used in the aggregated constraint */
      char *isint; /* char isint[1+m+n]; */
      /* isint[k], 1 <= k <= m+n, is a flag that means that variable
         x[k] is integer (otherwise, continuous) */
      double *lb; /* double lb[1+m+n]; */
      /* lb[k], 1 <= k <= m+n, is lower bound of x[k]; -DBL_MAX means
         that x[k] has no lower bound */
      int *vlb; /* int vlb[1+m+n]; */
      /* vlb[k] = k', 1 <= k <= m+n, is the number of integer variable,
         which defines variable lower bound x[k] >= lb[k] * x[k']; zero
         means that x[k] has simple lower bound */
      double *ub; /* double ub[1+m+n]; */
      /* ub[k], 1 <= k <= m+n, is upper bound of x[k]; +DBL_MAX means
         that x[k] has no upper bound */
      int *vub; /* int vub[1+m+n]; */
      /* vub[k] = k', 1 <= k <= m+n, is the number of integer variable,
         which defines variable upper bound x[k] <= ub[k] * x[k']; zero
         means that x[k] has simple upper bound */
      /*--------------------------------------------------------------*/
      /* current (fractional) point to be separated */
      double *x; /* double x[1+m+n]; */
      /* x[k] is current value of auxiliary (1 <= k <= m) or structural
         (m+1 <= k <= m+n) variable */
      /*--------------------------------------------------------------*/
      /* aggregated constraint sum a[k] * x[k] = b, which is a linear
         combination of original constraints transformed to equalities
         by introducing auxiliary variables */
      int agg_cnt;
      /* number of rows (original constraints) used to build aggregated
         constraint, 1 <= agg_cnt <= MAXAGGR */
      int *agg_row; /* int agg_row[1+MAXAGGR]; */
      /* agg_row[k], 1 <= k <= agg_cnt, is the row number used to build
         aggregated constraint */
      IOSVEC *agg_vec; /* IOSVEC agg_vec[1:m+n]; */
      /* sparse vector of aggregated constraint coefficients, a[k] */
      double agg_rhs;
      /* right-hand side of the aggregated constraint, b */
      /*--------------------------------------------------------------*/
      /* bound substitution flags for modified constraint */
      char *subst; /* char subst[1+m+n]; */
      /* subst[k], 1 <= k <= m+n, is a bound substitution flag used for
         variable x[k]:
         '?' - x[k] is missing in modified constraint
         'L' - x[k] = (lower bound) + x'[k]
         'U' - x[k] = (upper bound) - x'[k] */
      /*--------------------------------------------------------------*/
      /* modified constraint sum a'[k] * x'[k] = b', where x'[k] >= 0,
         derived from aggregated constraint by substituting bounds;
         note that due to substitution of variable bounds there may be
         additional terms in the modified constraint */
      IOSVEC *mod_vec; /* IOSVEC mod_vec[1:m+n]; */
      /* sparse vector of modified constraint coefficients, a'[k] */
      double mod_rhs;
      /* right-hand side of the modified constraint, b' */
      /*--------------------------------------------------------------*/
      /* cutting plane sum alpha[k] * x[k] <= beta */
      IOSVEC *cut_vec; /* IOSVEC cut_vec[1:m+n]; */
      /* sparse vector of cutting plane coefficients, alpha[k] */
      double cut_rhs;
      /* right-hand size of the cutting plane, beta */
};

/***********************************************************************
*  NAME
*
*  ios_mir_init - initialize MIR cut generator
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  void *ios_mir_init(glp_tree *tree);
*
*  DESCRIPTION
*
*  The routine ios_mir_init initializes the MIR cut generator assuming
*  that the current subproblem is the root subproblem.
*
*  RETURNS
*
*  The routine ios_mir_init returns a pointer to the MIR cut generator
*  working area. */

static void set_row_attrib(glp_tree *tree, struct MIR *mir)
{     /* set global row attributes */
      glp_prob *mip = tree->mip;
      int m = mir->m;
      int k;
      for (k = 1; k <= m; k++)
      {  GLPROW *row = mip->row[k];
         mir->skip[k] = 0;
         mir->isint[k] = 0;
         switch (row->type)
         {  case GLP_FR:
               mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
            case GLP_LO:
               mir->lb[k] = row->lb, mir->ub[k] = +DBL_MAX; break;
            case GLP_UP:
               mir->lb[k] = -DBL_MAX, mir->ub[k] = row->ub; break;
            case GLP_DB:
               mir->lb[k] = row->lb, mir->ub[k] = row->ub; break;
            case GLP_FX:
               mir->lb[k] = mir->ub[k] = row->lb; break;
            default:
               xassert(row != row);
         }
         mir->vlb[k] = mir->vub[k] = 0;
      }
      return;
}

static void set_col_attrib(glp_tree *tree, struct MIR *mir)
{     /* set global column attributes */
      glp_prob *mip = tree->mip;
      int m = mir->m;
      int n = mir->n;
      int k;
      for (k = m+1; k <= m+n; k++)
      {  GLPCOL *col = mip->col[k-m];
         switch (col->kind)
         {  case GLP_CV:
               mir->isint[k] = 0; break;
            case GLP_IV:
               mir->isint[k] = 1; break;
            default:
               xassert(col != col);
         }
         switch (col->type)
         {  case GLP_FR:
               mir->lb[k] = -DBL_MAX, mir->ub[k] = +DBL_MAX; break;
            case GLP_LO:
               mir->lb[k] = col->lb, mir->ub[k] = +DBL_MAX; break;
            case GLP_UP:
               mir->lb[k] = -DBL_MAX, mir->ub[k] = col->ub; break;
            case GLP_DB:
               mir->lb[k] = col->lb, mir->ub[k] = col->ub; break;
            case GLP_FX:
               mir->lb[k] = mir->ub[k] = col->lb; break;
            default:
               xassert(col != col);
         }
         mir->vlb[k] = mir->vub[k] = 0;
      }
      return;
}

static void set_var_bounds(glp_tree *tree, struct MIR *mir)
{     /* set variable bounds */
      glp_prob *mip = tree->mip;
      int m = mir->m;
      GLPAIJ *aij;
      int i, k1, k2;
      double a1, a2;
      for (i = 1; i <= m; i++)
      {  /* we need the row to be '>= 0' or '<= 0' */
         if (!(mir->lb[i] == 0.0 && mir->ub[i] == +DBL_MAX ||
               mir->lb[i] == -DBL_MAX && mir->ub[i] == 0.0)) continue;
         /* take first term */
         aij = mip->row[i]->ptr;
         if (aij == NULL) continue;
         k1 = m + aij->col->j, a1 = aij->val;
         /* take second term */
         aij = aij->r_next;
         if (aij == NULL) continue;
         k2 = m + aij->col->j, a2 = aij->val;
         /* there must be only two terms */
         if (aij->r_next != NULL) continue;
         /* interchange terms, if needed */
         if (!mir->isint[k1] && mir->isint[k2])
            ;
         else if (mir->isint[k1] && !mir->isint[k2])
         {  k2 = k1, a2 = a1;
            k1 = m + aij->col->j, a1 = aij->val;
         }
         else
         {  /* both terms are either continuous or integer */
            continue;
         }
         /* x[k2] should be double-bounded */
         if (mir->lb[k2] == -DBL_MAX || mir->ub[k2] == +DBL_MAX ||
             mir->lb[k2] == mir->ub[k2]) continue;
         /* change signs, if necessary */
         if (mir->ub[i] == 0.0) a1 = - a1, a2 = - a2;
         /* now the row has the form a1 * x1 + a2 * x2 >= 0, where x1
            is continuous, x2 is integer */
         if (a1 > 0.0)
         {  /* x1 >= - (a2 / a1) * x2 */
            if (mir->vlb[k1] == 0)
            {  /* set variable lower bound for x1 */
               mir->lb[k1] = - a2 / a1;
               mir->vlb[k1] = k2;
               /* the row should not be used */
               mir->skip[i] = 1;
            }
         }
         else /* a1 < 0.0 */
         {  /* x1 <= - (a2 / a1) * x2 */
            if (mir->vub[k1] == 0)
            {  /* set variable upper bound for x1 */
               mir->ub[k1] = - a2 / a1;
               mir->vub[k1] = k2;
               /* the row should not be used */
               mir->skip[i] = 1;
            }
         }
      }
      return;
}

static void mark_useless_rows(glp_tree *tree, struct MIR *mir)
{     /* mark rows which should not be used */
      glp_prob *mip = tree->mip;
      int m = mir->m;
      GLPAIJ *aij;
      int i, k, nv;
      for (i = 1; i <= m; i++)
      {  /* free rows should not be used */
         if (mir->lb[i] == -DBL_MAX && mir->ub[i] == +DBL_MAX)
         {  mir->skip[i] = 1;
            continue;
         }
         nv = 0;
         for (aij = mip->row[i]->ptr; aij != NULL; aij = aij->r_next)
         {  k = m + aij->col->j;
            /* rows with free variables should not be used */
            if (mir->lb[k] == -DBL_MAX && mir->ub[k] == +DBL_MAX)
            {  mir->skip[i] = 1;
               break;
            }
            /* rows with integer variables having infinite (lower or
               upper) bound should not be used */
            if (mir->isint[k] && mir->lb[k] == -DBL_MAX ||
                mir->isint[k] && mir->ub[k] == +DBL_MAX)
            {  mir->skip[i] = 1;
               break;
            }
            /* count non-fixed variables */
            if (!(mir->vlb[k] == 0 && mir->vub[k] == 0 &&
                  mir->lb[k] == mir->ub[k])) nv++;
         }
         /* rows with all variables fixed should not be used */
         if (nv == 0)
         {  mir->skip[i] = 1;
            continue;
         }
      }
      return;
}

void *ios_mir_init(glp_tree *tree)
{     /* initialize MIR cut generator */
      glp_prob *mip = tree->mip;
      int m = mip->m;
      int n = mip->n;
      struct MIR *mir;
#if _MIR_DEBUG
      xprintf("ios_mir_init: warning: debug mode enabled\n");
#endif
      /* allocate working area */
      mir = xmalloc(sizeof(struct MIR));
      mir->m = m;
      mir->n = n;
      mir->skip = xcalloc(1+m, sizeof(char));
      mir->isint = xcalloc(1+m+n, sizeof(char));
      mir->lb = xcalloc(1+m+n, sizeof(double));
      mir->vlb = xcalloc(1+m+n, sizeof(int));
      mir->ub = xcalloc(1+m+n, sizeof(double));
      mir->vub = xcalloc(1+m+n, sizeof(int));
      mir->x = xcalloc(1+m+n, sizeof(double));
      mir->agg_row = xcalloc(1+MAXAGGR, sizeof(int));
      mir->agg_vec = ios_create_vec(m+n);
      mir->subst = xcalloc(1+m+n, sizeof(char));
      mir->mod_vec = ios_create_vec(m+n);
      mir->cut_vec = ios_create_vec(m+n);
      /* set global row attributes */
      set_row_attrib(tree, mir);
      /* set global column attributes */
      set_col_attrib(tree, mir);
      /* set variable bounds */
      set_var_bounds(tree, mir);
      /* mark rows which should not be used */
      mark_useless_rows(tree, mir);
      return mir;
}

/***********************************************************************
*  NAME
*
*  ios_mir_gen - generate MIR cuts
*
*  SYNOPSIS
*
*  #include "glpios.h"
*  void ios_mir_gen(glp_tree *tree, void *gen, IOSPOOL *pool);
*
*  DESCRIPTION
*
*  The routine ios_mir_gen generates MIR cuts for the current point and
*  adds them to the cut pool. */

static void get_current_point(glp_tree *tree, struct MIR *mir)
{     /* obtain current point */
      glp_prob *mip = tree->mip;
      int m = mir->m;
      int n = mir->n;
      int k;
      for (k = 1; k <= m; k++)
         mir->x[k] = mip->row[k]->prim;
      for (k = m+1; k <= m+n; k++)
         mir->x[k] = mip->col[k-m]->prim;
      return;
}

#if _MIR_DEBUG
static void check_current_point(struct MIR *mir)
{     /* check current point */
      int m = mir->m;
      int n = mir->n;
      int k, kk;
      double lb, ub, eps;
      for (k = 1; k <= m+n; k++)
      {  /* determine lower bound */