/*
  This file is a part of the Dylp LP distribution.

        Copyright (C) 2005 -- 2007 Lou Hafer

        School of Computing Science
        Simon Fraser University
        Burnaby, B.C., V5A 1S6, Canada
        lou@cs.sfu.ca

  This code is licensed under the terms of the Eclipse Public License (EPL).
*/

/*
  This file contains the top-level routines for the dylp dynamic simplex
  package.

  The basic flow of the dynamic simplex algorithm is as follows:
   
    0) Initialisation: Initialise the constraint system with equalities and
       some suitable selection of inequalities. Go to 4.
    1) Purge Variables: Based on some threshold for reduced cost, deactivate
       variables which are unlikely to improve the objective.
    2) Generate Variables: Generate new variables not part of the original
       problem statement (column generation). These are added to the set of
       inactive variables.
    3) Add Variables: Price out the inactive variables. If any of them have
       suitable reduced cost, add them to the set of active variables and go
       to 4, otherwise go to 5.
    4) Primal: Optimise with primal simplex. Go to 1.
    5) Purge Constraints: If the objective strictly improved (indicating that
       we've moved to a new extreme point), deactivate any slack constraints.
    6) Generate Constraints: Generate new constraints, not part of the original
       problem statement (cut generation). Equalities go directly into the
       active set, inequalities go into the inactive set.
    7) Add Constraints: Check the inactive constraints for violations. If any
       violated constraints are found, add them to the active constraints.
       If no violated constraints are found, we're optimal over all constraints
       and variables.
    8) Dual: Optimise with dual simplex. Go to 1. If the user has forced use
       of primal phase I in place of dual phase II, a silent punt occurs
       immediately.
  
  If the caller specifies the fullsys option, the full constraint system is
  loaded and all constraint and variable deactivation/activation is skipped.

  The dual simplex is somewhat weaker than the primal simplex -- it has no
  phase I, and has only antidegen lite. Consequently, the dual will punt to
  the primal phase I when it finds itself in trouble due to loss of
  feasibility or stalling. It also punts to primal phase I when it needs to
  add variables but can't find any that are dual feasible.

  Steps 2 and 6 are implemented as stubs, as they are highly problem specific.
  dylp requires that any variables or constraints generated by a call must
  be inserted in the original constraint system. They can be left to be picked
  up in steps 3 or 7, or the generation routines can force them into the active
  sets as they are generated.

  Inactive constraints are assumed to be loose. They're pretty simple to
  activate -- the row is added to dy_sys with a logical, which is made
  basic.  If the constraint contains other variables, beyond the currrent
  active set, they remain inactive (this avoids any potential trouble in
  terms of maintaining dual feasibility).  The only point to be made is that
  for maximum efficiency all constraints in dy_sys are billed as
  architectural constraints. (This prevents consys from working harder to
  keep the architectural and cut constraints as two distinct sets.)

  Inactive variables require a bit more bookkeeping when they're inactive.
  Since we have to be able to construct a basic feasible (i.e., extreme
  point) solution, inactive variables must be at bound. Texts invariably
  assume the simple case of a lower bound of 0 and no upper bound; in this
  case no special action is needed. When a variable x<j> has non-zero lower
  and/or upper bounds, we'll need to correct b<i> in dy_sys whenever a<i,j>
  is non-zero, since the simplex routines have no other way of seeing the
  effect of the bound.

  It's worth mentioning again that dylp and consys use a nonstandard convention
  for variable indices. For a constraint system with m constraints and n
  architectural variables, logicals occupy indices 1..m and architectural
  variables occupy indices m+1..m+n.
*/

#define DYLP_INTERNAL

#include "dylp.h"

static char sccsid[] UNUSED = "@(#)dylp.c	4.7	10/15/05" ;
static char svnid[] UNUSED = "$Id: dylp.c 407 2010-12-31 20:48:48Z lou $" ;

/*
  To avoid passing an excessive number of parameters around inside the dylp
  implementation, the following globals are used to communicate between modules.
  These are defined in dylp.h, but their use is limited by the DYLP_INTERNAL
  conditional compilation variable.
*/

lp_struct *dy_lp = NULL ;
consys_struct *dy_sys = NULL ;
lpopts_struct *dy_opts = NULL ;
lptols_struct *dy_tols = NULL ;
#ifdef DYLP_STATISTICS
lpstats_struct *dy_stats = NULL ;
#endif
int *dy_actvars = NULL,
    *dy_actcons = NULL,
    *dy_origvars = NULL,
    *dy_origcons = NULL,
    *dy_basis = NULL,
    *dy_var2basis = NULL,
    *dy_brkout = NULL,
    *dy_degenset = NULL,
    *dy_ddegenset = NULL ;
flags *dy_status = NULL ;
double *dy_x = NULL,
       *dy_xbasic = NULL,
       *dy_y = NULL,
       *dy_gamma = NULL,
       *dy_cbar = NULL,
       *dy_rho = NULL ;
bool *dy_frame = NULL ;

/*
  In an object-oriented environment, there may well be multiple client
  objects in existence. For efficiency, they should be able to use hot
  starts, but that requires the ability to track who `owns' the solver,
  in the sense that the problem loaded into the data structures (and basis
  factorisation) is the problem belonging to that client. Hence this
  variable:

  dy_owner	ID of the client whose problem is loaded.

  Ownership is checked in dylp as part of the startup sequence and set in
  dy_finishup. Various other routines (in particular, the ones that produce
  solution vectors and tableau variables in the external reference frame)
  check as needed.

  Beware! The basic principle is that dylp checks only as far as it needs to
  in order to protect itself. It will not prevent one client from forcibly
  evicting another!

  It's up to the client to check whether it owns the solver and invoke
  the actions necessary to gain ownership: specifically, unload the other
  client's problem (or ask the other client to unload) and then load in
  its own. This implies that clients in a hostile environment should keep
  enough information to warm start (or be willing to do a cold start),
  as some other client could grab the solver between any two calls by a
  single client.

  A call with context cxUNLOAD will unload retained structures. A call with
  context cxLOAD or cxUSERPIV will do either a cold or warm start, just to
  the point where dylp is ready to begin (or resume) pivoting.

  In object-oriented environments, a handy ID is the value of `this'
  for the client object (for example, so that one client can ask another
  client to give up ownership).
*/

void *dy_owner = NULL ;

/*
  For any of the above to be useful, a client must be able to check ownership.
*/

void *dy_getOwner () { return (dy_owner) ; }

/*
  Startup control type.
*/

typedef enum {startCOLD, startWARM, startHOT} start_enum ;



static void updateOptsAndTols (lpopts_struct *client_opts,
			       lptols_struct *client_tols)
/*
  A nearly trivial routine to update options and tolerances. Pulled out of
  dylp startup sequence to keep down the clutter from comments and provide a
  single point of change. Most of the options and tolerances are read-only,
  but dylp already adjusts a few on the fly and this is likely to grow.

  The general strategy is simple enough: If we have retained structures,
  create local copies of each structure on the stack and use them to hold the
  current settings.  Copy the client's options and tolerances into the dylp
  structures, then bring back the fields that dylp modifies. The bulk copy is
  overkill but robust against changes in the structures.

  Options modified by dylp:
    dpsel.strat

  Tolerances modified by dylp:
    pfeas
    dfeas

  Parameters:
    client_opts:	options supplied by client
    client_tols:	tolerances supplied by client

  Returns: undefined
*/

