[Fwd: Bug in glpk integer programmer]

[Andrew Makhorin]

-------- Forwarded Message --------
From: s meagher <sjmeagher@gmail.com>
To: bug-glpk@gnu.org
Subject: Bug in glpk integer programmer
Date: Fri, 16 Oct 2020 00:19:14 +1100

> Hi
> 
> The following code demonstrates a possible bug in glpk. I fully
> appreciate that in fact the bug may be that I do understand glpk
> properly and have not configured it correctly. If this is the case,
> please accept my apologies for taking up your time.
> 
> Kind regards
> 
> Steve Meagher
> 
> /*    Possible integer programming bug in glpk v4.65
> 
>       (c) Stephen Meagher 2020
> 
>       The following code was written to demonstrate a bug I believed I
> found for integer programming in glpk. Instead, I seem to have
> discovered a different issue - which may be just because I'm using
> glpk incorrectly.
>       
>       The following integer programme
>       
>               objective function: f(x_1,x_2) = 0
>               row 1                           x_1 + x_2 = 3
>               row 2                   0<=     x_1 - x_2 <= 2
>               
> clearly has feasible solution x1 = 2, x2 = 1, but if I set it up as
>       
>               objective function: f(x_1,x_2) = 0
>               row 1                           x_1 + x_2 = 3
>               row 2                   0<=     x_1 - x_2 <=
> 2*(1.05)
>               
> then glpk tells me the solution is x1 = 2.55 and x2 = 0.45, which is
> wrong for two reasons. (Note, even if I set it up with 2 instead of
> 2*(1.05), glpk gives the wrong answer when configured as follows).
>       
>       I have compiled this programme on a MacBook Pro running macOS
> Catalina 10.15.6 using the command: 
>       
>       gcc -std=c18 glpkipbug.c -lglpk
>       
>       I have the gnu big num library installed and use it to configure
> glpk when I installed it.
>       
>       Here is the fulloutput:
>               
>       GLPK Simplex Optimizer, v4.65
> 2 rows, 2 columns, 4 non-zeros
>       0: obj =   0.000000000e+00 inf =   3.000e+00 (1)
>       2: obj =   0.000000000e+00 inf =   0.000e+00 (0)
> OPTIMAL LP SOLUTION FOUND
> GLPK Integer Optimizer, v4.65
> 2 rows, 2 columns, 4 non-zeros
> 2 integer variables, none of which are binary
> Preprocessing...
> 2 rows, 2 columns, 4 non-zeros
> 2 integer variables, none of which are binary
> Scaling...
>  A: min|aij| =  1.000e+00  max|aij| =  1.000e+00  ratio =  1.000e+00
> Problem data seem to be well scaled
> Constructing initial basis...
> Size of triangular part is 2
> Solving LP relaxation...
> GLPK Simplex Optimizer, v4.65
> 2 rows, 2 columns, 4 non-zeros
>       2: obj =   0.000000000e+00 inf =   3.000e+00 (1)
>       3: obj =   0.000000000e+00 inf =   0.000e+00 (0)
> OPTIMAL LP SOLUTION FOUND
> Integer optimization begins...
> Long-step dual simplex will be used
> +     3: mip =     not found yet >=              -inf        (1; 0)
> +     4: >>>>>   0.000000000e+00 >=   0.000000000e+00   0.0% (1; 0)
> +     4: mip =   0.000000000e+00 >=     tree is empty   0.0% (0; 1)
> INTEGER OPTIMAL SOLUTION FOUND
> Intopt success
> Solution as x1 = 2.550000, x2 = 0.450000
>               
> 
> (If you're still reading: The original bug occurred when I set certain
> row bounds to be inequalities. Then the rows which were equalities
> were no longer satisfied. However, this bug disappeared if I set the
> bounds to be integers instead of real numbers. Furthermore, I wrote
> that programme in C++ and not C.)
> 
> */
> 
> 
> #include <glpk.h>
> #include <stdio.h>
> 
> /* Defines and demontrates the bug */
> void bug() {
>       glp_prob *lp;
>       glp_iocp parm;
>       int ia[5];
>       int ja[5];
>       double arr[5];
>       lp = glp_create_prob();
>       glp_set_obj_dir(lp,GLP_MIN);
>       glp_init_iocp(&parm);
>       parm.presolve = GLP_ON;
>       glp_add_rows(lp,2);
>       glp_set_row_name(lp,1,"r1");
>       glp_set_row_bnds(lp,1,GLP_FX,3.0,3.0);
>       glp_set_row_name(lp,2,"r2");
>               
>       glp_set_row_bnds(lp,2,GLP_DB, 0,2*(1.05));      // non-integer 
> bounds
>       
>       glp_add_cols(lp,2);
>       glp_set_col_name(lp,1,"x1");
>       glp_set_obj_coef(lp,1,0.0);
>       glp_set_col_bnds(lp,1,GLP_LO,0.0,0.0);
>       glp_set_col_kind(lp,1,GLP_IV);
>       glp_set_col_name(lp,2,"x2");
>       glp_set_obj_coef(lp,2,0.0);
>       glp_set_col_bnds(lp,2,GLP_LO,0.0,0.0);
>       glp_set_col_kind(lp,2,GLP_IV);
>       ia[1] = 1; ja[1] = 1; arr[1] = 1; ia[2] = 1; ja[2] = 2; arr[2] =
> 1;
>       ia[3] = 2; ja[3] = 1; arr[3] = 1; ia[4] = 2; ja[4] = 2; arr[4] =
> -1;
>       glp_load_matrix(lp,4,ia,ja,arr);
>       glp_simplex(lp,NULL);
>       if (glp_intopt(lp, &parm) != 0) {
>               printf("Intopt error\n");       
>       } else {
>               printf("Intopt success\n");     
>       }
>       double x1 = glp_get_col_prim(lp,1);
>       double x2 = glp_get_col_prim(lp,2);
>       printf("Solution as x1 = %lf, x2 = %lf\n",x1,x2);
>       glp_delete_prob(lp);
> }
> 
> int main() {
>       bug();
>       return 0;
> }
> 
> 
>