aboutsummaryrefslogtreecommitdiffstats
path: root/test/monniaux/glpk-4.65/src/intopt/gmigen.c
blob: 627682cbd3741f53201993492c2e92b8045ecef0 (plain)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
/* gmigen.c (Gomory's mixed integer cuts generator) */

/***********************************************************************
*  This code is part of GLPK (GNU Linear Programming Kit).
*
*  Copyright (C) 2002-2018 Andrew Makhorin, Department for Applied
*  Informatics, Moscow Aviation Institute, Moscow, Russia. All rights
*  reserved. E-mail: <mao@gnu.org>.
*
*  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 "env.h"
#include "prob.h"

/***********************************************************************
*  NAME
*
*  glp_gmi_gen - generate Gomory's mixed integer cuts
*
*  SYNOPSIS
*
*  int glp_gmi_gen(glp_prob *P, glp_prob *pool, int max_cuts);
*
*  DESCRIPTION
*
*  This routine attempts to generate Gomory's mixed integer cuts for
*  integer variables, whose primal values in current basic solution are
*  integer infeasible (fractional).
*
*  On entry to the routine the basic solution contained in the problem
*  object P should be optimal, and the basis factorization should be
*  valid.
*
*  The cutting plane inequalities generated by the routine are added to
*  the specified cut pool.
*
*  The parameter max_cuts specifies the maximal number of cuts to be
*  generated. Note that the number of cuts cannot exceed the number of
*  basic variables, which is the number of rows in the problem object.
*
*  RETURNS
*
*  The routine returns the number of cuts that have been generated and
*  added to the cut pool. */

#define f(x) ((x) - floor(x))
/* compute fractional part of x */

struct var { int j; double f; };

static int CDECL fcmp(const void *p1, const void *p2)
{     const struct var *v1 = p1, *v2 = p2;
      if (v1->f > v2->f) return -1;
      if (v1->f < v2->f) return +1;
      return 0;
}

int glp_gmi_gen(glp_prob *P, glp_prob *pool, int max_cuts)
{     int m = P->m;
      int n = P->n;
      GLPCOL *col;
      struct var *var;
      int i, j, k, t, len, nv, nnn, *ind;
      double frac, *val, *phi;
      /* sanity checks */
      if (!(P->m == 0 || P->valid))
         xerror("glp_gmi_gen: basis factorization does not exist\n");
      if (!(P->pbs_stat == GLP_FEAS && P->dbs_stat == GLP_FEAS))
         xerror("glp_gmi_gen: optimal basic solution required\n");
      if (pool->n != n)
         xerror("glp_gmi_gen: cut pool has wrong number of columns\n");
      /* allocate working arrays */
      var = xcalloc(1+n, sizeof(struct var));
      ind = xcalloc(1+n, sizeof(int));
      val = xcalloc(1+n, sizeof(double));
      phi = xcalloc(1+m+n, sizeof(double));
      /* build the list of integer structural variables, which are
       * basic and have integer infeasible (fractional) primal values
       * in optimal solution to specified LP */
      nv = 0;
      for (j = 1; j <= n; j++)
      {  col = P->col[j];
         if (col->kind != GLP_IV)
            continue;
         if (col->type == GLP_FX)
            continue;
         if (col->stat != GLP_BS)
            continue;
         frac = f(col->prim);
         if (!(0.05 <= frac && frac <= 0.95))
            continue;
         /* add variable to the list */
         nv++, var[nv].j = j, var[nv].f = frac;
      }
      /* sort the list by descending fractionality */
      qsort(&var[1], nv, sizeof(struct var), fcmp);
      /* try to generate cuts by one for each variable in the list, but
       * not more than max_cuts cuts */
      nnn = 0;
      for (t = 1; t <= nv; t++)
      {  len = glp_gmi_cut(P, var[t].j, ind, val, phi);
         if (len < 1)
            goto skip;
         /* if the cut inequality seems to be badly scaled, reject it
          * to avoid numerical difficulties */
         for (k = 1; k <= len; k++)
         {  if (fabs(val[k]) < 1e-03)
               goto skip;
            if (fabs(val[k]) > 1e+03)
               goto skip;
         }
         /* add the cut to the cut pool for further consideration */
         i = glp_add_rows(pool, 1);
         glp_set_row_bnds(pool, i, GLP_LO, val[0], 0);
         glp_set_mat_row(pool, i, len, ind, val);
         /* one cut has been generated */
         nnn++;
         if (nnn == max_cuts)
            break;
skip:    ;
      }
      /* free working arrays */
      xfree(var);
      xfree(ind);
      xfree(val);
      xfree(phi);
      return nnn;
}

/* eof */