{ lptols_struct lcl_tols ;
  lpopts_struct lcl_opts ;

# if defined(DYLP_PARANOIA) || MALLOC_DEBUG == 2
  const char *rtnnme = "updateOptsAndTols" ;
# endif
  
/*
  Allocate dylp's structures, if they don't exist, or make copies if they do.
  It should be the case that we have structures iff dy_owner is non-null.
*/
# ifdef DYLP_PARANOIA
  if ((dy_owner != NULL && (dy_tols == NULL || dy_opts == NULL)) ||
      (dy_owner == NULL && (dy_tols != NULL || dy_opts != NULL)))
  { errmsg(1,rtnnme,__LINE__) ;
    return ; }
# endif

  if (dy_owner != NULL)
  { memcpy(&lcl_tols,dy_tols,sizeof(lptols_struct)) ;
    memcpy(&lcl_opts,dy_opts,sizeof(lpopts_struct)) ; }
  else
  { dy_tols = (lptols_struct *) MALLOC(sizeof(lptols_struct)) ;
    dy_opts = (lpopts_struct *) MALLOC(sizeof(lpopts_struct)) ; }
/*
  Copy the client's structures into dylp's structures.
*/
  memcpy(dy_tols,client_tols,sizeof(lptols_struct)) ;
  memcpy(dy_opts,client_opts,sizeof(lpopts_struct)) ;
/*
  If we had retained structures, copy over the fields that dylp modifies. If
  we're starting afresh, make sure we have a sane default in place for pfeas
  and dfeas.
*/
  if (dy_owner != NULL)
  { dy_opts->dpsel.strat = lcl_opts.dpsel.strat ;
    dy_tols->pfeas = lcl_tols.pfeas ;
    dy_tols->dfeas = lcl_tols.dfeas ; }
  else
  { if (dy_tols->pfeas <= 0)
    { dy_tols->pfeas = dy_tols->pfeas_scale*dy_tols->zero ; }
    if (dy_tols->dfeas <= 0)
    { dy_tols->dfeas = dy_tols->dfeas_scale*dy_tols->cost ; } }

  return ; }



static dyphase_enum addcon_nextphase (int actcnt)

/*
  This routine determines the appropriate next state after the addition of
  violated constraints in phase dyADDCON.  It's just a bit too complicated to
  leave laying out in the open. We look at the active (just finished) simplex
  and its return code, the target simplex, and the number of constraints that
  were activated (the parameter actcnt).

  Parameter:
    actcnt:	the number of constraints just activated; may be negative if
		there's been an error

  Returns: appropriate next phase, one of dyPRIMAL1, dyDUAL, dyDONE, or dyINV.
*/

{ dyphase_enum retval = dyINV ;
  flags chkflgs = ladPRIMFEAS|ladPFQUIET ;
  const char *rtnnme = "addcon_nextphase" ;

# ifndef NDEBUG
  if (dy_opts->print.major >= 2)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
	"\n    Entering simplex %s (%s), loadable %d, activated %d.",
	dy_prtlpphase(dy_lp->simplex.active,TRUE),dy_prtlpret(dy_lp->lpret),
	dy_lp->sys.cons.loadable,actcnt) ; }
# endif
	      
  if (actcnt < 0) return (dyINV) ;

  switch (dy_lp->simplex.active)
  { case dyPRIMAL2:
    { switch (dy_lp->lpret)
      { case lpOPTIMAL:
	{ if (actcnt == 0)
	  { retval = dyDONE ; }
	  else
	  { retval = dyPURGECON ; }
	  break ; }
	case lpUNBOUNDED:
	{ if (actcnt == 0)
	  { if (dy_lp->sys.cons.loadable > 0)
	    { retval = dyFORCEFULL ; }
	    else
	    { retval = dyDONE ; } }
	  else
	  { retval = dy_lp->simplex.next ; }
	  break ; }
	case lpSWING:
	{ if (actcnt == 0)
	  { retval = dyPRIMAL2 ;
	    dy_lp->simplex.next = dyPRIMAL2 ; }
	  else
	  { retval = dy_lp->simplex.next ; }
	  break ; }
        case lpFORCEDUAL:
	{ if (actcnt == 0)
	  { retval = dyFORCEFULL ; }
	  else
	  { retval = dyPURGECON ; }
	  break ; }
	default:
	{ break ; } }
      break ; }
    case dyPRIMAL1:
    { switch (dy_lp->lpret)
      { case lpUNBOUNDED:
	{ if (actcnt == 0)
	  { retval = dyFORCEFULL ; }
	  else
	  { retval = dyPRIMAL1 ;
	    dy_lp->simplex.next = dyPRIMAL1 ; }
	  break ; }
	case lpSWING:
	{ if (actcnt == 0)
	  { retval = dyPRIMAL1 ;
	    dy_lp->simplex.next = dyPRIMAL1 ; }
	  else
	  { retval = dy_lp->simplex.next ; }
	  break ; }
        case lpFORCEDUAL:
	{ if (actcnt == 0)
	  { retval = dyFORCEFULL ; }
	  else
	  { retval = dyPURGECON ; }
	  break ; }
        case lpPUNT:
        case lpSTALLED:
	{ if (actcnt == 0)
	  { retval = dyFORCEDUAL ; }
	  else
	  { retval = dyPRIMAL1 ;
	    dy_lp->simplex.next = dyPRIMAL1 ; }
	  break ; }
	default:
	{ break ; } }
      break ; }
    case dyDUAL:
    { switch (dy_lp->lpret)
      { case lpOPTIMAL:
	{ if (actcnt == 0)
	  { retval = dyDONE ; }
	  else
	  { retval = dyGENVAR ; }
	  break ; }
	case lpSWING:
	{ if (actcnt == 0)
	  { if (dy_lp->sys.vars.loadable <= 0)
	    { retval = dyPRIMAL1 ;
	      dy_lp->simplex.next = dyPRIMAL1 ; }
	    else
	    { retval = dyADDVAR ; } }
	  else
	  { retval = dy_lp->simplex.next ; }
	  break ; }
	case lpPUNT:
	case lpSTALLED:
	{ if (actcnt == 0)
	  { retval = dyFORCEPRIMAL ; }
	  else
	  { retval = dyDUAL ;
	    dy_lp->simplex.next = dyDUAL ; }
	  break ; }
	default:
	{ break ; } }
      break ; }
    default:
    { break ; } }

  if (retval == dyDUAL && dy_opts->usedual == FALSE)
  { retval = dyPRIMAL1 ;
    dy_lp->simplex.next = dyPRIMAL1 ; }

/*
  The call to dy_accchk is needed to ensure that the primal infeasibility
  information is consistent before primal1 calls initp1obj.
*/
  if (dy_lp->simplex.next == dyPRIMAL1)
  { if (dy_accchk(&chkflgs) != dyrOK)
    { retval = dyINV ; } }

  if (retval == dyINV)
  { errmsg(435,rtnnme,
	   dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,
	   actcnt,"constraints",dy_prtlpphase(dy_lp->simplex.active,TRUE),
	   dy_prtlpret(dy_lp->lpret),
	   dy_prtlpphase(dy_lp->simplex.next,TRUE)) ; }
  
# ifndef NDEBUG
  if (dy_opts->print.major >= 2)
  { dyio_outfmt(dy_logchn,dy_gtxecho,"\n    Leaving phase %s, simplex %s.",
		dy_prtlpphase(retval,TRUE),
		dy_prtlpphase(dy_lp->simplex.next,TRUE)) ; }
# endif

  return (retval) ; }



static dyphase_enum addvar_nextphase (int actcnt)

