commit f056c1f9981e12c95a6057a277ab70c4262834d5
Author: Chris Matrakidis <cmatraki@gmail.com>
Date:   Mon Jan 9 14:05:12 2017 +0200

    dual simplex intrnal API for pcost

diff --git a/src/glpios09.c b/src/glpios09.c
index 727761c..1e2bde3 100644
--- a/src/glpios09.c
+++ b/src/glpios09.c
@@ -25,6 +25,7 @@
 #include "env.h"
 #include "glpios.h"
 #include "bfd.h"
+#include "simplex.h"
 
 /***********************************************************************
 *  NAME
@@ -414,6 +415,8 @@ struct csa
       /* copy of row bind values to avoid identical factorizations */
       int *c_bind;
       /* copy of column bind values to avoid identical factorizations */
+      void *duals;
+      /* internal state for dual simplex reuse */
 };
 
 void *ios_pcost_init(glp_tree *tree)
@@ -431,6 +434,7 @@ void *ios_pcost_init(glp_tree *tree)
       }
       csa->work = NULL;
       csa->bfd = NULL;
+      csa->duals = NULL;
       return csa;
 }
 
@@ -488,14 +492,18 @@ static double eval_degrad(glp_tree *T, int j, double bnd)
       }
       /* fix column x[j] at specified value */
       glp_set_col_bnds(lp, j, GLP_FX, bnd, bnd);
+      /* initialise solver state if empty or update bound */
+      if (csa->duals == NULL)
+         csa->duals = spy_dual_init(lp, 1);
+      else
+         spy_dual_update_col_bnd(lp, j, csa->duals);
       /* try to solve resulting LP */
       glp_init_smcp(&parm);
       parm.msg_lev = GLP_MSG_OFF;
       parm.meth = GLP_DUAL;
       parm.it_lim = 30;
       parm.out_dly = 1000;
-      parm.meth = GLP_DUAL;
-      ret = glp_simplex(lp, &parm);
+      ret = spy_dual_main(lp, &parm, csa->duals);
       if (ret == 0 || ret == GLP_EITLIM)
       {  if (glp_get_prim_stat(lp) == GLP_NOFEAS)
          {  /* resulting LP has no primal feasible solution */
@@ -531,6 +539,7 @@ static double eval_degrad(glp_tree *T, int j, double bnd)
       /* restore column x[j] type and bounds */
       glp_set_col_bnds(lp, j, P->col[j]->type, P->col[j]->lb,
          P->col[j]->ub);
+      spy_dual_update_col_bnd(lp, j, csa->duals);
       return degrad;
 }
 
@@ -721,6 +730,11 @@ done: *_next = sel;
          xfree(csa->r_bind);
          xfree(csa->head);
       }
+      if (csa->duals !=NULL)
+      {  /* clean up solver state for next time */
+         spy_dual_free(csa->duals);
+         csa->duals = NULL;
+      }
       return jjj;
 }
 
diff --git a/src/simplex/simplex.h b/src/simplex/simplex.h
index 4820588..3fa998f 100644
--- a/src/simplex/simplex.h
+++ b/src/simplex/simplex.h
@@ -34,6 +34,22 @@ int spx_primal(glp_prob *P, const glp_smcp *parm);
 int spy_dual(glp_prob *P, const glp_smcp *parm);
 /* driver to dual simplex method */
 
+#define spy_dual_init _glp_spy_dual_init
+struct csa *spy_dual_init(glp_prob *P, int delta);
+/* initialise dual simplex state to allow reuse */
+
+#define spy_dual_main _glp_spy_dual_main
+int spy_dual_main(glp_prob *P, const glp_smcp *parm, struct csa *csa);
+/* dual simplex using previously initialised state */
+
+#define spy_dual_free _glp_spy_dual_free
+void spy_dual_free(struct csa *csa);
+/* free dual simplex state */
+
+#define spy_dual_update_col_bnd _glp_spy_dual_update_col_bnd
+void spy_dual_update_col_bnd(glp_prob *P, int j, struct csa *csa);
+/* update column bounds in dual simplex state */
+
 #endif
 
 /* eof */
diff --git a/src/simplex/spydual.c b/src/simplex/spydual.c
index e41baf2..5e1c828 100644
--- a/src/simplex/spydual.c
+++ b/src/simplex/spydual.c
@@ -146,6 +146,12 @@ struct csa
 #endif
       int p_stat, d_stat;
       /* primal and dual solution statuses */
+      int *map; /* int map[1+P->m+P->n]; */
+      /* mapping of original LP variables to working LP ones */
+      int *daeh; /* int map[1+n]; */
+      /* array to hold inverse of the working basis header */
+      double *delta; /* double delta[1+P->n]; */
+      /* if SHIFT is set, holds shift amount (NULL if not used) */
       /*--------------------------------------------------------------*/
       /* control parameters (see struct glp_smcp) */
       int msg_lev;
