[Bug-glpk] [PATCH] Fix pivot column selection in rational simplex.

[Gabriele Farina]

This patch fixes a glitch in the pivot column selection routine of
the rational simplex (`ssx_chuzc`).

Currently, the pivot column selection works by finding the column that
has the largest reduced cost. The reduced costs are compared
*as floating-point numbers*. The glitch occurs when a reduced cost
is so small that it is converted to the floating-point number `0.0`,

which causes an assertion to be thrown.

In order to not significantly degrade performance, this patch adopts
a hybrid approach: first, a comparison of reduced costs is attempted
using floating-point numbers. If the floating-point number is zero,
the reduced costs are compared as rational numbers instead.

I could not figure out where I should put a test case to make sure we

don't regress. Can you please point me in the right direction?

Thanks,

Gabriele

---
 src/draft/glpssx01.c | 21 +++++++++++++++++----
 1 file changed, 17 insertions(+), 4 deletions(-)

diff --git a/src/draft/glpssx01.c b/src/draft/glpssx01.c
index 9b70444..063671c 100644
--- a/src/draft/glpssx01.c
+++ b/src/draft/glpssx01.c
@@ -461,10 +461,14 @@ void ssx_chuzc(SSX *ssx)
       int *Q_col = ssx->Q_col;
       int *stat = ssx->stat;
       mpq_t *cbar = ssx->cbar;
+
+      mpq_t abs_best; mpq_init(abs_best);
+      mpq_t abs_cur; mpq_init(abs_cur);
+
       int j, k, s, q, q_dir;
       double best, temp;
       /* nothing is chosen so far */
-      q = 0, q_dir = 0, best = 0.0;
+      q = 0, q_dir = 0;
       /* look through the list of non-basic variables */
       for (j = 1; j <= n; j++)
       {  k = Q_col[m+j]; /* x[k] = xN[j] */
@@ -474,12 +478,21 @@ void ssx_chuzc(SSX *ssx)
          {  /* reduced cost of xN[j] indicates possible improving of
                the objective function */
             temp = fabs(mpq_get_d(cbar[j]));
-            xassert(temp != 0.0);
-            if (q == 0 || best < temp)
-               q = j, q_dir = - s, best = temp;
+            xassert(temp >= 0.0 && best >= 0.0);
+            if (temp == 0.0 && (q > 0 && best == 0.0)) {
+              mpq_abs(abs_cur, cbar[j]);
+              mpq_abs(abs_best, cbar[q]);
+              if (q == 0 || mpq_cmp(abs_best, abs_cur) < 0)
+                 q = j, q_dir = - s, best = temp;
+            } else {
+               if (q == 0 || best < temp)
+                 q = j, q_dir = - s, best = temp;
+            }
          }
       }
       ssx->q = q, ssx->q_dir = q_dir;
+      mpq_clear(abs_best);
+      mpq_clear(abs_cur);
       return;
 }

--
2.17.1