/*
  This routine determines the appropriate next state after the appropriate
  variable activation routine has been called in state dyADDVAR. It's just a
  bit too complicated to leave laying out in the open. We look at the active
  (just finished) simplex and its return code, the target simplex, and the
  number of variables that were activated (the parameter actcnt).

  Parameter:
    actcnt:	the number of variables just activated (may be negative if
		there's been an error

  Returns: appropriate next phase; one of dyPRIMAL1, dyPRIMAL2, dyDUAL,
	   dyPURGECON, dyFORCEDUAL, dyDONE, or dyINV.
*/
{ dyphase_enum retval = dyINV ;
  flags chkflgs = ladPRIMFEAS|ladPFQUIET ;
  const char *rtnnme = "addvar_nextphase" ;

  if (actcnt < 0) return (dyINV) ;

  switch (dy_lp->simplex.active)
  { case dyPRIMAL1:
    { switch (dy_lp->lpret)
      { case lpINFEAS:
	{ if (actcnt == 0)
	  { retval = dyDONE ; }
	  else
	  { retval = dyPRIMAL1 ;
	    dy_lp->simplex.next = dyPRIMAL1 ; }
	  break ; }
        case lpPUNT:
        case lpSTALLED:
	{ if (actcnt == 0)
	  { if (dy_lp->sys.cons.loadable > 0)
	    { retval = dyGENCON ; }
	    else
	    { retval = dyFORCEDUAL ; } }
	  else
	  { retval = dyPRIMAL1 ;
	    dy_lp->simplex.next = dyPRIMAL1 ; }
	  break ; }
	case lpFORCEDUAL:
	{ retval = dyDUAL ;
	  dy_lp->simplex.next = dyDUAL ;
	  break ; }
	default:
	{ break ; } }
      break ; }
    case dyPRIMAL2:
    { switch (dy_lp->lpret)
      { case lpOPTIMAL:
	{ if (dy_lp->simplex.next == dyPRIMAL2)
	  { if (actcnt == 0)
	    { retval = dyGENCON ;
	      dy_lp->simplex.next = dyDUAL ; }
	    else
	    { retval = dyPRIMAL2 ; } }
	  else
	  if (dy_lp->simplex.next == dyDUAL)
	  { retval = dyDUAL ; }
	  break ; }
	case lpPUNT:
        case lpSTALLED:
	{ if (actcnt == 0)
	  { retval = dyFORCEDUAL ; }
	  else
	  { retval = dyPRIMAL2 ;
	    dy_lp->simplex.next = dyPRIMAL2 ; }
	  break ; }
	case lpFORCEDUAL:
	{ retval = dyDUAL ;
	  dy_lp->simplex.next = dyDUAL ;
	  break ; }
	default:
	{ break ; } }
      break ; }
    case dyDUAL:
    { switch (dy_lp->lpret)
      { case lpOPTIMAL:
	{ if (dy_lp->simplex.next == dyPRIMAL2)
	  { if (actcnt == 0)
	    { retval = dyGENCON ;
	      dy_lp->simplex.next = dyDUAL ; }
	    else
	    { retval = dyPURGEVAR ; } }
	  else
	  if (dy_lp->simplex.next == dyDUAL)
	  { retval = dyDUAL ; }
	  break ; }
	case lpINFEAS: /* dual unbounded */
	{ retval = dyPRIMAL1 ;
	  dy_lp->simplex.next = dyPRIMAL1 ;
	  break ; }
	case lpSWING:
	{ if (actcnt == 0 && dy_lp->simplex.next != dyPRIMAL1)
	  { dy_lp->simplex.next = dyPRIMAL1 ;
	    retval = dyADDVAR ; }
	  else
	  { retval = dy_lp->simplex.next ; }
	  break ; }
	case lpLOSTFEAS:
	{ retval = dy_lp->simplex.next ;
	  break ; }
        case lpFORCEPRIMAL:
	{ retval = dyPURGEVAR ;
	  break ; }
	default:
	{ break ; } }
      break ; }
    default:
    { break ; } }

/*
  The call to dy_accchk is needed to ensure that the primal infeasibility
  information is consistent before primal1 calls dy_initp1obj.
*/
  if (dy_lp->simplex.next == dyPRIMAL1)
  { if (dy_accchk(&chkflgs) != dyrOK)
    { retval = dyINV ; } }

  if (retval == dyINV)
  { errmsg(435,rtnnme,
	   dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,
	   actcnt,"variables",dy_prtlpphase(dy_lp->simplex.active,TRUE),
	   dy_prtlpret(dy_lp->lpret),
	   dy_prtlpphase(dy_lp->simplex.next,TRUE)) ; }
  
  return (retval) ; }



static dyphase_enum initial_activation (lpprob_struct *orig_lp)

/*
  This routine handles requests for constraint and/or variable activation
  prior to starting simplex iterations. Either alone could be fit into the
  standard state cycle for dylp, but when both are specified it only makes
  sense to do constraint activation first, and that doesn't fit well with
  the standard state cycle.

  Parameter:
    orig_lp:	the original lp

  Returns: dyDUAL, dyPRIMAL1, or dyPRIMAL2 if all goes well, dyINV if there's
	   an error.
*/

{ int xindx,xipos,conresult,varresult ;
  flags xistatus,calcflgs ;
  consys_struct *orig_sys ;
  dyret_enum retval ;
  const char *rtnnme = "initial_activation" ;

  orig_sys = orig_lp->consys ;
  conresult = -1 ;
  varresult = -1 ;
/*
  If the client has asked for constraint activation, do that first. If we
  activate any constraints, we've lost primal feasibility, but there's still
  hope for dual feasibility.
*/
  if (flgon(orig_lp->ctlopts,lpctlINITACTCON))
  {
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.conmgmt >= 1)
    { dyio_outfmt(dy_logchn,dy_gtxecho,
		  "\n  [%s](%s)%d: pre-activating violated constraints, ",
		  dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),
		  dy_lp->tot.iters) ;
      dyio_outfmt(dy_logchn,dy_gtxecho,"obj = %g ...",dy_lp->z) ; }
#   endif
    conresult = dy_activateCons(orig_sys,FALSE) ;
    if (conresult > 0)
    { calcflgs = ladPRIMFEAS|ladPFQUIET|ladDUALFEAS|ladDFQUIET ;
      retval = dy_accchk(&calcflgs) ;
      if (retval != dyrOK)
      { errmsg(304,rtnnme,dy_sys->nme,
	       dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters) ;
	return (dyINV) ; }
      if (flgoff(calcflgs,ladPRIMFEAS))
      { dy_lp->simplex.next = dyPRIMAL2 ; }
      else
      if (flgoff(calcflgs,ladDUALFEAS))
      { dy_lp->simplex.next = dyDUAL ; }
      else
      { dy_lp->simplex.next = dyPRIMAL1 ; } } }
  else
  { conresult = 0 ; }
/*
  Now do variable activation, if requested. If we're targetting dual simplex,
  there's a little work involved to identify the infeasible basic primal
  variables. For each one, set ubnd.ndx to tell dualaddvars which row to
  check.
*/
  if (conresult >= 0 && flgon(orig_lp->ctlopts,lpctlINITACTVAR))
  {
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.varmgmt >= 1)
    { dyio_outfmt(dy_logchn,dy_gtxecho,
		  "\n%s: [%s](%s)%d: pre-activating variables, ",
		  rtnnme,dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),
		  dy_lp->tot.iters) ;
      dyio_outfmt(dy_logchn,dy_gtxecho,"%s rules, obj = %g ...",
		  dy_prtlpphase(dy_lp->simplex.next,TRUE),dy_lp->z) ; }
#   endif
    if (dy_lp->simplex.next == dyPRIMAL1 || dy_lp->simplex.next == dyPRIMAL2)
    { varresult = dy_activateVars(orig_sys,NULL) ; }
    else
    { for (xipos = 1 ; xipos <= dy_sys->concnt ; xipos++)
      { xindx = dy_basis[xipos] ;
	xistatus = dy_status[xindx] ;
	if (flgon(xistatus,vstatBLLB|vstatBUUB))
	{ if (flgon(xistatus,vstatBUUB)) xindx = -xindx ;
	  dy_lp->ubnd.ndx = xindx ;
	  varresult = dy_dualaddvars(orig_sys) ;
	  if (varresult < 0) break ; } }
      dy_lp->ubnd.ndx = 0 ; } }
/*
  Figure out the appropriate return value and we're done. Unless something's
  gone wrong, we want to head for the initial simplex phase.
*/
  if (conresult < 0 || varresult < 0)
    return (dyINV) ;
  else
    return (dy_lp->simplex.next) ; }



static bool commonstart (start_enum start)

/*
  This routine contains common initialisation actions for the antidegeneracy
  and pivot rejection algorithms, and DSE and PSE pricing. The allocation is
  unneeded for a hot start, but various other bits and pieces are necessary.

  Parameters:
    start:	type of startup sequence (hot, warm, cold)

  Returns: TRUE if the initialisation completes without error, FALSE
	   otherwise.
*/