@@ -1770,23 +1776,49 @@ fini:
 
 int spy_dual(glp_prob *P, const glp_smcp *parm)
 {     /* driver to dual simplex method */
-      struct csa csa_, *csa = &csa_;
-      SPXLP lp;
-#if USE_AT
-      SPXAT at;
-#else
-      SPXNT nt;
-#endif
-      SPYSE se;
-      int ret, *map, *daeh;
-      /* build working LP and its initial basis */
+      struct csa *csa;
+      int ret;
+      csa = spy_dual_init(P, 0);
+      ret = spy_dual_main(P, parm, csa);
+      spy_dual_free(csa);
+      /* return to calling program */
+      return ret;
+}
+
+/***********************************************************************
+*  spy_dual_init - initialise dual simplex state to allow reuse
+*
+*  This routine initialises the dual simplex state and returns a pointer
+*  to the csa structure that can be used in subsequent calls to dual
+*  simplex functions. */
+
+static double * spy_init_delta(glp_prob *P, int *map)
+{     int m = P->m;
+      int n = P->n;
+      int k;
+      double *delta = talloc(1+m+n, double);
+      for (k = 1; k <= n; k++)
+      {  GLPCOL *col = P->col[k];
+         if (map[m+k] == 0 || col->type == GLP_FR)
+            delta[k] = 0;
+         else if (map[m+k] > 0)
+            delta[k] = col->lb / col->sjj;
+         else
+            delta[k] = col->ub / col->sjj;
+      }
+      return delta;
+}
+
+struct csa *spy_dual_init(glp_prob *P, int use_delta)
+{     /* initialise dual simplex state to allow reuse */
+      struct csa *csa;
+      csa = talloc(1, struct csa);
       memset(csa, 0, sizeof(struct csa));
-      csa->lp = &lp;
+      csa->lp = talloc(1, SPXLP);
       spx_init_lp(csa->lp, P, EXCL);
       spx_alloc_lp(csa->lp);
-      map = talloc(1+P->m+P->n, int);
-      spx_build_lp(csa->lp, P, EXCL, SHIFT, map);
-      spx_build_basis(csa->lp, P, map);
+      csa->map = talloc(1+P->m+P->n, int);
+      spx_build_lp(csa->lp, P, EXCL, SHIFT, csa->map);
       switch (P->dir)
       {  case GLP_MIN:
             csa->dir = +1;
@@ -1798,46 +1830,13 @@ int spy_dual(glp_prob *P, const glp_smcp *parm)
             xassert(P != P);
       }
       csa->orig_b = talloc(1+csa->lp->m, double);
-      memcpy(csa->orig_b, csa->lp->b, (1+csa->lp->m) * sizeof(double));
 #if PERTURB
       csa->orig_c = talloc(1+csa->lp->n, double);
-      memcpy(csa->orig_c, csa->lp->c, (1+csa->lp->n) * sizeof(double));
 #endif
       csa->orig_l = talloc(1+csa->lp->n, double);
-      memcpy(csa->orig_l, csa->lp->l, (1+csa->lp->n) * sizeof(double));
       csa->orig_u = talloc(1+csa->lp->n, double);
-      memcpy(csa->orig_u, csa->lp->u, (1+csa->lp->n) * sizeof(double));
-#if USE_AT
-      /* build matrix A in row-wise format */
-      csa->at = &at;
-      csa->nt = NULL;
-      spx_alloc_at(csa->lp, csa->at);
-      spx_build_at(csa->lp, csa->at);
-#else
-      /* build matrix N in row-wise format for initial basis */
-      csa->at = NULL;
-      csa->nt = &nt;
-      spx_alloc_nt(csa->lp, csa->nt);
-      spx_init_nt(csa->lp, csa->nt);
-      spx_build_nt(csa->lp, csa->nt);
-#endif
-      /* allocate and initialize working components */
-      csa->phase = 0;
       csa->beta = talloc(1+csa->lp->m, double);
-      csa->beta_st = 0;
       csa->d = talloc(1+csa->lp->n-csa->lp->m, double);
-      csa->d_st = 0;
-      switch (parm->pricing)
-      {  case GLP_PT_STD:
-            csa->se = NULL;
-            break;
-         case GLP_PT_PSE:
-            csa->se = &se;
-            spy_alloc_se(csa->lp, csa->se);
-            break;
-         default:
-            xassert(parm != parm);
-      }
 #if 0 /* 30/III-2016 */
       csa->list = talloc(1+csa->lp->m, int);
       csa->trow = talloc(1+csa->lp->n-csa->lp->m, double);
@@ -1856,6 +1855,66 @@ int spy_dual(glp_prob *P, const glp_smcp *parm)
       fvs_alloc_vec(&csa->wrow, csa->lp->n-csa->lp->m);
       fvs_alloc_vec(&csa->wcol, csa->lp->m);
 #endif
+      csa->daeh = talloc(1+csa->lp->n, int);
+      if (SHIFT && use_delta != 0)
+         csa->delta = spy_init_delta(P, csa->map);
+      return csa;
+}
+
+/***********************************************************************
+*  spy_dual_main - dual simplex using previously initialised state
+*
+*  This routine is the main body of the dual simplex method. The csa
+*  structure is assumed to have been initialised by a call to
+*  spy_dual_init, however repeated calls to spy_dual_main can be done
+*  reusing the same csa structure.
+*
+*  On exit the return values are as in spy_dual described earlier */
+
+int spy_dual_main(glp_prob *P, const glp_smcp *parm, struct csa *csa)
+{     /* dual simplex using previously initialised state */
+      int ret;
+      /* build working LP initial basis */
+      spx_build_basis(csa->lp, P, csa->map);
+      memcpy(csa->orig_b, csa->lp->b, (1+csa->lp->m) * sizeof(double));
+#if PERTURB
+      memcpy(csa->orig_c, csa->lp->c, (1+csa->lp->n) * sizeof(double));
+#endif
+      memcpy(csa->orig_l, csa->lp->l, (1+csa->lp->n) * sizeof(double));
+      memcpy(csa->orig_u, csa->lp->u, (1+csa->lp->n) * sizeof(double));
+#if USE_AT
+      /* build matrix A in row-wise format */
+      if (csa->at == NULL)
+      {  csa->at = talloc(1, SPXAT);
+         spx_alloc_at(csa->lp, csa->at);
+         spx_build_at(csa->lp, csa->at);
+      }
+#else
+      /* build matrix N in row-wise format for initial basis */
+      if (csa->nt == NULL)
+      {  csa->nt = talloc(1, SPXNT);
+         spx_alloc_nt(csa->lp, csa->nt);
+         spx_init_nt(csa->lp, csa->nt);
+      }
+      spx_build_nt(csa->lp, csa->nt);
+#endif
+      /* initialize working components */
+      csa->phase = 0;
+      csa->beta_st = 0;
+      csa->d_st = 0;
+      switch (parm->pricing)
+      {  case GLP_PT_STD:
+            break;
+         case GLP_PT_PSE:
+            if (csa->se == NULL)
+            {  csa->se = talloc(1, SPYSE);
+               spy_alloc_se(csa->lp, csa->se);
+            }
+            csa->se->valid = 0;
+            break;
+         default:
+            xassert(parm != parm);
+      }
       /* initialize control parameters */
       csa->msg_lev = parm->msg_lev;
       csa->dualp = (parm->meth == GLP_DUALP);
@@ -1876,7 +1935,8 @@ int spy_dual(glp_prob *P, const glp_smcp *parm)
          case GLP_RT_HAR:
             break;
          case GLP_RT_FLIP:
-            csa->bp = talloc(1+csa->lp->n-csa->lp->m, SPYBP);
+            if (csa->bp == NULL)
+               csa->bp = talloc(1+csa->lp->n-csa->lp->m, SPYBP);
             break;
          default:
             xassert(parm != parm);
@@ -1927,16 +1987,14 @@ int spy_dual(glp_prob *P, const glp_smcp *parm)
          goto skip;
       /* convert working LP basis to original LP basis and store it to
        * problem object */
-      daeh = talloc(1+csa->lp->n, int);
-      spx_store_basis(csa->lp, P, map, daeh);
+      spx_store_basis(csa->lp, P, csa->map, csa->daeh);
       /* compute simplex multipliers for final basic solution found by
        * the solver */
       spx_eval_pi(csa->lp, csa->work);
       /* convert working LP solution to original LP solution and store
        * it to problem object */
-      spx_store_sol(csa->lp, P, SHIFT, map, daeh, csa->beta, csa->work,
-         csa->d);
-      tfree(daeh);
+      spx_store_sol(csa->lp, P, SHIFT, csa->map, csa->daeh, csa->beta,
+         csa->work, csa->d);
       /* save simplex iteration count */
       P->it_cnt = csa->it_cnt;
       /* report auxiliary/structural variable causing unboundedness */
@@ -1949,16 +2007,26 @@ int spy_dual(glp_prob *P, const glp_smcp *parm)
          xassert(1 <= k && k <= csa->lp->n);
          /* convert to number of original variable */
          for (kk = 1; kk <= P->m + P->n; kk++)
-         {  if (abs(map[kk]) == k)
+         {  if (abs(csa->map[kk]) == k)
             {  P->some = kk;
                break;
             }
          }
          xassert(P->some != 0);
       }
-skip: /* deallocate working objects and arrays */
+skip: /* return to calling program */
+      return ret >= 0 ? ret : GLP_EFAIL;
+}
+
+/***********************************************************************
+*  spy_dual_free - free dual simplex state
+*
+*  This routine clears and deallocates the csa structure. */
+
+void spy_dual_free(struct csa *csa)
+{     /* free dual simplex state */
       spx_free_lp(csa->lp);
-      tfree(map);
+      tfree(csa->map);
       tfree(csa->orig_b);
 #if PERTURB
       tfree(csa->orig_c);
@@ -1995,8 +2063,117 @@ skip: /* deallocate working objects and arrays */
       fvs_free_vec(&csa->wrow);
       fvs_free_vec(&csa->wcol);
 #endif
-      /* return to calling program */
-      return ret >= 0 ? ret : GLP_EFAIL;
+      tfree(csa->daeh);
+      if (csa->delta != NULL)
+         tfree(csa->delta);
+      tfree(csa->lp);
+      if (csa->at != NULL)
+         tfree(csa->at);
+      if (csa->nt != NULL)
+         tfree(csa->nt);
+      if (csa->se != NULL)
+         tfree(csa->se);
+      tfree(csa);
+      return;
+}
+
+/***********************************************************************
+*  spy_dual_update_col_bnd - update column bounds in dual simplex state
+*
+*  This routine updates the csa structure for a changed bound in a
+*  structural variable of the original LP so that spy_dual_main can be
+*  called again without needing a new call to csa_dual_init. */
+
+void spy_dual_update_col_bnd(glp_prob *P, int j, struct csa *csa)
+{     /* update column bounds in dual simplex state */
+      SPXLP *lp = csa->lp;
+      int *map = csa->map;
+      int m = lp->m;
+      int n = lp->n;
+      double *l = lp->l;
+      double *u = lp->u;
+      int k;
+      double delta;
+      GLPCOL *col = P->col[j];
+      GLPAIJ *aij;
+      k = abs(map[m+j]);
+      xassert(k != 0);
+      /* update bounds of variable x[k] in working LP */
+      switch (col->type)
+      {  case GLP_FR:
+            l[k] = -DBL_MAX, u[k] = +DBL_MAX;
+            break;
+         case GLP_LO:
+            l[k] = col->lb / col->sjj, u[k] = +DBL_MAX;
+            break;
+         case GLP_UP:
+            l[k] = -DBL_MAX, u[k] = col->ub / col->sjj;
+            break;
+         case GLP_DB:
+            l[k] = col->lb / col->sjj, u[k] = col->ub / col->sjj;
+            xassert(l[k] != u[k]);
+            break;
+         case GLP_FX:
+            l[k] = u[k] = col->lb / col->sjj;
+            break;
+         default:
+            xassert(col != col);
+      }
+      if (SHIFT)
+      {  /* shift bounds of variable x[k] */
+         if (l[k] == -DBL_MAX && u[k] == +DBL_MAX)
+         {  /* x[k] is unbounded variable */
+            delta = 0.0;
+         }
+         else if (l[k] != -DBL_MAX && u[k] == +DBL_MAX)
+         {  /* shift lower bound to zero */
+            delta = l[k];
+            l[k] = 0.0;
+         }
+         else if (l[k] == -DBL_MAX && u[k] != +DBL_MAX)
+         {  /* shift upper bound to zero */
+            map[m+j] = -k;
+            delta = u[k];
+            u[k] = 0.0;
+         }
+         else if (l[k] != u[k])
+         {  /* x[k] is double bounded variable */
+            if (fabs(l[k]) <= fabs(u[k]))
+            {  /* shift lower bound to zero */
+               delta = l[k];
+               l[k] = 0.0, u[k] -= delta;
+            }
+            else
+            {  /* shift upper bound to zero */
+               map[m+j] = -k;
+               delta = u[k];
+               l[k] -= delta, u[k] = 0.0;
+            }
+            xassert(l[k] != u[k]);
+         }
+         else
+         {  /* shift fixed value to zero */
+            delta = l[k];
+            l[k] = u[k] = 0.0;
+         }
+         /* update all constraints and the objective function of
+          * working LP */
+         if (delta != csa->delta[j])
+         {  int ptr, end;
+            int *A_ptr = lp->A_ptr;
+            int *A_ind = lp->A_ind;
+            double *A_val = lp->A_val;
+            double *b = lp->b;
+            double *c = lp->c;
+            ptr = A_ptr[k];
+            end = A_ptr[k+1];
+            for (; ptr < end; ptr++)
+               b[A_ind[ptr]] -= A_val[ptr] * (delta - csa->delta[j]);
+            c[0] += c[k] * (delta - csa->delta[j]);
+            csa->delta[j] = delta;
+         }
+      }
+      return;
 }
 
 /* eof */