Adding a stability failure limit to simplex.

timothy-king · Feb 1, 2014 · a35b8ed · a35b8ed
1 parent 147e2f7
commit a35b8ed
Show file tree

Hide file tree

Showing 4 changed files with 24 additions and 1 deletion.
diff --git a/src/glpapi06.c b/src/glpapi06.c
@@ -488,6 +488,7 @@ void glp_init_smcp(glp_smcp *parm)
       parm->out_frq = 500;
       parm->out_dly = 0;
       parm->presolve = GLP_OFF;
+      parm->stability_lmt = 200;
       return;
 }
 

diff --git a/src/glpk.h b/src/glpk.h
@@ -130,6 +130,7 @@ typedef struct
       int out_frq;            /* spx.out_frq */
       int out_dly;            /* spx.out_dly (milliseconds) */
       int presolve;           /* enable/disable using LP presolver */
+      int stability_lmt;      /* maximum number of check stability failures before stopping */
       double foo_bar[36];     /* (reserved) */
 } glp_smcp;
 

diff --git a/src/glpspx01.c b/src/glpspx01.c
@@ -238,6 +238,9 @@ struct csa
       double *work2; /* double work2[1+m]; */
       double *work3; /* double work3[1+m]; */
       double *work4; /* double work4[1+m]; */
+
+      /** Things Tim has added. */
+      int stability_failures;
 };
 
 static const double kappa = 0.10;
@@ -400,6 +403,8 @@ static void init_csa(struct csa *csa, glp_prob *lp)
       csa->refct = 0;
       memset(&refsp[1], 0, (m+n) * sizeof(char));
       for (j = 1; j <= n; j++) gamma[j] = 1.0;
+
+      csa->stability_failures = 0;
       return;
 }
 
@@ -2647,6 +2652,11 @@ int spx_primal(glp_prob *lp, const glp_smcp *parm)
          if (check_stab(csa, parm->tol_bnd))
          {  /* there are excessive bound violations due to round-off
                errors */
+            csa->stability_failures++;
+            if (csa->stability_failures >= parm->stability_lmt){
+              ret = GLP_EINSTAB;
+              goto done;
+            }
             if (parm->msg_lev >= GLP_MSG_ERR)
                xprintf("Warning: numerical instability (primal simplex,"
                   " phase %s)\n", csa->phase == 1 ? "I" : "II");

diff --git a/src/glpspx02.c b/src/glpspx02.c
@@ -288,6 +288,9 @@ struct csa
       double *work2; /* double work2[1+m]; */
       double *work3; /* double work3[1+m]; */
       double *work4; /* double work4[1+m]; */
+
+      /* count the number of stability failures */
+      int stability_failures;
 };
 
 static const double kappa = 0.10;
@@ -501,6 +504,8 @@ static void init_csa(struct csa *csa, glp_prob *lp)
       csa->refct = 0;
       memset(&refsp[1], 0, (m+n) * sizeof(char));
       for (i = 1; i <= m; i++) gamma[i] = 1.0;
+
+      csa->stability_failures = 0;
       return;
 }
 
@@ -2741,7 +2746,13 @@ int spx_dual(glp_prob *lp, const glp_smcp *parm)
          /* make sure that the current basic solution remains dual
             feasible */
          if (check_stab(csa, parm->tol_dj) != 0)
-         {  if (parm->msg_lev >= GLP_MSG_ERR)
+         {
+            csa->stability_failures++;
+            if (csa->stability_failures >= parm->stability_lmt){
+              ret = GLP_EINSTAB;
+              goto done;
+            }
+            if (parm->msg_lev >= GLP_MSG_ERR)
                xprintf("Warning: numerical instability (dual simplex, p"
                   "hase %s)\n", csa->phase == 1 ? "I" : "II");
 #if 1