{ const char *rtnnme = "commonstart" ;

/*
  Create and attach dy_brkout and dy_degenset (primal anti-degeneracy
  algorithm structures) and dy_ddegenset (dual anti-degeneracy structure).
*/
  if (start != startHOT)
  { if (consys_attach(dy_sys,CONSYS_COL,
		      sizeof(int),(void **) &dy_brkout) == FALSE)
    { errmsg(100,rtnnme,dy_sys->nme,"breakout vector") ;
      return (FALSE) ; }
    if (consys_attach(dy_sys,CONSYS_COL,
		      sizeof(int),(void **) &dy_degenset) == FALSE)
    { errmsg(100,rtnnme,dy_sys->nme,"primal degenerate set vector") ;
      return (FALSE) ; }
    if (consys_attach(dy_sys,CONSYS_ROW,
		      sizeof(int),(void **) &dy_ddegenset) == FALSE)
    { errmsg(100,rtnnme,dy_sys->nme,"dual degenerate set vector") ;
      return (FALSE) ; }
      }
  dy_lp->degen = 0 ;
/*
  Allocate the initial pivot rejection list.
*/
  if (start != startHOT)
  { dy_initpivrej(dy_sys->varcnt/10) ; }
/*
  Create and attach dy_gamma, and dy_frame (PSE structures), and dy_rho (DSE
  structure).
*/
  if (start != startHOT)
  { if (consys_attach(dy_sys,CONSYS_ROW,
		      sizeof(double),(void **) &dy_gamma) == FALSE)
    { errmsg(100,rtnnme,dy_sys->nme,"column norm vector") ;
      return (FALSE) ; }
    if (consys_attach(dy_sys,CONSYS_ROW,
		      sizeof(bool),(void **) &dy_frame) == FALSE)
    { errmsg(100,rtnnme,dy_sys->nme,"reference frame vector") ;
      return (FALSE) ; }
    if (consys_attach(dy_sys,CONSYS_COL,
		      sizeof(double),(void **) &dy_rho) == FALSE)
    { errmsg(100,rtnnme,dy_sys->nme,"basis inverse row norm") ;
      return (FALSE) ; } }
/*
  Set up for steepest edge pricing. If we're headed into primal simplex, scan
  the status array and initialise the reference frame to the current nonbasic
  variables. If we're headed into dual simplex, we need the basis inverse row
  norms; call dy_dseinit to deal with it.
*/
  if (dy_lp->simplex.next == dyDUAL)
  { dy_dseinit() ;
    dy_lp->simplex.init_dse = FALSE ;
    dy_lp->simplex.init_pse = TRUE ; }
  else
  { dy_pseinit() ;
    dy_lp->simplex.init_pse = FALSE ;
    dy_lp->simplex.init_dse = TRUE ; }

  return (TRUE) ; }



lpret_enum dylp (lpprob_struct *orig_lp, lpopts_struct *orig_opts,
	         lptols_struct *orig_tols, lpstats_struct *orig_stats)

/*
  This is the top-level routine for dylp. It orchestrates the main flow through
  the 8 step dynamic simplex algorithm outlined at the beginning of the file.

  Note that orig_sys should NOT have logical variables associated with it when
  passed to dylp. They will be added in dy_sys, but orig_sys doesn't need them
  and they just get in the way.

  If orig_opts->context == cxUNLOAD, dylp will free any retained data structures
  and return. The return code will always be lpINV.
  If orig_opts->context == cxLOAD or cxUSERPIV, dylp will do as much as is
  necessary to prepare for pivoting, then return.
  
  Parameters:
    orig_lp:	(i) The original problem --- at least the constraint system,
		    and possibly an initial basis.
		(o) The solution to the problem, including status and basis
		    vectors and values for the primal and dual variables, as
		    appropriate to the problem status.
    orig_opts:	(read only) dylp option values. The global dy_opts is used to
		pass these around to the various modules of dylp.
    orig_tols:	dylp tolerance values. Copied to an internal structure as
		some of these are written during execution.

      The orig_stats structure is used only if DYLP_STATISTICS is defined.

    orig_stats: (i) A statistics structure; may be null, in which case no
		    statistics are collected.
		(o) Statistics on the dylp run.

  Returns: Any of lpOPTIMAL, lpUNBOUNDED, or lpINFEAS can be returned when
	   dylp runs to completion. The remaining codes (see dylp.h for
	   details) indicate a forced stop (lpITERLIM) or some extraordinary
	   event.
*/

{ int cnt ;
  dyphase_enum phase ;
  double tol ;
  lpret_enum lpresult ;
  flags checks ;
  consys_struct *orig_sys ;
  const char *rtnnme = "dylp" ;

  start_enum start ;

  /* dy_force.c */
  dyphase_enum dy_forcePrimal2Dual(consys_struct *orig_sys) ;
  dyphase_enum dy_forceDual2Primal(consys_struct *orig_sys) ;
  dyphase_enum dy_forceFull(consys_struct *orig_sys) ;

# ifdef DYLP_PARANOIA
  if (orig_lp == NULL)
  { errmsg(2,rtnnme,"orig_lp") ;
    return (lpINV) ; }
  if (orig_opts == NULL)
  { errmsg(2,rtnnme,"orig_opts") ;
    return (lpINV) ; }
# endif
/*
  The first possibility is that this call is solely for the purpose of
  freeing the problem data structures. The indication is a context of
  cxUNLOAD in orig_lp->context.  Note that there may be multiple independent
  clients (objects) using dylp, and clients can make multiple calls to free
  data structures or attempt to free data structures even though dylp has
  never been called and the client doesn't own the solver.

  Note that dylp does not prevent one client from forcibly evicting another!
*/
  if (orig_opts->context == cxUNLOAD)
  { dy_finishup(orig_lp,dyINV) ;
    return (lpINV) ; }

# ifdef DYLP_PARANOIA
/*
  Check the rest of the parameters now, we'll need them unless something goes
  wrong. Paranoid checks don't attempt to clean up.
*/
  if (orig_lp->consys == NULL)
  { errmsg(2,rtnnme,"orig_sys") ;
    return (lpINV) ; }
  if (orig_tols == NULL)
  { errmsg(2,rtnnme,"orig_tols") ;
    return (lpINV) ; }
# endif

/*
  We're here to do actual work. If the constraint system is marked as corrupt,
  free any retained structures, signal an error, and bail out immediately.
  Attempting cleanup is risky here, and in some other errors below, but the
  alternative is to leak a lot of allocated space. Limit the exposure from now
  on by only attempting cleanup if a client owns the solver. (Note that
  dy_finishup will clear dy_owner.)
*/
  orig_lp->phase = dyINV ;
  orig_sys = orig_lp->consys ;
  if (flgon(orig_sys->opts,CONSYS_CORRUPT))
  { errmsg(115,rtnnme,orig_sys->nme) ;
    if (dy_owner != NULL)
    { 
#     ifndef DYLP_NDEBUG
      if (orig_opts->print.major >= 1)
	dyio_outfmt(dy_logchn,dy_gtxecho,
	  "\n  Attempting cleanup of retained structures, owner %#08x.",
	  dy_owner) ;
#     endif
      orig_lp->lpret = lpFATAL ;
      dy_finishup(orig_lp,dyINV) ; }
    return (lpFATAL) ; }
/*
  Next we need to check the forcewarm and forcecold options, and set start
  accordingly. Cold dominates warm dominates hot, for convenience of use and
  debugging. (Forcing a cold start puts a sort of limited firewall between
  dylp and the client).

  If the caller is trying for a hot start, do some quick checks --- the
  client should own the solver (hence we have retained data structures), and
  the previous run should have ended cleanly (i.e., with a result of optimal,
  unbounded, infeasible, or iteration limit).  Iteration limit is (sort of)
  a special case -- this arises frequently when a B&C code is using dylp
  for strong branching; the iteration limit will be deliberately inadequate.

  If we're paranoid, do a few additional consistency checks.
*/
# ifdef DYLP_PARANOIA
  if ((flgoff(orig_lp->ctlopts,lpctlDYVALID) && dy_owner != NULL) ||
      (flgon(orig_lp->ctlopts,lpctlDYVALID) && dy_owner == NULL))
  { errmsg(1,rtnnme,__LINE__) ;
    return (lpINV) ; }
# endif

  if (orig_opts->forcecold == TRUE)
  { start = startCOLD ; }
  else
  if (orig_opts->forcewarm == TRUE)
  { start = startWARM ; }
  else
  { start = startHOT ;
    if (dy_owner != orig_lp->owner)
    { errmsg(396,rtnnme,orig_sys->nme,orig_lp->owner,dy_owner,"hot start") ;
      return (lpINV) ; }
    if (!(orig_lp->lpret == lpOPTIMAL || orig_lp->lpret == lpUNBOUNDED ||
	  orig_lp->lpret == lpINFEAS || orig_lp->lpret == lpITERLIM))
    { errmsg(395,rtnnme,orig_sys->nme,dy_prtlpret(orig_lp->lpret)) ;
      return (lpINV) ; } }

# ifndef DYLP_NDEBUG
  if (orig_opts->print.major >= 1)
    dyio_outfmt(dy_logchn,dy_gtxecho,"\n  %s start for lp %s.",
	        (start == startHOT)?"hot":((start == startWARM)?"warm":"cold"),
	        orig_sys->nme) ;
# endif
# ifdef DYLP_PARANOIA
/*
  In the context of an B&C code, it's common to run an LP to check a
  solution. If many or all variables are fixed, a presolve phase may give a
  0x0 system. Warn about it if we're paranoid.
*/
  if (orig_sys->concnt < 1 || orig_sys->varcnt < 1)
  { warn(351,rtnnme,orig_sys->nme,dy_prtlpphase(dyINV,TRUE),0,
	 orig_sys->concnt,orig_sys->varcnt) ; }
  if (flgon(orig_sys->opts,CONSYS_LVARS))
  { errmsg(123,rtnnme,orig_sys->nme) ;
    return (lpINV) ; }
# endif

/*
  Let's get to it. Do a cleanup and initialise the dy_lp control structure
  unless we're doing a hot start. In any event, mark the previous results as
  invalid, as we're about to start manipulating them.
*/
  orig_lp->phase = dyINIT ;
  orig_lp->lpret = lpINV ;
  lpresult = lpINV ;
  if (start != startHOT)
  { dy_finishup(orig_lp,dyINV) ;
    dy_lp = (lp_struct *) CALLOC(1,sizeof(lp_struct)) ;
    dy_lp->p1obj.installed = FALSE ;
    dy_lp->p1obj.infvars = NULL ;
    dy_lp->p1obj.p1obj = NULL ;
    dy_lp->p1obj.p2obj = NULL ; }
  dy_lp->phase = dyINIT ;
  dy_lp->lpret = lpINV ;
  clrflg(orig_lp->ctlopts,lpctlDYVALID) ;
/*
  If we're doing a warm or cold start, the first order of business is to
  establish the environment -- options and tolerances that'll be used by dylp
  as it works. For a hot start, these are already in place, but we allow the
  user to change them (give `em more than enough rope, I say ... ). Same for
  context.
*/
  updateOptsAndTols(orig_opts,orig_tols) ;
# ifdef DYLP_STATISTICS
  dy_stats = orig_stats ;
# endif
/*
  Initialise the local constraint system, if required. Scaling occurs here,
  if allowed. The original system is hidden away and orig_lp->consys is
  replaced by a scaled copy. See dy_scaling for details.
*/
  if (dy_initlclsystem(orig_lp,(start == startHOT)?TRUE:FALSE) != TRUE)
  { errmsg(406,rtnnme,orig_sys->nme) ;
    orig_lp->lpret = lpFATAL ;
    dy_finishup(orig_lp,dy_lp->phase) ;
    return (lpFATAL) ; }
  orig_sys = orig_lp->consys ;
/*
  Pick the appropriate startup routine for a hot, warm, or cold start.

  With a cold start, the user hasn't supplied an initial basis. Call
  dy_startup and dy_crash to get going. With a warm start, the user's given
  us a basis, but we need to call dy_warmstart to build the dylp data
  structures.  If the previous problem was left active, we use dy_hotstart.

  On completion, the active constraint system and cross-reference structures
  will be built, a basis will be established and factored, the status vector
  will be valid, and values of the primal and dual variables and reduced costs
  will be valid.

  If the problem is prima facie infeasible (lower and upper bounds cross for
  some variable or constraint) this is reported via dy_lp->lpret, and we are
  immediately done. Similarly for unboundedness (an unconstrained variable
  with nonzero objective coefficient).
*/
  switch (start)
  { case startHOT:
    { if (dy_hotstart(orig_lp) != dyrOK)
      { errmsg(371,rtnnme,
	       orig_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),0,"hot start") ;
	dy_lp->lpret = lpFATAL ; }
      break ; }
    case startWARM:
    { if (dy_warmstart(orig_lp) != dyrOK)
      { errmsg(371,rtnnme,
	       orig_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),0,"warm start") ;
	dy_lp->lpret = lpFATAL ; }
      break ; }
    case startCOLD:
    { if (dy_coldstart(orig_sys) != dyrOK)
      { errmsg(371,rtnnme,
	       orig_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),0,"cold start") ;
	dy_lp->lpret = lpFATAL ; }
      else
      if (dy_crash() != dyrOK)
      { errmsg(302,rtnnme,dy_sys->nme,
	       dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,"crash") ;
	dy_lp->lpret = lpFATAL ; }
      break ; } }
  if (dy_lp->lpret != lpINV)
  { orig_lp->lpret = dy_lp->lpret ;
    if (orig_lp->lpret == lpINFEAS || orig_lp->lpret == lpUNBOUNDED)
    { dy_lp->phase = dyDONE ;
      orig_lp->phase = dyDONE ;
#ifndef DYLP_NDEBUG
      if (dy_opts->print.major >= 1)
      { dyio_outfmt(dy_logchn,dy_gtxecho,"\n\n%s (%s): prima facie %s.",
	    rtnnme,dy_sys->nme,
	    ((orig_lp->lpret == lpINFEAS)?"infeasibility":"unboundedness")) ; }
#   endif
    }
    dy_finishup(orig_lp,dy_lp->phase) ;
    return (orig_lp->lpret) ; }
/*
  Do a little more setup prior to invoking a simplex. commonstart handles
  initial setup for the antidegeneracy and pivot rejection algorithms,
  PSE and DSE pricing, and reduced costs.
*/
  phase = dy_lp->simplex.next ;
  if (dy_opts->usedual == FALSE && phase == dyDUAL) phase = dyPRIMAL1 ;
  dy_lp->simplex.active = dyINV ;
  dy_lp->simplex.init_dse = TRUE ;
  dy_lp->simplex.init_pse = TRUE ;
  if (commonstart(start) == FALSE)
  { errmsg(371,rtnnme,
	   orig_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),0,"common start") ;
    orig_lp->lpret = lpFATAL ;
    dy_finishup(orig_lp,dy_lp->phase) ;
    return (lpFATAL) ; }
/*
  lastz.cd and lastz.vd control whether we do constraint and variable
  purging, respectively.  Objective change is not, in general, monotonic, so
  we're just looking to get sufficient change to be sure we've moved since
  the last purge. ubnd.ndx is the index of the primal variable fingered for
  primal unboundedness.
*/
  dy_lp->lastz.cd = -dy_tols->inf ;
  dy_lp->lastz.vd = -dy_tols->inf ;
  dy_lp->ubnd.ndx = 0 ;
/*
  Various special startup activities. These are mutually exclusive, and need
  to be checked in the order given.

  Is this one of those pathological cases with no constraints? If so, make it
  look like we're optimal. No constraints doesn't mean no variables, so
  head for dyGENVAR. 
*/
  if (dy_sys->concnt == 0)
  { dy_lp->simplex.active = dyPRIMAL2 ;
    dy_lp->simplex.next = dyPRIMAL2 ;
    dy_lp->lpret = lpOPTIMAL ;
    phase = dyGENVAR ; }
/*
  Do we want to do an initial variable purge, to try and cut down the number
  of variables? (This is particularly aimed at large set covering problems
  where all constraints are equalities and there are thousands of columns.
  dylp's normal cold start procedure will load the full system.)

  We need to rerun dy_initp1obj solely because deletion of variables can cause
  other variables to move, and this might invalidate indices in the infvars
  vector.
*/
  else
  if (start == startCOLD && dy_opts->fullsys == FALSE &&
      dy_sys->archvcnt > dy_opts->coldvars)
  { cnt = dy_deactivateVars(orig_sys) ;
    if (cnt < 0)
    { errmsg(371,rtnnme,orig_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),0,
	     "initial variable deactivation") ; }
    if (phase == dyPRIMAL1 && cnt > 0)
    { if (dy_initp1obj() == FALSE)
      { errmsg(318,rtnnme,dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),
	       dy_lp->tot.iters,"initialise") ;
	phase = dyINV ; } } }
/*
  Is the user requesting initial constraint and/or variable activation? This
  is a warm or hot start activity (cold start has its own rules to decide how
  much of the system to activate). A fatal error here will return dyINV for
  the phase.
*/
  else
  if (start == startWARM || start == startHOT)
  { if (dy_opts->fullsys == FALSE &&
	flgon(orig_lp->ctlopts,lpctlINITACTCON|lpctlINITACTVAR))
    { phase = initial_activation(orig_lp) ; } }
# ifdef DYLP_STATISTICS
  if (dy_stats != NULL) dy_stats->ini_simplex = phase ;
# endif
/*
  Open the main loop for dylp.  Generally, there are two minor cycles ---
  purge/generate/add variables, primal simplex; and purge/generate/add
  constraints, dual simplex --- within an outer cycle which alternates
  between the primal and dual subcycles. There are, of course, complications
  in handling unbounded and infeasible problems. phase tracks the current
  algorithm phase.
*/
  while (phase != dyDONE && phase != dyINV)
  { dy_lp->phase = phase ;
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.major >= 1)
      dyio_outfmt(dy_logchn,dy_gtxecho,
		  "\n\n%s (%s): entering phase %s, iter %d.",rtnnme,
		  dy_sys->nme,dy_prtlpphase(phase,FALSE),dy_lp->tot.iters) ;
#   endif
#   ifdef DYLP_STATISTICS
    if (dy_stats != NULL) dy_stats->phasecnts[phase]++ ;
#   endif

    switch (phase)
    {
/*
  Under normal circumstances, we'll enter dyPRIMAL[1|2] from dyINIT to solve
  the initial LP. We'll also come here from dyADDVAR, to reoptimise after
  adding variables, and from dyADDCONS, to reoptimise after adding
  constraints to an unbounded primal. There are a number of more exotic error
  recovery paths which lead back to primal simplex: primal phase I is the
  fallback when we can't recover primal or dual feasibility.

  Depending on where we've been, it may be necessary to reset the PSE
  reference frame before starting simplex iterations (simplex.init_pse ==
  TRUE). In particular, if we've been running dual simplex, or have added
  constraints, we need to do a reset. Since we'll change the basis, we'll
  need to to a DSE reset if we ever return to dual simplex.

  If the result of the LP is lpOPTIMAL, we'll do a gen/add/purge constraints
  sequence, then dual simplex. If the result is lpINFEAS, we'll do a
  generate/add variables sequence, returning here if new variables were
  added.

  If the result is unbounded, we'll do a generate/add constraints sequence,
  returning here to try again.  It's unlikely that phase I will go unbounded,
  but we're only dealing with a subset of the constraints, and it happens on
  occasion.

  If we return a punt or stall, it's an indication that we're in trouble down
  in the simplex. As a first cut, head off to check the variables. Maybe we
  can activate some desireable candidates that will allow us to pivot past
  this point.

  Anything else is an error.
*/
      case dyPRIMAL1:
      case dyPRIMAL2:
      {
	dy_lp->simplex.active = phase ;
	if (dy_lp->simplex.init_pse == TRUE)
	{ dy_pseinit() ;
	  dy_lp->simplex.init_pse = FALSE ; }
	dy_lp->simplex.init_dse = TRUE ;
	lpresult = dy_primal() ;
	dy_lp->simplex.next = dy_lp->phase ; 
	if (dy_opts->context == cxLOAD)
	{ phase = dyDONE ; }
	else
        { switch (lpresult)
	  { case lpOPTIMAL:
	    { phase = dyGENVAR ;
	      dy_lp->simplex.next = dyPRIMAL2 ;
	      break ; }
	    case lpINFEAS:
	    { phase = dyGENVAR ;
	      break ; }
	    case lpUNBOUNDED:
	    { if (dy_sys->concnt < orig_sys->concnt)
	      { phase = dyGENCON ; }
	      else
	      { phase = dyDONE ; }
  #	    ifdef DYLP_PARANOIA
	      if (dy_lp->ubnd.ndx == 0)
	      { errmsg(1,rtnnme,__LINE__) ;
		phase = dyINV ;
		break ; }
  #	    endif
	      break ; }
	    case lpSWING:
	    { phase = dyGENCON ;
  #	    ifdef DYLP_PARANOIA
	      if (dy_lp->ubnd.ndx == 0)
	      { errmsg(1,rtnnme,__LINE__) ;
		phase = dyINV ;
		break ; }
  #	    endif
	      break ; }
	    case lpPUNT:
	    case lpSTALLED:
	    { phase = dyGENVAR ;
	      break ; }
	    case lpACCCHK:
	    { phase = dyFORCEFULL ;
	      break ; }
	    default:
	    { if (!(dy_opts->context == cxBANDC && lpresult == lpITERLIM))
	      { errmsg(353,rtnnme,orig_sys->nme,"primal",
		       dy_prtlpret(lpresult)) ; }
	      phase = dyDONE ;
	      break ; } } }
	break ; }
/*
  In the best case, dual simplex has reported optimal, we've added variables
  with a  gen/add variables sequence, and now we're doing a purge before
  returning to primal simplex.  Call dy_deactivateVars to remove variables
  with sufficiently lousy reduced costs (in the sense that it's highly
  unlikely they'll return to the basis in the optimal solution).

  In the not so good case, we've been forced out of dual simplex for some
  reason and may or may not be primal feasible.

  There's no guarantee of monotonic change in the objective from one purge
  opportunity to the next. We just need to know we've moved.
*/
      case dyPURGEVAR:
      { tol = dy_tols->purge*(1.0+fabs(dy_lp->z)) ;
	if (!withintol(dy_lp->z,dy_lp->lastz.vd,tol))
	{
#	  ifndef DYLP_NDEBUG
	  if (dy_opts->print.varmgmt >= 1)
	  { dyio_outfmt(dy_logchn,dy_gtxecho,
		        "\n  [%s](%s)%d: purging variables, obj = %g ...",
		        dy_sys->nme,dy_prtlpphase(phase,TRUE),dy_lp->tot.iters,
		        dy_lp->z) ; }
#	  endif
	  cnt = dy_deactivateVars(orig_sys) ;
	  dy_lp->lastz.vd = dy_lp->z ;
	  if (cnt < 0)
	  { phase = dyINV ; }
	  else
	  { phase = dy_lp->simplex.next ; } }
	else
	{
#	  ifndef DYLP_NDEBUG
	  if (dy_opts->print.varmgmt >= 2)
	  { dyio_outfmt(dy_logchn,dy_gtxecho,
		        "\n  [%s](%s)%d: variable purge skipped, obj = %g.",
		        dy_sys->nme,dy_prtlpphase(phase,TRUE),dy_lp->tot.iters,
		        dy_lp->z) ; }
#	  endif
	  phase = dy_lp->simplex.next ; }
	break ; }
/*
  This phase should call a routine that generates new variables, not already
  present in orig_sys. This would be used for column generation algorithms
  and similar sorts of things. At the end of the routine, orig_sys should be
  modified such that the generated variables are part of the set of inactive
  variables, where they will be priced and added if advantageous during the
  dyADDVAR phase. See comments for dyADDVAR.
*/
      case dyGENVAR:
      { phase = dyADDVAR ;
	break ; }
/*
  This phase scans the inactive variables and adds the ones that look
  useful.  The definition of useful varies, depending on the target simplex
  phase specified by dy_lp->simplex.next --- the simplex that will
  (eventually) run after dy_activateVars completes.

  If we're targetting primal simplex, we're looking for variables that price
  out as favourable (non-optimal) for their current status. The only
  complication is whether we're heading for primal phase I or II, and
  dy_activateVars will make sure it's using the correct reduced costs.

  If we're here because the dual is unbounded, we're really looking for
  columns which, when considered as dual constraints, can bound the dual.
  dy_dualaddvars will look for useful variables, considering progressively
  more exotic possibilities.

  Where we go from here depends on how we got here and how many variables
  were activated. See the comments with addvar_nextphase. If all is going
  well and we activated variables of the desired type, we'll head for the
  target simplex phase. If we didn't find variables, the possibilities are
  many and varied.

  Why would no variables be activated? A straightforward cause is there are
  no variables of the appropriate type (e.g., we're optimal or infeasible).
  Then there are things that dy_dualaddvars (more accurately, dual simplex)
  just can't cope with (e.g., no variables can be activated without losing
  dual feasibility).

  There are a number of other exotic things that can go wrong. See the
  comments with addvar_nextphase.
*/
      case dyADDVAR:
      { if (dy_lp->sys.vars.loadable <= 0)
	{ 
#	  ifndef DYLP_NDEBUG
	  if (dy_opts->print.varmgmt >= 1)
	  { dyio_outfmt(dy_logchn,dy_gtxecho,
		        "\n  [%s](%s)%d: no loadable variables; skipping.",
		        dy_sys->nme,dy_prtlpphase(phase,TRUE),
			dy_lp->tot.iters) ; }
#	  endif
	  phase = addvar_nextphase(0) ;
	  break ; }
	else
	{ if (dy_lp->simplex.next == dyDUAL && dy_lp->lpret == lpINFEAS)
	  { if (dy_opts->dualadd > 0)
	    {
#	      ifndef DYLP_NDEBUG
	      if (dy_opts->print.varmgmt >= 1)
	      { dyio_outfmt(dy_logchn,dy_gtxecho,
			    "\n%s: [%s](%s)%d: activating variables, ",rtnnme,
			    dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),
			    dy_lp->tot.iters) ;
		dyio_outfmt(dy_logchn,dy_gtxecho,"%s rules, obj = %g ...",
			    dy_prtlpphase(dy_lp->simplex.next,TRUE),
			    dy_lp->z) ; }
#	      endif
	      cnt = dy_dualaddvars(orig_sys) ; }
	    else
	    {  cnt = 0 ; }
	    if (cnt > 0)
	    { phase = dyDUAL ; }
	    else
	    if (cnt == 0)
	    { dy_lp->simplex.next = dyPRIMAL1 ; 
	      checks = ladPRIMFEAS|ladPFQUIET ;
	      if (dy_accchk(&checks) != dyrOK)
	      { phase = dyINV ; } }
	    else
	    { phase = dyINV ; } }
	  if (phase == dyADDVAR)
	  { 
#	    ifndef DYLP_NDEBUG
	    if (dy_opts->print.varmgmt >= 1)
	    { dyio_outfmt(dy_logchn,dy_gtxecho,
			  "\n%s: [%s](%s)%d: activating variables, ",
			  rtnnme,dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),
			  dy_lp->tot.iters) ;
	      dyio_outfmt(dy_logchn,dy_gtxecho,"%s rules, obj = %g ...",
			  dy_prtlpphase(dy_lp->simplex.next,TRUE),dy_lp->z) ; }
#	    endif
	    cnt = dy_activateVars(orig_sys,NULL) ;
	    phase = addvar_nextphase(cnt) ; } }
	break ; }
/*
  To arrive here, we've solved to optimality with primal simplex, added
  constraints with a  gen/add constraints, and are now purging loose
  constraints (dual variables with unfavourable dual reduced costs) before
  heading for dual simplex.

  This phase calls dy_deactivateCons to scan the active constraints and
  remove any that are loose. This is done only when there's been
  strict degradation of the objective (i.e., we've cut off the previous
  optimum point) to minimise the chance of cycling. This will not, however,
  guarantee that a constraint will never reactivate. A loose constraint can
  become tight again when variables are activated and change value.
*/
      case dyPURGECON:
      { tol = dy_tols->purge*(1.0+fabs(dy_lp->z)) ;
	if (dy_lp->z-dy_lp->lastz.cd > tol)
	{
#	  ifndef DYLP_NDEBUG
	  if (dy_opts->print.conmgmt >= 1)
	  { dyio_outfmt(dy_logchn,dy_gtxecho,
		        "\n  [%s](%s)%d: purging constraints, obj = %g ...",
		        dy_sys->nme,dy_prtlpphase(phase,TRUE),dy_lp->tot.iters,
		        dy_lp->z) ; }
#	  endif
	  cnt = dy_deactivateCons(orig_sys) ;
	  if (cnt < 0)
	    phase = dyINV ;
	  else
	    phase = dyDUAL ;
	  dy_lp->lastz.cd = dy_lp->z ;
# 	  ifdef DYLP_PARANOIA
	  if (dy_chkdysys(orig_sys) == FALSE) phase = dyINV ;
# 	  endif
	}
	else
	{
#	  ifndef DYLP_NDEBUG
	  if (dy_opts->print.conmgmt >= 2)
	  { dyio_outfmt(dy_logchn,dy_gtxecho,
		        "\n  [%s](%s)%d: constraint purge skipped, obj = %g.",
		        dy_sys->nme,dy_prtlpphase(phase,TRUE),dy_lp->tot.iters,
		        dy_lp->z) ; }
#	  endif
	  phase = dyDUAL ; }
	break ; }
/*
  This phase should call a routine that generates new constraints --- cutting
  planes, for example. As with dyGENVAR, these constraints should be added to
  the original system, so that they can be found in phase dyADDCON. See the
  comments with dyADDCON.
*/
      case dyGENCON:
      { phase = dyADDCON ;
	break ; }
/*
  In the common case, we're adding constraints before starting dual
  simplex. Also normal, and all too common when working with a partial
  constraint system, we're here to look for constraints to bound the primal.

  dy_activateBndCons scans the inactive constraints for constraints that
  could bound the problem but are not violated.  dy_activateCons scans the
  inactive constraints for violated constraints. When trying to bound a
  problem, it's fairly common that there are no nonviolated bounding
  constraints. In this case, we'll go for violated constraints. (Hence the
  seemingly redundant test around dy_activateCons.)

  There's no sense in even bothering to scan if there are no constraints to
  load.

  Inactive variables referenced by activated constraints can themselves be
  activated. Care is taken to preserve dual feasibility (primal feasibility
  cannot be lost).
*/
      case dyADDCON:
      { if (dy_lp->sys.cons.loadable <= 0)
	{ phase = addcon_nextphase(0) ;
#	  ifndef DYLP_NDEBUG
	  if (dy_opts->print.conmgmt >= 1)
	  { dyio_outfmt(dy_logchn,dy_gtxecho,
		        "\n  [%s](%s)%d: no loadable constraints; skipping.",
		        dy_sys->nme,dy_prtlpphase(phase,TRUE),
			dy_lp->tot.iters) ; }
#	  endif
	  break ; }
	if (dy_lp->lpret == lpSWING ||
	    (dy_lp->lpret == lpUNBOUNDED &&
	      (dy_lp->simplex.next == dyPRIMAL2 ||
	       dy_lp->simplex.next == dyPRIMAL1)))
	{
#         ifndef DYLP_NDEBUG
	  if (dy_opts->print.conmgmt >= 1)
	  { dyio_outfmt(dy_logchn,dy_gtxecho,
		        "\n  [%s](%s)%d: activating bounding constraints, ",
		        dy_sys->nme,dy_prtlpphase(phase,TRUE),
			dy_lp->tot.iters) ;
	    dyio_outfmt(dy_logchn,dy_gtxecho,"obj = %g ...",dy_lp->z) ; }
#         endif
	  cnt = dy_activateBndCons(orig_sys) ;
	  if (cnt > 0)
	  { phase = dy_lp->simplex.next ; }
	  else
	  if (cnt == 0)
	  { dy_lp->simplex.next = dyPRIMAL1 ; }
	  else
	  { phase = dyINV ; } }
	if (phase == dyADDCON)
	{
#         ifndef DYLP_NDEBUG
	  if (dy_opts->print.conmgmt >= 1)
	  { dyio_outfmt(dy_logchn,dy_gtxecho,
		        "\n  [%s](%s)%d: activating violated constraints, ",
		        dy_sys->nme,dy_prtlpphase(phase,TRUE),
			dy_lp->tot.iters) ;
	    dyio_outfmt(dy_logchn,dy_gtxecho,"obj = %g ...",dy_lp->z) ; }
#         endif
	  cnt = dy_activateCons(orig_sys,TRUE) ;
	  phase = addcon_nextphase(cnt) ; }
	break ; }
/*
  Dual simplex has no phase I, and must start from a dual feasible basis.  It
  can happen that the startup routines generate a dual feasible basis in
  dyINIT.  Most often, dual feasibility is obtained by running primal simplex
  to optimality.

  If dylp is running normally, we arrive here after a gen/add/purge
  constraints sequence and will use dual simplex to reoptimise. We can also
  arrive here after adding variables (dual constraints) to an unbounded
  dual.  There are a number of more exotic error recovery paths.

  Bottom line, though, is we don't get here without dual feasibility.  That
  limits the results to optimal, dual unbounded (translated to primal
  infeasible), or some sort of problem.

  If we're optimal, we'll do a gen/add/purge variables sequence, then primal
  simplex.

  If we come back dual unbounded (seen here as primal infeasible) we need to
  add constraints that will bound the dual, i.e., we need to add primal
  variables in such a way that we can return to dual simplex. Off to gen/add
  variables. If we can't add bounding constraints without losing dual
  feasibility, dylp will fall back to primal phase I.

  If we return a punt or stall, we'll try to add some dual variables (primal
  constraints). In the event of loss of dual feasibility, we'll try to
  force dual feasibility by dropping the offending variables, again falling
  back to primal phase I if we're unsuccessful.
*/
      case dyDUAL:
      { dy_lp->simplex.active = dyDUAL ;
	if (dy_lp->simplex.init_dse == TRUE)
	{ dy_dseinit() ;
	  dy_lp->simplex.init_dse = FALSE ; }
	dy_lp->simplex.init_pse = TRUE ;
	lpresult = dy_dual() ;
	dy_lp->simplex.next = dyDUAL ;
	if (dy_opts->context == cxLOAD)
	{ phase = dyDONE ; }
	else
	{ switch (lpresult)
	  { case lpOPTIMAL:
	    { phase = dyGENVAR ;
	      dy_lp->simplex.next = dyPRIMAL2 ;
	      break ; }
	    case lpINFEAS: /* dual unbounded */
	    { phase = dyGENVAR ;
	      break ; }
	    case lpLOSTFEAS:
	    { phase = dyFORCEDUAL ;
	      break ; }
	    case lpPUNT:
	    case lpSTALLED:
	    case lpSWING:
	    { phase = dyGENCON ;
	      break ; }
	    case lpACCCHK:
	    { phase = dyFORCEFULL ;
	      break ; }
	    default:
	    { if (!(dy_opts->context == cxBANDC && lpresult == lpITERLIM))
	      { errmsg(353,rtnnme,orig_sys->nme,"dual",
		       dy_prtlpret(lpresult)) ; }
	      phase = dyDONE ;
	      break ; } } }
	break ; }
      case dyFORCEDUAL:
      { phase = dy_forcePrimal2Dual(orig_sys) ;
	break ; }
      case dyFORCEPRIMAL:
      { phase = dy_forceDual2Primal(orig_sys) ;
	break ; }
      case dyFORCEFULL:
      { if (dy_lp->sys.forcedfull == TRUE ||
	    (dy_lp->sys.cons.loadable <= 0 && dy_lp->sys.vars.loadable <= 0))
	{ dy_lp->lpret = lpFORCEFULL ;
	  phase = dyDONE ; }
	else
	{ phase = dy_forceFull(orig_sys) ;
	  dy_lp->sys.forcedfull = TRUE ; }
	break ; }
      default:
      { phase = dyINV ;
	errmsg(1,rtnnme,__LINE__) ;
	break ; } } }
/*
  End of main loop on dylp phase. phase has where we're going, and should be
  one of dyDONE or dyINV; the latter means something on the order of internal
  confusion. dy_lp->phase records the phase just completed.

  First order of business is a final constraint and/or variable purge. Do this
  only if we have an optimal solution. For constraints, cut the purge level
  back to 0 (i.e., purge strictly loose constraints).
*/
  if (phase == dyDONE &&
      dy_lp->lpret == lpOPTIMAL && dy_opts->context != cxLOAD)
  { if (dy_opts->finpurge.vars == TRUE)
    {
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.major >= 1)
      { dyio_outfmt(dy_logchn,dy_gtxecho,
		    "\n\n%s (%s): entering phase %s (final), iter %d.",
		    rtnnme,dy_sys->nme,dy_prtlpphase(dyPURGEVAR,FALSE),
		    dy_lp->tot.iters) ; }
#     endif
      cnt = dy_deactivateVars(orig_sys) ;
      if (cnt < 0)
      { errmsg(371,rtnnme,
	       orig_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.pivs,
	       "final variable deactivation") ; } }
    if (dy_opts->finpurge.cons == TRUE)
    { 
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.major >= 1)
      { dyio_outfmt(dy_logchn,dy_gtxecho,
		    "\n\n%s (%s): entering phase %s (final), iter %d.",
		    rtnnme,dy_sys->nme,dy_prtlpphase(dyPURGECON,FALSE),
		    dy_lp->tot.iters) ; }
#     endif
      dy_opts->con.deactlvl = 0 ;
      cnt = dy_deactivateCons(orig_sys) ;
      if (cnt < 0)
      { errmsg(371,rtnnme,
	       orig_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.pivs,
	       "final constraint deactivation") ; } }
  }
/*
  If we're infeasible, and the phase I objective is still in place (the
  normal situation), swap it out for the original objective and recalculate
  duals and reduced costs. Otherwise our dual variable information is
  all wrong. We need to pretend the phase is dyPRIMAL2 for this to work.
*/
  if (phase == dyDONE && dy_lp->lpret == lpINFEAS)
  { phase = dyPRIMAL2 ;
    if (dy_swapobjs(dyPRIMAL2) == FALSE)
    { phase = dyDONE ;
      errmsg(318,rtnnme,dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),
	     dy_lp->tot.iters,"remove") ;
      dy_lp->lpret = lpFATAL ; }
    dy_calcduals() ;
    if (dy_calccbar() == FALSE)
    { phase = dyDONE ;
      errmsg(384,rtnnme,dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),
	     dy_lp->tot.iters) ;
      dy_lp->lpret = lpFATAL ; }
    phase = dyDONE ; }
/*
  Call dy_finishup to assemble the final answer (as best it can) and clean up
  the working environment.
*/
  orig_lp->phase = phase ;
# ifdef DYLP_STATISTICS
  if (dy_stats != NULL) dy_finalstats(dy_stats) ;
# endif
  dy_finishup(orig_lp,dy_lp->phase) ;

  return (orig_lp->lpret) ; }
 
