/*
  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 routines to populate the active constraint system and
  select an initial basis for an lp. Both of these activities are discussed in
  more detail below. For here, it suffices to say that dylp offers a great deal
  of flexibility in establishing the initial constraint system and selecting
  the initial basis.
*/

#define DYLP_INTERNAL

#include "dylp.h"

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



/*
  Cold start routines.

  dy_coldstart and subsidiary routines select a set of constraints to populate
  the initial active system. There are a number of options. But first a bit of
  general explanation.

  Ideally, we'd like to be able to guess the constraints that will be tight at
  optimum and load exactly those constraints. If we can also guess the right
  basic variables, we have the answer. It's not that easy, of course. See the
  comments in the second half of the file for an explanation of how dylp
  chooses the initial set of basic variables.

  Back to the subject at hand. Recall that we're dealing with equalities and
  <= constraints. Also recall that we're minimising. Then the normals a<i> to
  the inequality constraints point out of the feasible region, and at optimum
  the normal c to the objective will point into the feasible region. Hence
  ideal `alignment' with the objective occurs when the angle between a<i> and
  c is 180 degrees. When the angle is 0 degrees, the constraint is forming
  the far side of the polytope.

  For convenience, call the inequalities with 180 <= angle(a<i>,c) < 90 the
  `near' group, those with angle(a<i>,c) = 90 the `perp' group, and those
  with 90 < angle(a<i>,c) <= 0 the `far' group.  The implementation creates a
  sorted list of inequalities, in nonincreasing order using angle(a<i>,c),
  and marks the boundaries between the near, perp, and far groups.

  To populate the active system, equalities are loaded first. All equalities
  are loaded, always.  Empty constraints are never loaded. For the
  inequalities, it's possible to load one or two specified angular ranges,
  with a specified sampling rate. This seems to allow as much flexibility as
  is useful.

  NOTE: The routines that construct the initial basis use knowledge of the
	load order, and knowledge of how consys augments columns, to scan the
	columns in the proper order when selecting architectural variables
	for basis positions. If you change one, you must change both. See
	ib_archvselect.

  Note that specifying a range of 180 - 0 degrees with sampling rate 1.0
  (i.e., activate all the original constraints) is >not< quite the same as
  the fullsys option. 100% activation still allows initial variable
  deactivation (hence multiple variable activation/deactivation phases) and
  final variable and constraint deactivation; fullsys suppresses these.

  Recall that angle(a<i>,c) = (180/pi)*arccos(dot(a<i>,c)/(||a<i>||||c||))
*/


/*
  Utility structures

  angle_struct

  Field		Definition
  -----		----------
  ndx		constraint index
  angle		angle between a<ndx> and c

  ineq_struct

  Field		Definition
  -----		----------
  cnt		Number of inequalities
  perp		index in angles of first inequality with
		angle(a<i>,c) = 90 degrees
  far		index in angles of first inequality with
		angle(a<i>,c) < 90 degrees
  angles	inequalities, sorted by angle (indexed from 0)
*/

typedef struct { int ndx ;
		 double angle ; } angle_struct ;

typedef struct { int cnt ;
		 int perp ;
		 int far ;
		 angle_struct *angles ; } ineq_struct ;


static int near_perp_far (const void *elem1, const void *elem2)

{ const angle_struct *c1,*c2 ;
  c1 = (const angle_struct *) elem1 ;
  c2 = (const angle_struct *) elem2 ;

  if (c1->angle > c2->angle)
    return (-1) ;
  else
  if (c1->angle < c2->angle)
    return (1) ;
  else
    return (0) ; }


static bool cold_sortcons (consys_struct *orig_sys,
			   int **p_eqs, ineq_struct **p_ineqs, int **p_noload)

/*
  This routine separates the constraints into equalities and inequalities,
  dropping empty and nonbinding constraints. It also checks the bounds of
  ranged constraints for prima facie infeasibility. Depending on options,
  various additional work is performed:
    * If fullsys is not specified, the inequalities are sorted according to
      their angle from the objective function.
    * If full statistics are enabled, constraint angles are loaded into the
      statistics structure and a histogram is compiled.
    * If debug printing is enabled, the requested information is dumped to the
      log.

  Parameters:
    orig_sys:	The original constraint system
    p_eqs:	(i) empty array of int; assumed to be sufficiently large;
		    allocated if NULL
		(o) eqs[0] set to number of equalities
		    eqs[1 .. eqs[0]] set to indices of equality constraints
    p_ineqs:	(i) empty ineq_struct; assumed to be sufficiently large;
		    allocated if NULL
		(o) filled in with inequality information as described in
		    comments at head of file
    p_noload	(i) empty array of int; assumed to be sufficiently large;
		    allocated if NULL
		(o) noload[0] set to number of constraints ineligible for
		    activation; noload[1 .. noload[0]] set to indices of
		    ineligible constraints

  Return value: TRUE if all goes as planned, FALSE on error.
		dy_lp->lpret is set to lpINFEAS for prima facie infeasibility
*/

{ int i,ndx,m,n,eqcnt,ineqcnt,noloadcnt,nearcnt,perpcnt,farcnt ;
  double cnorm,ainorm,aidotc,pi180,anglei ;
  double *c ;
  int *eqs,*noload ;
  ineq_struct *ineqs ;
  angle_struct *angles ;
  contyp_enum *ctyp ;
  double *rhs,*rhslow ;

  bool retval,need_angles ;

  pkvec_struct *ai ;

  const char *rtnnme = "cold_sortcons" ;

# ifdef DYLP_STATISTICS
  int k ;
# endif

# ifdef DYLP_PARANOIA
  if (orig_sys == NULL)
  { errmsg(2,rtnnme,"orig_sys") ;
    return (FALSE) ; }
  if (p_eqs == NULL)
  { errmsg(2,rtnnme,"&eqs") ;
    return (FALSE) ; }
  if (p_ineqs == NULL)
  { errmsg(2,rtnnme,"&ineqs") ;
    return (FALSE) ; }
  if (p_noload == NULL)
  { errmsg(2,rtnnme,"&noload") ;
    return (FALSE) ; }
# endif

/*
  Allocate the equality, inequality, and noload data structures, if the
  client didn't supply them.
*/
  if (*p_eqs == NULL)
  { eqs = (int *) MALLOC((orig_sys->concnt+1)*sizeof(int)) ; }
  else
  { eqs = *p_eqs ; }
  if (*p_ineqs == NULL)
  { ineqs = CALLOC(1,sizeof(ineq_struct)) ;
    ineqs->angles =
	(angle_struct *) MALLOC(orig_sys->concnt*sizeof(angle_struct)) ; }
  else
  { ineqs = *p_ineqs ; }
  if (*p_noload == NULL)
  { noload = (int *) MALLOC((orig_sys->concnt+1)*sizeof(int)) ; }
  else
  { noload = *p_noload ; }
# ifdef DYLP_PARANOIA
  if (ineqs->angles == NULL)
  { errmsg(2,rtnnme,"angle array") ;
    return (FALSE) ; }
# endif
  angles = ineqs->angles ;

  m = orig_sys->concnt ;
  n = orig_sys->varcnt ;
  retval = TRUE ;
/*
  Scan orig_sys, sorting the constraints into eqs and ineqs->angles and
  discarding empty and nonbinding constraints and checking range constraints
  for prima facie infeasibility.
*/
  ai = pkvec_new(0) ;
  eqcnt = 0 ;
  ineqcnt = 0 ;
  noloadcnt = 0 ;
  ctyp = orig_sys->ctyp ;
  rhs = orig_sys->rhs ;
  rhslow = orig_sys->rhslow ;

  for (i = 1 ; i <= m ; i++)
  { if (ctyp[i] == contypNB)
    { dy_lp->sys.cons.unloadable++ ;
      dy_lp->sys.cons.loadable-- ;
      noload[++noloadcnt] = i ;
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.setup >= 4)
      { dyio_outfmt(dy_logchn,dy_gtxecho,
		    "\n\tskipping nonbinding constraint %s (%d).",
		    consys_nme(orig_sys,'c',i,FALSE,NULL),i) ; }
#     endif
      continue ; }
    if (ctyp[i] == contypRNG)
    { if (rhs[i] < rhslow[i])
      { 
#	ifndef DYLP_NDEBUG
	if (dy_opts->print.setup >= 1)
	{ dyio_outfmt(dy_logchn,dy_gtxecho,
		      "\n  Prima facie infeasibility for %s (%d),",
		      consys_nme(orig_sys,'c',i,FALSE,NULL),i) ;
	  dyio_outfmt(dy_logchn,dy_gtxecho," rhslow = %g > rhs = %g.",
		      rhslow[i],rhs[i]) ; }
#	endif
	dy_lp->lpret = lpINFEAS ; } }
    if (consys_getrow_pk(orig_sys,i,&ai) == FALSE)
    { errmsg(122,rtnnme,orig_sys->nme,"row",
	     consys_nme(orig_sys,'c',i,TRUE,NULL),i) ;
      retval = FALSE ;
      break ; }
    if (ai->cnt == 0)
    { dy_lp->sys.cons.unloadable++ ;
      dy_lp->sys.cons.loadable-- ;
      noload[++noloadcnt] = i ;
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.setup >= 4)
      { dyio_outfmt(dy_logchn,dy_gtxecho,
		    "\n\tskipping empty constraint %s (%d).",
		    consys_nme(orig_sys,'c',i,FALSE,NULL),i) ; }
#     endif
    }
    else
    { if (ctyp[i] == contypEQ)
      { eqs[++eqcnt] = i ; }
      else
      { angles[ineqcnt++].ndx = i ; } } }

  if (ai != NULL) pkvec_free(ai) ;
  eqs[0] = eqcnt ;
  ineqs->cnt = ineqcnt ;
  noload[0] = noloadcnt ;

# ifndef DYLP_NDEBUG
  if (dy_opts->print.setup >= 2)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
		"\n    found %d equalities, %d inequalities",
	        eqcnt,ineqcnt) ;
    if (noloadcnt != 0)
    { dyio_outfmt(dy_logchn,dy_gtxecho,", discarded %d empty/nonbinding rows",
		  noloadcnt) ; }
    dyio_outchr(dy_logchn,dy_gtxecho,'.') ; }
# endif
# ifdef DYLP_PARANOIA
  if (dy_lp->sys.cons.loadable+dy_lp->sys.cons.unloadable != orig_sys->concnt)
  { errmsg(1,rtnnme,__LINE__) ;
    retval = FALSE ; }
# endif

  if (retval == FALSE)
  { if (*p_eqs == NULL && eqs != NULL) FREE(eqs) ;
    if (*p_ineqs == NULL && ineqs != NULL)
    { if (ineqs->angles != NULL) FREE(angles) ;
      FREE(ineqs) ; }
    if (*p_noload == NULL && noload != NULL) FREE(noload) ;
    return (FALSE) ; }
/*
  Now, how much more work do we need to do? Decide if we need to calculate
  angles.
*/
  need_angles = TRUE ;
  if (dy_opts->fullsys == TRUE) need_angles = FALSE ;
# ifdef DYLP_STATISTICS
  need_angles = TRUE ;
# endif
# ifndef DYLP_NDEBUG
  if (dy_opts->print.setup >= 2) need_angles = TRUE ;
# endif
/*
  If we need the angles, get down to it. Calculate
      angle(a<i>,c) = (180/pi)*arccos(dot(a<i>,c)/(||a<i>||||c||))
  for each inequality, then sort the list of inequalities.
*/
  if (need_angles == TRUE)
  { c = orig_sys->obj ;
    nearcnt = 0 ;
    perpcnt = 0 ;
    farcnt = 0 ;
    cnorm = exvec_2norm(c,n) ;
    pi180 = 180/acos(-1.0) ;
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.setup >= 4)
    { dyio_outfmt(dy_logchn,dy_gtxecho,"\n\t||c|| = %.4f",cnorm) ;
      dyio_outfmt(dy_logchn,dy_gtxecho,
	     "\n\tConstraint\t\t||a||\t(c/||c||).(a/||a||)\tangle") ; }
#   endif

    for (ndx = 0 ; ndx < ineqcnt ; ndx++)
    { i = angles[ndx].ndx ;
      aidotc = consys_dotrow(orig_sys,i,c) ;
      if (withintol(aidotc,0,dy_tols->zero))
      { 
#	ifndef DYLP_NDEBUG
	ainorm = consys_2normrow(orig_sys,i) ;
#	endif
	anglei = 90 ; }
      else
      { ainorm = consys_2normrow(orig_sys,i) ;
	anglei = pi180*acos(aidotc/(ainorm*cnorm)) ; }
      angles[ndx].angle = anglei ;
      if (anglei > 90)
      { nearcnt++ ; }
      else
      if (anglei < 90)
      { farcnt++ ; }
      else
      { perpcnt++ ; }
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.setup >= 4)
      { dyio_outfmt(dy_logchn,dy_gtxecho,
		    "\n\t%-10s (%3d) %12.4f %18.10f%15.6f",
		    consys_nme(orig_sys,'c',i,FALSE,NULL),i,
		    ainorm,aidotc/(ainorm*cnorm),anglei) ; }
#     endif
    }
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.setup >= 2)
    { dyio_outfmt(dy_logchn,dy_gtxecho,
		  "\n    inequality partition %d near, %d perp, %d far.",
		  nearcnt,perpcnt,farcnt) ; }
#   endif
    ineqs->perp = nearcnt ;
    ineqs->far = nearcnt+perpcnt ;
    qsort(&angles[0],ineqcnt,sizeof(angle_struct),near_perp_far) ; }

# ifdef DYLP_STATISTICS
/*
  If we're collecting full statistics, record the angle for each inequality
  and calculate a histogram (36 5 degree bins and a dedicated 90 degree
  bin).  This could have been folded into the previous loop, but it's more
  clear as a separate loop. Bins below 90 degrees include their lower bound,
  bins above 90 degrees include their upper bound.
*/
  if (dy_stats != NULL)
  { for (ndx = 0 ; ndx < ineqcnt ; ndx++)
    { i = angles[ndx].ndx ;
      anglei = angles[ndx].angle ;
      dy_stats->cons.angle[i] = anglei ;
      if (anglei == 90)
      { dy_stats->angle.hist[DYSTATS_HISTBINS-1]++ ; }
      else
      { if (anglei < 90)
	{ k = (int) floor(anglei/5) ; }
	else
	{ k = (int) ceil(anglei/5)-1 ; }
	k = minn(k,(DYSTATS_HISTBINS-1)) ;
	dy_stats->angle.hist[k]++ ; }
      if (anglei > dy_stats->angle.max)
      { dy_stats->angle.max = (float) anglei ; }
      if (anglei < dy_stats->angle.min)
      { dy_stats->angle.min = (float) anglei ; } } }
# endif

/*
  Set up the return values and we're done.
*/
  if (*p_eqs == NULL) *p_eqs = eqs ;
  if (*p_ineqs == NULL) *p_ineqs = ineqs ;
  if (*p_noload == NULL) *p_noload = noload ;

  return (TRUE) ; }



static bool cold_createdysys (consys_struct *orig_sys, int eqcnt, int ineqcnt)
/*
  This routine creates the dy_sys constraint system and attaches the
  translation vectors that are used to move between the original (orig_sys)
  and active (dy_sys) constraint systems.
  
  While unrelated to creating dy_sys, it turns out that this routine is a
  convenient place for another pre-loading activity, grooming the bounds on
  variables in orig_sys and assigning an initial status to the inactive
  variables. Where lb<j> > ub<j>, indicate prima facie infeasibility.
  Where |ub<j> - lb<j>| < dy_tols->pfeas, the bounds are forced to
  equality and the variable is marked as fixed (NBFX). (This sort of thing
  can happen when a client program is tweaking the bounds on variables.) This
  may or may not be visible to the client --- depends on whether dylp has
  made a local copy of the constraint system.  Variables which are not fixed
  are arbitrarily set to the bound that best matches their objective
  coefficient.

  What's a good estimate of maximum size? If the fullsys option is specified,
  just use the original constraint system size. Otherwise, take the attitude
  that some (user) specified fraction of the original constraints and
  variables (active.cons, active.vars) could be active at any one time while
  we're solving the LP. This is purely for efficiency -- if we're undersized,
  the system will automatically expand.

  Given that all equalities are loaded, initcons.frac is interpreted as the
  initial fraction of inequalities that should be loaded. Again, for
  efficiency active.cons should be greater than initcons.frac. dy_sys is
  created with logicals enabled.

  Parameters:
    orig_sys:	the original constraint system
    eqcnt:	the number of equalities in orig_sys
    ineqcnt:	the number of inequalities in orig_sys

  Returns: TRUE if the system is created without error, FALSE otherwise.
*/

{ int m_sze,n_sze ;

  char nmebuf[50] ;

  flags parts = CONSYS_OBJ|CONSYS_VUB|CONSYS_VLB|CONSYS_RHS|CONSYS_RHSLOW|
		CONSYS_VTYP|CONSYS_CTYP,
	opts = CONSYS_LVARS|CONSYS_WRNATT ;

  const char *rtnnme = "cold_createdysys" ;

/*
  Settle the appropriate size for dy_sys, and create the consys structure.
  The variable capacity is set to accommodate logicals.
*/
  dyio_outfxd(nmebuf,-((int) (sizeof(nmebuf)-1)),
	      'l',"%s[actv]",orig_sys->nme) ;
  if (dy_opts->fullsys == TRUE)
  { m_sze = orig_sys->concnt ;
    n_sze = orig_sys->archvcnt+m_sze ; }
  else
  { m_sze = eqcnt+((int) (orig_sys->concnt*dy_opts->active.cons)) ;
    n_sze = (int) (orig_sys->archvcnt*dy_opts->active.vars+m_sze) ; }
# ifndef DYLP_NDEBUG
  if (dy_opts->print.setup >= 1)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
		"\n    creating constraint system %s (%d x %d)",
	        nmebuf,m_sze,n_sze) ; }
# endif
  dy_sys = consys_create(nmebuf,parts,opts,m_sze,n_sze,dy_tols->inf) ;
  if (dy_sys == NULL)
  { errmsg(152,rtnnme,nmebuf) ;
    return (FALSE) ; }
/*
  Hang a set of translation vectors onto each system: origcons and origvars
  on orig_sys, and actcons and actvars on dy_sys. dy_origvars is cleared to
  0 as it's attached, and this is important to indicate that the original
  variables have no predefined status.

  Note that we won't leak memory here, even if there's an error. Any attached
  vectors will be deleted when dy_sys is deleted.
*/
  dy_actvars = NULL ;
  if (consys_attach(dy_sys,CONSYS_ROW,
		    sizeof(int),(void **) &dy_actvars) == FALSE)
  { errmsg(100,rtnnme,dy_sys->nme,"active -> original variable map") ;
    return (FALSE) ; }
  dy_actcons = NULL ;
  if (consys_attach(dy_sys,CONSYS_COL,
		    sizeof(int),(void **) &dy_actcons) == FALSE)
  { errmsg(100,rtnnme,dy_sys->nme,"active -> original constraint map") ;
    return (FALSE) ; }
  dy_origvars = NULL ;
  if (consys_attach(orig_sys,CONSYS_ROW,
		    sizeof(int),(void **) &dy_origvars) == FALSE)
  { errmsg(100,rtnnme,orig_sys->nme,"original -> active variable map") ;
    return (FALSE) ; }
  dy_origcons = NULL ;
  if (consys_attach(orig_sys,CONSYS_COL,
		    sizeof(int),(void **) &dy_origcons) == FALSE)
  { errmsg(100,rtnnme,orig_sys->nme,"original -> active constraint map") ;
    return (FALSE) ; }

  return (TRUE) ; }


static bool cold_scanvars (consys_struct *orig_sys)
/*
  Assign a status to all variables in orig_sys. There are lots of reasons
  for doing this right off the top:

    * We can use origvars == 0 as a paranoid check from here on out.

    * It's not uncommon for people to fix variables in an MPS model. This takes
      them out of consideration right from the start. Assign status NBFX and
      mark them as ineligible for activation.

    * dylp's intended use is in branch-and-cut codes, which tend to update
      bounds on a regular basis. This can result in bounds which are not
      precisely equal, but are within the primal feasibility tolerance of one
      another. The variable is fixed, for all practical purposes, but it's a
      `dirty' fix. Push it to one bound or another and mark it ineligible for
      activation.

    * A careless client may hand us a model with prima facie infeasibility,
      lb<j> > ub<j>. We need to detect this and report it --- the main body
      of dylp does not react well to this situation.

    * More legitimate (but not by much), the client may hand us a model with
      variables that are only present in the objective. We can set these to
      bound based on the objective coefficient and mark them as ineligible
      for activation. If the appropriate bound is not finite, we can report
      unboundedness.

    * If there are no variables with upper and lower bounds (`flippable', in
      dual multipivot) then we might as well turn multipivot off.

  Parameters:
    orig_sys:	the original constraint system

  Returns: TRUE if all goes well, FALSE otherwise.
	   If prima facie infeasibility is detected, dy_lp->lpret is set to
	   lpINFEAS.
	   If prima facie unboundedness is detected, dy_lp->lpret is set to
	   lpUNBOUNDED.
*/

{ int j,flippable ;
  flags statj ;
  double *vlb,*vub,*obj ;
  double lj,uj,cj ;

  bool retval ;

  pkvec_struct *aj ;

  char *rtnnme = "cold_scanvars" ;

  retval = TRUE ;
/*
  Prep, then open up a loop to do the scan.
*/
  vlb = orig_sys->vlb ;
  vub = orig_sys->vub ;
  obj = orig_sys->obj ;
  flippable = 0 ;

  aj = pkvec_new(0) ;

  for (j = 1 ; j <= orig_sys->varcnt ; j++)
  { lj = vlb[j] ;
    uj = vub[j] ;
    cj = obj[j] ;
/*
  Sort into cases. Variables with two finite bounds first. If we have a
  `dirty' fix, set the variable to one of the bounds and mark it ineligible
  for activation. If the bounds are crossed, note prima facie infeasibility.
  Otherwise, note that the variable is flippable and assign a status based on
  the objective coefficient.
*/
    if (lj > -dy_tols->inf && uj < dy_tols->inf)
    { if (atbnd(lj,uj) && lj != uj)
      { if (cj < 0)
	{ statj = vstatNBUB|vstatNOLOAD ; }
	else
	{ statj = vstatNBLB|vstatNOLOAD ; }
#	ifndef DYLP_NDEBUG
	if (dy_opts->print.setup >= 3)
	{ dyio_outfmt(dy_logchn,dy_gtxecho,"\n\tDirty fixed variable %s (%d)",
		      consys_nme(orig_sys,'v',j,0,0),j) ;
	  dyio_outfmt(dy_logchn,dy_gtxecho,
		      " assigned status %s.",dy_prtvstat(statj)) ;
	  dyio_outfmt(dy_logchn,dy_gtxecho,
		      "\n\t  original lb = %g, ub = %g, diff = %g, tol = %g",
		      lj,uj,uj-lj,dy_tols->pfeas) ; }
#	endif
      }
      else
      if (uj < lj)
      { dy_lp->lpret = lpINFEAS ;
#       ifndef DYLP_NDEBUG
	if (dy_opts->print.setup >= 1)
	{ dyio_outfmt(dy_logchn,dy_gtxecho,
	      "\n\tPrima facie infeasibility for %s (%d), lb = %g > ub = %g.",
	      consys_nme(orig_sys,'v',j,0,0),j,lj,uj) ; }
#       endif
	statj = vstatNBLB|vstatNOLOAD ; }
      else
      if (lj == uj)
      { statj = vstatNBFX|vstatNOLOAD ;
#	ifndef DYLP_NDEBUG
	if (dy_opts->print.setup >= 3)
	{ dyio_outfmt(dy_logchn,dy_gtxecho,"\n\tFixed variable %s (%d) = %g",
		      consys_nme(orig_sys,'v',j,0,0),j,lj) ;
	  dyio_outfmt(dy_logchn,dy_gtxecho,
		      " assigned status %s.",dy_prtvstat(statj)) ; }
#	endif
      }
      else
      { flippable++ ;
	if (cj < 0)
	{ statj = vstatNBUB ; }
	else
	{ statj = vstatNBLB ; } } }
/*
  For variables with one finite bound, choose the finite bound as the initial
  status. If there are no finite bounds, we're left with NBFR.
*/
    else
    if (lj > -dy_tols->inf)
    { statj = vstatNBLB ; }
    else
    if (uj < dy_tols->inf)
    { statj = vstatNBUB ; }
    else
    { statj = vstatNBFR ; }
/*
  Check if the variable is not referenced by any constraint. Note that this
  isn't the definitive test. There's a variation of this pathology in which
  the client has left a nonbinding row in the system and the variable is
  referenced only by the nonbinding row.  Doing the necessary work here fails
  the `make the common case fast' test. Doing the necessary work later is
  still questionable. See comments in dy_varmgmt::scanPrimVarStdAct.

  If we pass the test for boundedness, the status is already correct.  Just
  mark it as ineligible for activation.
*/
    if (consys_getcol_pk(orig_sys,j,&aj) == FALSE)
    { errmsg(122,rtnnme,dy_sys->nme,"column",
	     consys_nme(dy_sys,'v',j,TRUE,NULL),j) ;
      retval = FALSE ;
      break ; }
    if (aj->cnt == 0)
    { if ((cj > 0 && lj <= -dy_tols->inf) || (cj < 0 && uj >= dy_tols->inf))
      { dy_lp->lpret = lpUNBOUNDED ;
	dy_lp->z = -j ;
#       ifndef DYLP_NDEBUG
	if (dy_opts->print.setup >= 1)
	{ dyio_outfmt(dy_logchn,dy_gtxecho,
		      "\n\tPrima facie unboundedness for %s (%d), ",
		      consys_nme(orig_sys,'v',j,0,0),j) ;
	  dyio_outfmt(dy_logchn,dy_gtxecho,"c = %g, lb = %g, ub = %g.",
		      cj,lj,uj) ; }
#       endif
      }
      else
      { setflg(statj,vstatNOLOAD) ;
#	ifndef DYLP_NDEBUG
	if (dy_opts->print.setup >= 3)
	{ dyio_outfmt(dy_logchn,dy_gtxecho,
		      "\n\tEmpty column for variable %s (%d),",
		      consys_nme(orig_sys,'v',j,0,0),j,lj) ;
	  dyio_outfmt(dy_logchn,dy_gtxecho,
		      " assigned status %s,",dy_prtvstat(statj)) ;
	  dyio_outfmt(dy_logchn,dy_gtxecho,
		      " c = %g, lb = %g, ub = %g.",cj,lj,uj) ; }
#	endif
      } }

    dy_origvars[j] = -((int) statj) ;
    if (flgon(statj,vstatNOLOAD) == TRUE)
    { dy_lp->sys.vars.unloadable++ ;
      dy_lp->sys.vars.loadable-- ; }
  }

  if (aj != NULL) pkvec_free(aj) ;
/*
  That's it for the scan loop.

  Disable dual multipivoting? Fairly arbitrarily, give 25% flippable variables
  as the threshold.
*/
# ifdef DYLP_STATISTICS
  if (dy_stats != NULL) dy_stats->dmulti.flippable = flippable ;
# endif
  uj = ((double) flippable)/orig_sys->varcnt ;
  if (uj < .25 && dy_opts->dpsel.flex == TRUE)
  { dy_opts->dpsel.flex = FALSE ;
    dy_opts->dpsel.strat = 0 ;
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.setup >= 2 || dy_opts->print.dual >= 2)
    { dyio_outfmt(dy_logchn,dy_gtxecho,
	    "\n    %d (%.1f%%) flippable variables; disabling dual multipivot.",
	    flippable,uj*100) ; }
#   endif
  }
# ifdef DYLP_PARANOIA
  if (dy_lp->sys.vars.unloadable+dy_lp->sys.vars.loadable != orig_sys->varcnt)
  { errmsg(1,rtnnme,__LINE__) ;
    retval = FALSE ; }
# endif

  return (retval) ; }





static bool cold_loadfull (consys_struct *orig_sys,
			    int *eqs, ineq_struct *ineqs)
/*
  This routine loads the full original constraint system into dy_sys. It's
  just easier to handle this separate from the more complicated business of
  loading a partial system.

  Parameters:
    orig_sys:	the original constraint system
    eqs:	indices of equality constraints; eq[0] contains the number of
		valid entries
    ineqs:	information on inequalities; only the number of inequalities
		and their indices are used here.

  Returns: TRUE if all constraints are loaded successfully, FALSE otherwise.
*/

{ int i,ndx,eqcnt,ineqcnt ;
  bool retval ;

  angle_struct *angles ;

  const char *rtnnme = "cold_loadfull" ;

  eqcnt = eqs[0] ;
  ineqcnt = ineqs->cnt ;
  retval = TRUE ;
/*
  Load any equalities first.
*/
  if (eqcnt > 0)
  { for (ndx = 1 ; ndx <= eqcnt ; ndx++)
    { i = eqs[ndx] ;
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.setup >= 3)
      { dyio_outfmt(dy_logchn,dy_gtxecho, "\n    activating %s %s (%d) ...",
		    consys_prtcontyp(orig_sys->ctyp[i]),
		    consys_nme(orig_sys,'c',i,FALSE,NULL),i) ; }
#     endif
#     ifdef DYLP_STATISTICS
      if (dy_stats != NULL) dy_stats->cons.init[i] = TRUE ;
#     endif
      if (dy_loadcon(orig_sys,i,TRUE,NULL) == FALSE)
      { errmsg(430,rtnnme,
	       dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,
	       "activate","constraint",
	       consys_nme(orig_sys,'c',i,TRUE,NULL),i) ;
	retval = FALSE ;
	break ; } }
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.setup >= 2)
    { dyio_outfmt(dy_logchn,dy_gtxecho,
		  "\n    transferred %d equalities ...",eqcnt) ; }
#   endif
    if (retval == FALSE)
    { return (FALSE) ; } }
/*
  And the inequalities.
*/
  if (ineqcnt > 0)
  { angles = ineqs->angles ;
    for (ndx = 0 ; ndx < ineqcnt ; ndx++)
    { i = angles[ndx].ndx ;
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.setup >= 3)
      { dyio_outfmt(dy_logchn,dy_gtxecho, "\n    activating %s %s (%d) ...",
		    consys_prtcontyp(orig_sys->ctyp[i]),
		    consys_nme(orig_sys,'c',i,FALSE,NULL),i) ; }
#     endif
#     ifdef DYLP_STATISTICS
      if (dy_stats != NULL) dy_stats->cons.init[i] = TRUE ;
#     endif
      if (dy_loadcon(orig_sys,i,TRUE,NULL) == FALSE)
      { errmsg(430,rtnnme,
	       dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,
	       "activate","constraint",
	       consys_nme(orig_sys,'c',i,TRUE,NULL),i) ;
	retval = FALSE ;
	break ; } }
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.setup >= 2)
    { dyio_outfmt(dy_logchn,dy_gtxecho,
		  "\n    transferred %d inequalities ...",ineqcnt) ; }
#   endif
  }

  return (retval) ; }



static bool cold_loadpartial (consys_struct *orig_sys,
			      int *eqs, ineq_struct *ineqs)
/*
  This routine loads some fraction of the original constraint system into
  dy_sys. As explained in the comments at the head of the file, the
  inequalities are sorted by angle(a<i>,c). The client can specify one or two
  angular intervals, and the upper and lower bounds of each interval can be
  closed or open. Both intervals are sampled at the specified fractional rate.

  Parameters:
    orig_sys:	the original constraint system
    eqs:	indices of equality constraints; eq[0] contains the number of
		valid entries
    ineqs:	information on inequalities; only the number of inequalities
		and their indices are used here.

  Returns: TRUE if all constraints are loaded successfully, FALSE otherwise.
*/

{ int i,ndx,eqcnt,ineqcnt,ineq_actvcnt ;
  bool retval ;
  angle_struct *angles ;

  int intcnt,intndx ;
  double albs[2],aubs[2],alb,aub,dblndx,incr ;

  const char *rtnnme = "cold_loadpartial" ;

  eqcnt = eqs[0] ;
  ineqcnt = ineqs->cnt ;
  retval = TRUE ;
/*
  Load any equalities first.
*/
  if (eqcnt > 0)
  { for (ndx = 1 ; ndx <= eqcnt ; ndx++)
    { i = eqs[ndx] ;
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.setup >= 3)
      { dyio_outfmt(dy_logchn,dy_gtxecho, "\n    activating %s %s (%d) ...",
		    consys_prtcontyp(orig_sys->ctyp[i]),
		    consys_nme(orig_sys,'c',i,FALSE,NULL),i) ; }
#     endif
#     ifdef DYLP_STATISTICS
      if (dy_stats != NULL) dy_stats->cons.init[i] = TRUE ;
#     endif
      if (dy_loadcon(orig_sys,i,TRUE,NULL) == FALSE)
      { errmsg(430,rtnnme,
	       dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,
	       "activate","constraint",
	       consys_nme(orig_sys,'c',i,TRUE,NULL),i) ;
	retval = FALSE ;
	break ; } }
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.setup >= 2)
    { dyio_outfmt(dy_logchn,dy_gtxecho,
		  "\n    transferred %d equalities ...",eqcnt) ; }
#   endif
    if (retval == FALSE)
    { return (FALSE) ; } }
/*
  If there are no inequalities, we're done.
*/
  if (ineqcnt == 0) return (TRUE) ;
/*
  Establish the number of intervals, angular boundaries, and fractional
  increment. Recall that constraints are sorted nonincreasing from 180
  degrees to 0 degrees. To approximate an open interval, we pull in the
  specified bound by the zero tolerance; for a closed interval, we push it
  out.
*/
# ifndef DYLP_NDEBUG
  if (dy_opts->print.setup >= 2)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
	        "\n    inequalities: sampling %.2f, %c %.4f %.4f %c",
	        dy_opts->initcons.frac,
	        (dy_opts->initcons.i1uopen == TRUE)?'(':'[',
	        dy_opts->initcons.i1u,dy_opts->initcons.i1l,
	        (dy_opts->initcons.i1lopen == TRUE)?')':']') ;
    if (dy_opts->initcons.i2valid == TRUE)
    { dyio_outfmt(dy_logchn,dy_gtxecho,", %c %.4f %.4f %c",
		  (dy_opts->initcons.i2uopen == TRUE)?'(':'[',
		  dy_opts->initcons.i2u,dy_opts->initcons.i2l,
		  (dy_opts->initcons.i2lopen == TRUE)?')':']') ; } }
# endif
  incr = 1/dy_opts->initcons.frac ;
  if (dy_opts->initcons.i1lopen == TRUE)
  { albs[0] = dy_opts->initcons.i1l+dy_tols->zero ; }
  else
  { albs[0] = dy_opts->initcons.i1l-dy_tols->zero ; }
  if (dy_opts->initcons.i1uopen == TRUE)
  { aubs[0] = dy_opts->initcons.i1u-dy_tols->zero ; }
  else
  { aubs[0] = dy_opts->initcons.i1u+dy_tols->zero ; }
  if (dy_opts->initcons.i2valid == TRUE)
  { if (dy_opts->initcons.i2lopen == TRUE)
    { albs[1] = dy_opts->initcons.i2l+dy_tols->zero ; }
    else
    { albs[1] = dy_opts->initcons.i2l-dy_tols->zero ; }
    if (dy_opts->initcons.i2uopen == TRUE)
    { aubs[1] = dy_opts->initcons.i2u-dy_tols->zero ; }
    else
    { aubs[1] = dy_opts->initcons.i2u+dy_tols->zero ; }
    intcnt = 2 ; }
  else
  { intcnt = 1 ; }
  angles = ineqs->angles ;
  ineq_actvcnt = 0 ;
/*
  The outer loop is just to avoid replicating code for each interval.
*/
  ndx = 0 ;
  for (intndx = 0 ; intndx < intcnt ; intndx++)
  { alb = albs[intndx] ;
    aub = aubs[intndx] ;
/*
  Scan to the upper boundary of the interval. Consider the possibility that
  we might run off the end of angles. An easy and all too common example: the
  intervals exclude 90 degrees, and all inequalities are at 90 degrees to the
  objective. This is an artifact of a modelling style in which the objective
  function is a single variable z, constrained to be equal to the real
  objective cx. Then z appears only in the single equality z - cx = 0.
*/
    while (ndx < ineqcnt && angles[ndx].angle > aub) ndx++ ;
    dblndx = ndx ;
/*
  Now the actual loading loop. We load the constraint, then move ndx to the
  next value according to the sampling fraction.
*/
    while (ndx < ineqcnt && angles[ndx].angle > alb) 
    { i = angles[ndx].ndx ;
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.setup >= 3)
      { dyio_outfmt(dy_logchn,dy_gtxecho,
		    "\n    activating %s %s (%d) %g off 90 ...",
		    consys_prtcontyp(orig_sys->ctyp[i]),
		    consys_nme(orig_sys,'c',i,FALSE,NULL),i,
		    (angles[ndx].angle-90)) ; }
#     endif
#     ifdef DYLP_STATISTICS
      if (dy_stats != NULL) dy_stats->cons.init[i] = TRUE ;
#     endif
      if (dy_loadcon(orig_sys,i,TRUE,NULL) == FALSE)
      { errmsg(430,rtnnme,
	       dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,
	       "activate","constraint",
	       consys_nme(orig_sys,'c',i,TRUE,NULL),i) ;
	retval = FALSE ;
	break ; }
      ineq_actvcnt++ ;
      dblndx += incr ;
      ndx = (int) dblndx ; } }
/*
  That's it. A little information and we're out of here.
*/
# ifndef DYLP_NDEBUG
  if (dy_opts->print.setup >= 2)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
	        "\n    transferred %d inequalities ...",ineq_actvcnt) ; }
# endif

  return (retval) ; }



dyret_enum dy_coldstart (consys_struct *orig_sys)

/*
  This routine is responsible for setting up the lp problem that'll be solved
  by dylp. It creates the lp and constraint system structures (dy_lp and
  dy_sys, respectively) and oversees the load of the initial constraint system.

  Parameters:
    orig_sys:	The original constraint system

  Returns: dyrOK if the setup completes without error, dyrFATAL otherwise.
*/

{ int i,j,ndx,n,eqcnt,ineqcnt,noloadcnt ;
  double *vlb,*vub,*obj ;
  double vlbj,vubj,objj ; 
  flags statj ;
  int *eqs,*noload ;
  ineq_struct *ineqs ;

# ifndef DYLP_NDEBUG
  int nbfxcnt = 0 ;
# endif

  bool retval ;
  
  const char *rtnnme = "dy_coldstart" ;

/*
  Initialise the statistics on loadable/unloadable variables and constraints.
*/
  dy_lp->sys.forcedfull = FALSE ;
  dy_lp->sys.vars.loadable = orig_sys->varcnt ;
  dy_lp->sys.vars.unloadable = 0 ;
  dy_lp->sys.cons.loadable = orig_sys->concnt ;
  dy_lp->sys.cons.unloadable = 0 ;
/*
  To get started, sort the constraints into equalities and inequalities. If
  we're loading a partial active system, the inequalities will also be sorted
  by angle from the objective function.
*/
  eqs = NULL ;
  ineqs = NULL ;
  noload = NULL ;
  retval = cold_sortcons(orig_sys,&eqs,&ineqs,&noload) ;
  if (retval == FALSE)
  { errmsg(312,rtnnme,
	  orig_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters) ;
    if (eqs != NULL) FREE(eqs) ;
    if (ineqs != NULL)
    { if (ineqs->angles != NULL) FREE(ineqs->angles) ;
      FREE(ineqs) ; }
    if (noload != NULL) FREE(noload) ;
    return (dyrFATAL) ; }
  eqcnt = eqs[0] ;
  ineqcnt = ineqs->cnt ;
/*
  Next, create dy_sys and attach the translation vectors that will allow us to
  move between the original and active systems.
*/
  retval = cold_createdysys(orig_sys,eqcnt,ineqcnt) ;
  if (retval == FALSE)
  { if (eqs != NULL) FREE(eqs) ;
    if (ineqs != NULL)
    { if (ineqs->angles != NULL) FREE(ineqs->angles) ;
      FREE(ineqs) ; }
    if (noload != NULL) FREE(noload) ;
    return (dyrFATAL) ; }
/*
  Do an initial scan of the variables. Assign initial status, deal with a few
  pathological cases, and detect prima facie infeasibility and unboundedness.
*/
  retval = cold_scanvars(orig_sys) ;
  if (retval == FALSE)
  { if (eqs != NULL) FREE(eqs) ;
    if (ineqs != NULL)
    { if (ineqs->angles != NULL) FREE(ineqs->angles) ;
      FREE(ineqs) ; }
    if (noload != NULL) FREE(noload) ;
    return (dyrFATAL) ; }
/*
  Mark the constraints that are ineligible for loading and dispose of noload.
*/
  noloadcnt = noload[0] ;
  if (noloadcnt > 0)
  { for (ndx = 1 ; ndx <= noloadcnt ; ndx++)
    { i = noload[ndx] ;
      MARK_UNLOADABLE_CON(i) ; } }
  if (noload != NULL) FREE(noload) ;
/*
  Transfer the required constraints. Once this is done,  we're finished with
  the lists of equalities and inequalities.

  There's a pathology we may need to deal with here. It can happen that the
  client is negligent and passes in a constraint system which contains
  nonbinding rows and variables which are referenced only in the nonbinding
  rows. If we're running in dynamic mode, these variables will be found during
  normal variable activation. But if we're running in fullsys mode, there will
  be no second chance. Check for the situation and deal with it.
*/
  if (dy_opts->fullsys == TRUE)
  { retval = cold_loadfull(orig_sys,eqs,ineqs) ;
    dy_lp->sys.forcedfull = TRUE ;
    if (retval == TRUE &&
	dy_sys->archvcnt+dy_lp->sys.vars.unloadable < orig_sys->varcnt)
    { 
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.setup >= 1)
      { dyio_outfmt(dy_logchn,dy_gtxecho,
	  "\n  system %s has %d variables referenced only by nonbinding rows.",
	  orig_sys->nme,
	  (orig_sys->varcnt-(dy_sys->archvcnt+dy_lp->sys.vars.unloadable))) ; }
#     endif
    } }
  else
  { retval = cold_loadpartial(orig_sys,eqs,ineqs) ;
    dy_lp->sys.forcedfull = FALSE ; }
  if (eqs != NULL) FREE(eqs) ;
  if (ineqs != NULL)
  { if (ineqs->angles != NULL) FREE(ineqs->angles) ;
    FREE(ineqs) ; }
  if (retval == FALSE)
  { errmsg(313,rtnnme,
	    dy_sys->nme,dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,
	    orig_sys->nme) ;
    return (dyrFATAL) ; }
/*
  Scan the variables in orig_sys once again and calculate the correction to
  the objective function. (loadcon handled any rhs adjustments for the
  constraints it loaded.)  Inactive free variables are assumed to be 0.
*/
  n = orig_sys->varcnt ;
  vlb = orig_sys->vlb ;
  vub = orig_sys->vub ;
  obj = orig_sys->obj ;
  dy_lp->inactzcorr = 0 ;
  for (j = 1 ; j <= n ; j++)
  { if (INACTIVE_VAR(j))
    { statj = (flags) (-dy_origvars[j]) ;
      clrflg(statj,vstatNOLOAD) ;
      vlbj = vlb[j] ;
      vubj = vub[j] ;
      objj = obj[j] ;
      switch (statj)
      { case vstatNBLB:
	case vstatNBFX:
	{ dy_lp->inactzcorr += objj*vlbj ;
#	  ifndef DYLP_NDEBUG
	  if (statj == vstatNBFX) nbfxcnt++ ;
#	  endif
	  break ; }
	case vstatNBUB:
	{ dy_lp->inactzcorr += objj*vubj ;
	  break ; } } } }
/*
  Paranoid checks and informational print statements. Apologies for abusing
  previously declared variables for the print loop.
*/
# ifdef DYLP_PARANOIA
  if (dy_chkdysys(orig_sys) == FALSE) return (dyrFATAL) ;
# endif
# ifndef DYLP_NDEBUG
  if (dy_opts->print.setup >= 1)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
	        "\n  system %s has %d constraints, (%d+%d) variables",
	        dy_sys->nme,dy_sys->concnt,dy_sys->archvcnt,dy_sys->logvcnt) ;
    dyio_outfmt(dy_logchn,dy_gtxecho,"\n  %d constraints, %d variables ",
	        orig_sys->concnt-dy_sys->concnt,
	        orig_sys->archvcnt-dy_sys->archvcnt) ;
    if (nbfxcnt > 0) dyio_outfmt(dy_logchn,dy_gtxecho," (%d fixed)",nbfxcnt) ;
    dyio_outfmt(dy_logchn,dy_gtxecho,
		" remain inactive in system %s.",orig_sys->nme) ;
    if (dy_lp->sys.cons.unloadable > 0 || dy_lp->sys.vars.unloadable > 0)
    { dyio_outfmt(dy_logchn,dy_gtxecho,
		  "\n  %d constraints, %d variables unloadable.",
		  dy_lp->sys.cons.unloadable,dy_lp->sys.vars.unloadable) ; }
    if (dy_opts->print.setup >= 6)
    { vubj = 0 ;
      for (j = 1 ; j <= n ; j++)
      { if (INACTIVE_VAR(j))
	{ statj = (flags)(-dy_origvars[j]) ;
	  switch (getflg(statj,vstatSTATUS))
	  { case vstatNBUB:
	    { vlbj = vub[j] ;
	      break ; }
	    case vstatNBLB:
	    { vlbj = vlb[j] ;
	      break ; }
	    case vstatNBFX:
	    case vstatNBFR:
	    { vlbj = 0 ;
	      break ; }
	    default:
	    { errmsg(433,rtnnme,dy_sys->nme,
		     dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,
		     "inactive",consys_nme(orig_sys,'v',j,TRUE,NULL),
		     j,dy_prtvstat(statj)) ;
	      return (dyrFATAL) ; } }
	  if (vlbj != 0)
	  { if (vubj == 0)
	    { dyio_outfmt(dy_logchn,dy_gtxecho,
			  "\n\tinactive variables with nonzero values:") ; }
	    vubj++ ;
	    dyio_outfmt(dy_logchn,dy_gtxecho,"\n\t%s (%d) = %g, status %s",
		        consys_nme(orig_sys,'v',j,FALSE,NULL),j,vlbj,
		        dy_prtvstat(statj)) ; } } }
      if (vubj == 0)
      { dyio_outfmt(dy_logchn,dy_gtxecho,
		    "\n\tall inactive variables are zero.") ; } } }
# endif

  return (dyrOK) ; }




/*
  Routines to select an initial basis.

  dylp offers three types of starting basis:

    logical:    Basis positions occupied by inequalities are covered with
		slack variables.  Basis positions occupied by equalities are
		covered with artificial variables.  This is the standard
		starting basis.

    slack:	Basis positions occupied by inequalities are covered with
		slack variables. For basis positions occupied by equalities,
		an attempt is made to select architectural variables before
		falling back on artificial variables.

    architectural: For all basis positions, an attempt is made to cover the
		position with an architectural variable, before falling back on
		a slack or artificial. The selection process is performed in a
		way that will use an architectural to cover an equality before
		considering inequalities.
  
  As Bixby notes in an ORSA JOC article, crash is often taken to mean trying
  to guess the optimal basis. That's part of what's happening here, but not
  all. Architectural is chosen as the default basis type on the grounds that
  the default algorithms used in dy_coldstart and subsidiary routines are
  trying to choose constraints that will be tight at optimum. So it makes
  sense to build an initial basis with as few basic logicals as possible.
  But just as important, dy_crash and subsidiary routines try to select
  variables with an eye for numerical stability and ease of factorization.

  The algorithms used here to construct the slack and architectural basis
  styles had their roots in Bixby's article:

    Bixby, R., "Implementing the Simplex Method: The Initial Basis",
    ORSA J. on Computing, 4(3), Summer, 1992, pp. 267-284.
*/
/*
  Some utility declarations used in building the initial basis.
  
  The ibrank structure holds information useful in selecting architectural
  variable for inclusion in the initial basis. The function ib_archvcomp is
  used by qsort to sort the list of eligible candidates.

  ibrank_struct

  bndcnt and nonzero are used for the initial sort of the architecturals.

  If the matrix has been scaled, ajmax = 1 by definition for every row and
  column. But if the matrix hasn't been scaled, it's useful to know the
  maximum.

  Field		Definition
  -----		----------
  ndx		The index of the variable
  bndcnt	The number of bounds
  nonzero	The number of nonzero coefficients in the column
  ajmax		The largest value in the column
*/

typedef struct {
  int ndx ;
  int bndcnt ;
  int nonzero ;
  double ajmax ; } ibrank_struct ;


static int ib_archvcomp (const void *elem1, const void *elem2)
/*
  This function is called by qsort to compare architectural variables.
  The comparison is calculated using the following criteria:
    * Number of bounds: free variables (no finite bounds) are best, followed
      by variables with one finite bound, then two finite bounds
    * Number of coefficients in the column: Fewer coefficients are preferable,
      as this will tend to minimize the work during factoring.
*/


{ const ibrank_struct *v1,*v2 ;

  v1 = (const ibrank_struct *) elem1 ;
  v2 = (const ibrank_struct *) elem2 ;

  if (v1->bndcnt < v2->bndcnt)
  { return (-1) ; }
  else
  if (v1->bndcnt > v2->bndcnt)
  { return (1) ; }
  else
  if (v1->nonzero < v2->nonzero)
  { return (-1) ; }
  else
  if (v1->nonzero > v2->nonzero)
  { return (1) ; }
  else
  { return (0) ; } }



static bool ib_archvrank (int *p_cnt, ibrank_struct **p_archvars)

/*
  This routine sorts the architectural variables in preparation for
  constructing the initial basis. It scans the architectural variables,
  compiling information on the eligible candidates in archvars, and returns a
  sorted list.

  The ranking criteria are explained above in the comments for ib_archvcomp
  (the comparison function passed to qsort).

  Parameters:
    p_cnt:	(o) the number of eligible variables returned in p_archvars
    p_archvars:	(i) empty crashrank_struct vector of sufficient size; if NULL,
		    one will be allocated
		(o) list of eligible architectural variables, sorted in order
		    of preference for inclusion in the basis; may be NULL if
		    there are no candidates (extremely unlikely)

  Returns: TRUE if the sort concludes without error, FALSE otherwise
*/

{ int j,m,n,eligible ;
  double *vlb,*vub ;
  double vlbj,vubj ;
  pkvec_struct *aj ;
  ibrank_struct *archvars ;
  bool scaled,retval ;

  const char *rtnnme = "ib_archvrank" ;

# ifndef DYLP_NDEBUG

  int freecnt,onebndcnt,twobndcnt ;

  freecnt = 0 ;
  onebndcnt = 0 ;
  twobndcnt = 0 ;

# endif

  retval = TRUE ;

  vlb = dy_sys->vlb ;
  vub = dy_sys->vub ;

/*
  Did the client supply a vector? If not, allocate one.
*/
  if (*p_archvars == NULL)
  { archvars =
	(ibrank_struct *) MALLOC(dy_sys->archvcnt*sizeof(ibrank_struct)) ; }
  else
  { archvars = *p_archvars ; }
/*
  Scan the architecturals, recording the necessary information in archvars.
  Fixed variables should not be loaded, so we should not see them here. If
  the constraint system has not been scaled, recover the maximum coefficient
  for the column.
*/
  n = dy_sys->varcnt ;
  m = dy_sys->logvcnt ;
  scaled = dy_isscaled() ;
  eligible = 0 ;
  aj = pkvec_new(0) ;
  for (j = m+1 ; j <= n ; j++)
  { if (consys_getcol_pk(dy_sys,j,&aj) == FALSE)
    { errmsg(122,rtnnme,dy_sys->nme,"column",
	     consys_nme(dy_sys,'v',j,TRUE,NULL),j) ;
      retval = FALSE ;
      break ; }
    archvars[eligible].ndx = j ;
    archvars[eligible].nonzero = aj->cnt ;
    if (scaled)
    { archvars[eligible].ajmax = 1.0 ; }
    else
    { archvars[eligible].ajmax = consys_infnormcol(dy_sys,j) ; }
    vlbj = vlb[j] ;
    vubj = vub[j] ;
    if (vlbj > -dy_tols->inf && vubj < dy_tols->inf)
    { archvars[eligible].bndcnt = 2 ;
#     ifdef DYLP_PARANOIA
      if (vlbj == vubj)
      { errmsg(1,rtnnme,__LINE__) ;
	retval = FALSE ;
	break ; }
#     endif
#     ifndef DYLP_NDEBUG
      twobndcnt++ ;
#     endif
    }
    else
    if (vlbj > -dy_tols->inf || vubj < dy_tols->inf)
    { archvars[eligible].bndcnt = 1 ;
#     ifndef DYLP_NDEBUG
      onebndcnt++ ;
#     endif
    }
    else
    { archvars[eligible].bndcnt = 0 ;
#     ifndef DYLP_NDEBUG
      freecnt++ ;
#     endif
    }
    eligible++ ; }

# ifndef DYLP_NDEBUG
  if (dy_opts->print.crash >= 4)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
	    "\n\t%d eligible architecturals: %d free, %d 1-bound, %d 2-bound",
	    eligible,freecnt,onebndcnt,twobndcnt) ; }
    # endif

  if (aj != NULL) pkvec_free(aj) ;
  if (retval == FALSE)
  { if (archvars != NULL) FREE(archvars) ;
    return (FALSE) ; }
/*
  And sort them. If there aren't any candidates, and we supplied the vector
  for archvars, release it.
*/
  if (eligible > 1)
  { qsort(&archvars[0],eligible,sizeof(ibrank_struct),ib_archvcomp) ; }
  else
  if (eligible == 0)
  { if (*p_archvars == NULL)
    { if (archvars != NULL) FREE(archvars) ;
      archvars = NULL ; } }
/*
  Return the sorted list.
*/
  *p_cnt = eligible ;
  *p_archvars = archvars ;

  return (TRUE) ; }




static int ib_archvselect (int cnt, ibrank_struct *vars)
/*
  This routine tries to populate the basis using the ranked architectural
  variables listed in vars. Equalities are examined first, because if we can't
  cover them with architectural variables, we'll have to use artificials.

  The variables in vars are already ranked according to how desireable it is
  to place them in the basis, taking into account bounds and number of
  nonzeros in a column. Here we'll try to allocate the variables to basis
  positions based on the pivot element a<ii>, and a notion that in an ideal
  world we'd like to construct a lower diagonal basis matrix.
  
  Why lower diagonal? Because we can turn it into a diagonal matrix, which is
  trivially easy to invert! To require that the columns we select conform
  strictly to the lower diagonal view, when we pick a pivot a<ii>, we would
  need to require that a<k,i> = 0, k < i.

  In practice, this is pretty restrictive. So we settle for requiring that
  elimination using multiples of rows a<k>, k < i, have minimal impact on our
  chosen pivot a<ii>.  We do this by requiring that coefficients a<k,i>, k <
  i, be small compared to the pivots a<kk>: specifically, a<k,i> < .1
  a<kk>. (If you're having trouble visualising this, you might want to go off
  and do a quick 3x3 example about now.)

  As we're selecting, we also try to pay a bit of attention to numerical
  stability.  If the matrix is scaled, the maximum coefficient in each row
  and column is 1.0, and we simply require |a<ij>| > .9. If the matrix is
  unscaled, well, we settle for |a<ij>/a<max,j>| > .9.  (Arguably we should
  also look at a<p,max>, but if you're really having stability trouble, you
  should consider allowing dylp to scale the matrix.)

  Note that we don't have to make marginal choices here --- we'll fill any
  vacant positions with the logical (slack or artificial) for the constraint.

  Parameters:
    cnt:	the number of variables in vars
    vars:	ranked vector of architectural variables

  Returns: the number of basis positions covered, or -1 if an error occurs.
*/

{ int i,j,m,ndx,pkndx,covered ;
  double ratio ;
  double *estpiv ;
  ibrank_struct *var ;

  bool scaled,select ;
  pkvec_struct *aj ;
  pkcoeff_struct *aij ;

  dyret_enum retval ;

  const char *rtnnme = "ib_archvselect" ;

# ifndef DYLP_NDEBUG
  int bndcnt ;
  double maxratio ;
# endif

  retval = 0 ;
  ratio = 0 ;
/*
  Allocate bookkeeping arrays.
*/
  m = dy_sys->concnt ;
  scaled = dy_isscaled() ;

  aj = NULL ;
  estpiv = (double *) MALLOC((m+1)*sizeof(double)) ;
  for (i = 1 ; i <= m ; i++) estpiv[i] = dy_tols->inf ;
  covered = 0 ;
/*
  Open a loop to scan vars and attempt to use each variable.
*/
  for (ndx = 0 ; ndx < cnt ; ndx++)
  { var = &vars[ndx] ;
    j = var->ndx ;
    if (consys_getcol_pk(dy_sys,j,&aj) == FALSE)
    { errmsg(122,rtnnme,dy_sys->nme,"column",
	     consys_nme(dy_sys,'v',j,TRUE,NULL),j) ;
      retval = -1 ;
      break ; }
#   ifndef DYLP_NDEBUG
    bndcnt = 0 ;
    if (dy_sys->vlb[j] > -dy_tols->inf) bndcnt++ ;
    if (dy_sys->vub[j] < dy_tols->inf) bndcnt++ ;
    maxratio = -1.0 ;
#   endif
/*
  Scan the column: We need a coefficient in an uncovered row (dy_basis[i] ==
  0) which satisfies the stability criteria (ratio > .9). The coefficients we
  scan before finding a suitable pivot can't do too much damage to our goal
  of lower diagonal (|aij->val| < .1*estpiv[i]).

  NOTE: Scan order is crucially important here. We need to consider
	equalities before inequalities in order to minimize eventual use of
	artificials to cover equalities. Because of the way dy_sys is loaded
	(equalities first) and the way consys augments columns (new
	coefficients go at the front) we need to scan from back to front in
	order to consider equalities first.
*/
    select = FALSE ;
    for (pkndx = aj->cnt-1 ; pkndx >= 0 ; pkndx--)
    { aij = &aj->coeffs[pkndx] ;
      i = aij->ndx ;
      ratio = fabs(aij->val) ;
      if (!scaled) ratio /= var->ajmax ;
      if (dy_basis[i] == 0)
      { if (ratio > .9)
	{ select = TRUE ;
	  break ; }
#     ifndef DYLP_NDEBUG
	else
	{ if (maxratio < ratio) maxratio = ratio ; }
#     endif
      }
      else
      { if (ratio > .1*estpiv[i]) break ; } }
/*
  Did we select this variable? If so, install it. If we have a full basis,
  break out of the search loop.
*/
    if (select == TRUE)
    { dy_basis[i] = j ;
      dy_var2basis[j] = i ;
      estpiv[i] = ratio ;
      covered++ ;
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.crash >= 4)
      { dyio_outfmt(dy_logchn,dy_gtxecho,"\n\t  adding %s (%d) (%d bounds)",
		    consys_nme(dy_sys,'v',j,FALSE,NULL),j,bndcnt) ;
	dyio_outfmt(dy_logchn,dy_gtxecho," to cover %s %s (%d),",
		    consys_prtcontyp(dy_sys->ctyp[i]),
		    consys_nme(dy_sys,'c',i,FALSE,NULL),i) ;
	dyio_outfmt(dy_logchn,dy_gtxecho," |a<%d,%d>/max(a<*,%d>)| = %g.",
		    i,j,j,ratio) ; }
#     endif
      if (covered == m) break ; }
#   ifndef DYLP_NDEBUG
    else
    if (dy_opts->print.crash >= 4)
    { dyio_outfmt(dy_logchn,dy_gtxecho,"\n\t  rejected %s (%d) (%d bounds)",
		  consys_nme(dy_sys,'v',j,FALSE,NULL),j,bndcnt) ;
      if (pkndx >= 0)
      { dyio_outfmt(dy_logchn,dy_gtxecho,"; lower diag violation at .1;") ;
	dyio_outfmt(dy_logchn,dy_gtxecho,
		    " a<%d,%d> = %g, estpiv<%d> = %g, ratio %g.",
		    i,j,ratio,i,estpiv[i],ratio/estpiv[i]) ; }
      else
      { if (maxratio > 0)
	{ dyio_outfmt(dy_logchn,dy_gtxecho," at %g < .9; no suitable pivots.",
		      maxratio) ; }
	else
	{ dyio_outfmt(dy_logchn,dy_gtxecho,
		      " no non-zeroes in uncovered rows.") ; } } }
#   endif
  }
  if (aj != NULL) pkvec_free(aj) ;
  if (estpiv != NULL) FREE(estpiv) ;
  if (retval < 0) return (retval) ;
/*
  We're done. Print a summary and return.
*/

# ifndef DYLP_NDEBUG
  { if (dy_opts->print.crash >= 2)
    { dyio_outfmt(dy_logchn,dy_gtxecho,
		  "\n    added %d architectural variables.",covered) ; } }
# endif

  return (covered) ; }



static int ib_slackselect (void)

/*
  This routine will populate the basis slots corresponding to inequalities
  using the logical (slack) variable for the inequality.

  Parameters: none

  Returns: the number of basis slots filled.
*/

{ int i,m,covered ;
  contyp_enum *ctyp ;

/*
  Walk the constraints, installing the logical for each inequality.
*/
  m = dy_sys->concnt ;
  ctyp = dy_sys->ctyp ;
  covered = 0 ;

  for (i = 1 ; i <= m ; i++)
  { if (ctyp[i] != contypEQ && dy_basis[i] == 0)
    { dy_basis[i] = i ;
      dy_var2basis[i] = i ;
      covered++ ;
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.crash >= 4)
      { dyio_outfmt(dy_logchn,dy_gtxecho,"\n\t  adding %s (%d)",
		    consys_nme(dy_sys,'v',i,FALSE,NULL),i) ;
	dyio_outfmt(dy_logchn,dy_gtxecho," to cover %s (%d).",
		    consys_nme(dy_sys,'c',i,FALSE,NULL),i) ; }
#     endif
    } }
# ifndef DYLP_NDEBUG
  if (dy_opts->print.crash >= 2)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
		"\n    added %d slack/surplus variables.",covered) ; }
# endif

  return (covered) ; }


static int ib_artifselect (void)

/*
  This routine will populate the basis slots corresponding to equalities using
  the logical (artificial) variable for the equality.

  Parameters: none

  Returns: the number of basis slots filled
*/

{ int i,m,covered ;
  contyp_enum *ctyp ;

/*
  Walk the constraints, installing the logical for each inequality.
*/
  m = dy_sys->concnt ;
  ctyp = dy_sys->ctyp ;
  covered = 0 ;

  for (i = 1 ; i <= dy_sys->concnt ; i++)
  { if (ctyp[i] == contypEQ && dy_basis[i] == 0)
    { dy_basis[i] = i ;
      dy_var2basis[i] = i ;
      covered++ ;
#     ifndef DYLP_NDEBUG
      if (dy_opts->print.crash >= 4)
      { dyio_outfmt(dy_logchn,dy_gtxecho,"\n\t  adding %s (%d)",
		    consys_nme(dy_sys,'v',i,FALSE,NULL),i) ;
	dyio_outfmt(dy_logchn,dy_gtxecho," to cover %s (%d).",
		    consys_nme(dy_sys,'c',i,FALSE,NULL),i) ; }
#     endif
    } }
# ifndef DYLP_NDEBUG
  if (dy_opts->print.crash >= 2)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
	        "\n    added %d artificial variables.",covered) ; }
# endif

  return (covered) ; }



static bool ib_populatebasis (void)

/*
  This routine is responsible for populating the basis (i.e., selecting basic
  variables for each basis position) using a mixture of architectural and
  logical variables. The type of basis is determined by the coldbasis option,
  as follows:
    ibLOGICAL:	(logical) uses only logical variables (slacks and artificials).
    ibSLACK:	(slack) uses slacks for inequalities, then covers as many
		equalities as possible with architecturals before falling back
		on artificials.
    ibARCH:	(architectural) uses architecturals first, trying to cover
		equalities before inequalities, and falling back on logicals
		for uncovered positions.

  Parameters: none

  Returns: TRUE if a basis is built, FALSE if something goes wrong.
*/

{ int m,basiscnt,archvcnt,slkcnt,artifcnt,iretval ;
  bool bretval ;

  ibrank_struct *archvars ;

  const char *rtnnme = "ib_populatebasis" ;

# ifndef DYLP_NDEBUG
  int i,j ;
# endif

# ifdef DYLP_PARANOIA
  if (!(dy_opts->coldbasis >= ibLOGICAL && dy_opts->coldbasis <= ibARCH))
  { errmsg(5,rtnnme,"initial basis type",dy_opts->coldbasis) ;
    return (FALSE) ; }
# endif

  m = dy_sys->concnt ;

# ifndef DYLP_NDEBUG
  if (dy_opts->print.crash >= 1)
  { dyio_outfmt(dy_logchn,dy_gtxecho,"\n  constructing ") ;
    switch (dy_opts->coldbasis)
    { case ibLOGICAL:
      { dyio_outfmt(dy_logchn,dy_gtxecho,"logical") ;
	break ; }
      case ibSLACK:
      { dyio_outfmt(dy_logchn,dy_gtxecho,"slack") ;
	break ; }
      case ibARCH:
      { dyio_outfmt(dy_logchn,dy_gtxecho,"architectural") ;
	break ; }
      default:
      { errmsg(1,rtnnme,__LINE__) ;
	return (FALSE) ; } }
    dyio_outfmt(dy_logchn,dy_gtxecho," basis for system %s, %d constraints.",
	        dy_sys->nme,m) ; }
# endif


  basiscnt = 0 ;

/*
  Logical and slack basis types will use slacks to cover off inequalities.
*/
  if (dy_opts->coldbasis == ibLOGICAL || dy_opts->coldbasis == ibSLACK)
  { slkcnt = ib_slackselect() ;
    basiscnt += slkcnt ; }
/*
  Slack and architectural basis types will use architectural variables. For a
  slack basis, we're just trying to cover equalities before falling back on
  artificial variables. For an architectural basis, we're covering equalities
  and inequalities, scanning the columns in a way that will see equalities
  before inequalities (hence preferentially use architecturals to cover
  equalities).

  Start by sorting them into their respective classes by number of bounds,
  then ordering within classes by the rating described in the header.

  Exclude fixed variables from the two-bound group. A simple equality test is
  all we need, as dy_coldstart has already groomed the bounds.
*/
  if (basiscnt < m &&
      (dy_opts->coldbasis == ibSLACK || dy_opts->coldbasis == ibARCH))
  { archvars = NULL ;
    bretval = ib_archvrank(&archvcnt,&archvars) ;
    if (bretval == FALSE)
    { errmsg(305,rtnnme,dy_sys->nme,
	     dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,
	     "rank","architectural") ;
      if (archvars != NULL) FREE(archvars) ;
      return (FALSE) ; }
    iretval = ib_archvselect(archvcnt,archvars) ;
    if (archvars != NULL) FREE(archvars) ;
    if (iretval < 0)
    { errmsg(305,rtnnme,dy_sys->nme,
	     dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,
	     "select","architectural") ;
      return (FALSE) ; }
    basiscnt += iretval ; }
/*
  Use logicals for any uncovered inequalities in an architectural basis.
*/
  if (basiscnt < m && dy_opts->coldbasis == ibARCH)
  { slkcnt = ib_slackselect() ;
    basiscnt += slkcnt ; }
/*
  If we have uncovered equalities (and that's all that can remain, at this
  point), give in and use artificial variables.
*/
  if (basiscnt < m)
  { artifcnt = ib_artifselect() ;
    basiscnt += artifcnt ; }

# ifndef DYLP_NDEBUG
  if (dy_opts->print.crash >= 4)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
		"\n\t    Pos'n Variable           Constraint") ;
    for (i = 1 ; i <= basiscnt ; i++)
    { j = dy_basis[i] ;
      dyio_outfmt(dy_logchn,dy_gtxecho,"\n\t     %3d  (%3d) %-15s",i,j,
		  consys_nme(dy_sys,'v',j,FALSE,NULL)) ;
      dyio_outfmt(dy_logchn,dy_gtxecho,"%-15s",
		  consys_nme(dy_sys,'c',i,FALSE,NULL)) ; } }
# endif

/*
  And if we still have uncovered rows, there's something seriously wrong.
*/
  if (basiscnt < m)
  { errmsg(301,rtnnme,basiscnt,m,dy_sys->nme) ;
    return (FALSE) ; }

  return (TRUE) ; }




dyret_enum dy_crash (void)

/*
  This routine coordinates the selection of an initial basis for an lp
  problem. The task of selecting a set of basic variables is hidden away in
  ib_populatebasis and subsidiary routines (see the comments at the head of
  this group of routines). What's visible here is the work associated with
  factoring the basis, establishing status, and calculating primal and dual
  variables and reduced costs.



  Parameters:  dy_lp:      lp problem dy_opts:    lp algorithm options
    dy_tols:    lp algorithm control and tolerances

  Returns: dyrOK if a basis is constructed; originates dyrFATAL on error, and
	   can relay dyrNUMERIC, dyrBSPACE, and dyrSINGULAR originated by
	   dy_factor.
*/

{ int vndx ;
  double *vub,*vlb,*obj ;
  flags calcflgs ;
  dyret_enum retval ;

  const char *rtnnme = "dy_crash" ;

  extern void dy_setfinalstatus(void) ;		/* dy_hotstart.c */

# ifndef DYLP_NDEBUG
  int cndx ;
# endif

# ifdef DYLP_PARANOIA
  if (dy_lp == NULL)
  { errmsg(2,rtnnme,"dy_lp") ;
    return (dyrFATAL) ; }
  if (dy_sys == NULL)
  { errmsg(2,rtnnme,"dy_sys") ;
    return (dyrFATAL) ; }
  if (dy_sys == NULL)
  { errmsg(2,rtnnme,consys_assocnme(dy_sys,CONSYS_MTX)) ;
    return (dyrFATAL) ; }
  if (dy_sys->vlb == NULL)
  { errmsg(2,rtnnme,consys_assocnme(dy_sys,CONSYS_VLB)) ;
    return (dyrFATAL) ; }
  if (dy_sys->obj == NULL)
  { errmsg(2,rtnnme,consys_assocnme(dy_sys,CONSYS_OBJ)) ;
    return (dyrFATAL) ; }
  if (dy_sys->ctyp == NULL)
  { errmsg(2,rtnnme,consys_assocnme(dy_sys,CONSYS_CTYP)) ;
    return (dyrFATAL) ; }
  if (flgoff(dy_sys->opts,CONSYS_LVARS))
  { errmsg(311,rtnnme,dy_sys->nme) ;
    return (dyrFATAL) ; }
  if (dy_opts == NULL)
  { errmsg(2,rtnnme,"dy_opts") ;
    return (dyrFATAL) ; }
  if (dy_tols == NULL)
  { errmsg(2,rtnnme,"dy_tols") ;
    return (dyrFATAL) ; }
# endif
/*
  Unpack a few arrays we'll use frequently, and attach the basis and inverse
  basis vectors to the constraint system. consys_attach will initialise them
  to 0.
*/
  obj = dy_sys->obj ;
  vlb = dy_sys->vlb ;
  vub = dy_sys->vub ;
  if (consys_attach(dy_sys,CONSYS_COL,
		    sizeof(int),(void **) &dy_basis) == FALSE)
  { errmsg(100,rtnnme,dy_sys->nme,"basis vector") ;
    return (dyrFATAL) ; }
  if (consys_attach(dy_sys,CONSYS_ROW,
		    sizeof(int),(void **) &dy_var2basis) == FALSE)
  { errmsg(100,rtnnme,dy_sys->nme,"inverse basis vector") ;
    return (dyrFATAL) ; }
/*
  Populate the basis.
*/
  if (ib_populatebasis() == FALSE)
  { errmsg(302,rtnnme,dy_sys->nme,
	   dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters,"populate") ;
    return (dyrFATAL) ; }
/*
  Factor the basis. We don't want any of the primal or dual variables
  calculated just yet. If this fails we're in deep trouble.
*/
  if (dy_sys->concnt > 0)
  {
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.crash >= 2)
      dyio_outfmt(dy_logchn,dy_gtxecho,"\n\tfactoring ...") ;
#   endif
    calcflgs = 0 ;
    retval = dy_factor(&calcflgs) ;
    switch (retval)
    { case dyrOK:
      case dyrPATCHED:
      { break ; }
      default:
      { errmsg(309,rtnnme,dy_sys->nme) ;
	return (retval) ; } } }
/*
  Attach and clear the vectors which will hold the status and values of
  variables, and the reduced cost.
*/
  if (consys_attach(dy_sys,CONSYS_ROW,
		    sizeof(flags),(void **) &dy_status) == FALSE)
  { errmsg(100,rtnnme,dy_sys->nme,"status vector") ;
    return (dyrFATAL) ; }
  if (consys_attach(dy_sys,CONSYS_COL,
		    sizeof(double),(void **) &dy_xbasic) == FALSE)
  { errmsg(100,rtnnme,dy_sys->nme,"basic variable vector") ;
    return (dyrFATAL) ; }
  if (consys_attach(dy_sys,CONSYS_ROW,
		    sizeof(double),(void **) &dy_x) == FALSE)
  { errmsg(100,rtnnme,dy_sys->nme,"primal variable vector") ;
    return (dyrFATAL) ; }
  if (consys_attach(dy_sys,CONSYS_COL,
		    sizeof(double),(void **) &dy_y) == FALSE)
  { errmsg(100,rtnnme,dy_sys->nme,"dual variable vector") ;
    return (dyrFATAL) ; }
  if (consys_attach(dy_sys,CONSYS_ROW,
		    sizeof(double),(void **) &dy_cbar) == FALSE)
  { errmsg(100,rtnnme,dy_sys->nme,"reduced cost vector") ;
    return (dyrFATAL) ; }
/*
  Might as well work toward a dual feasible start. Calculate the duals, then
  the reduced costs, so we can make intelligent decisions about the status of
  the nonbasic variables.
*/
  dy_calcduals() ;
  if (dy_calccbar() == FALSE)
  { errmsg(384,rtnnme,dy_sys->nme,
	   dy_prtlpphase(dy_lp->phase,TRUE),dy_lp->tot.iters) ;
    return (dyrFATAL) ; }
/*
  For nonbasic variables, we have to decide the proper status, based on
  number of bounds and the sign of the reduced cost (we're minimising,
  remember). Once we know the status, we can set a value in dy_x.  Nonbasic
  free variables are arbitrarily awarded a value of 0.  For basic variables,
  if we can't say they're fixed or free, assume strictly between bounds until
  we calculate their initial values. We won't bother to set dy_x and
  dy_xbasic until we've done the calculation.
*/

# ifndef DYLP_NDEBUG
  if (dy_opts->print.crash >= 2)
    dyio_outfmt(dy_logchn,dy_gtxecho,
	        "\n\testablishing initial status and reference frame ...") ;
# endif
  for (vndx = 1 ; vndx <= dy_sys->varcnt ; vndx++)
  { if (dy_var2basis[vndx] != 0)
    { if (vlb[vndx] == vub[vndx])
      { dy_status[vndx] = vstatBFX ; }
      else
      if (vlb[vndx] <= -dy_tols->inf && vub[vndx] >= dy_tols->inf)
      { dy_status[vndx] = vstatBFR ; }
      else
      { dy_status[vndx] = vstatB ; } }
    else
    { if (vlb[vndx] > -dy_tols->inf && vub[vndx] < dy_tols->inf)
      { if (vub[vndx] == vlb[vndx])
	{ dy_status[vndx] = vstatNBFX ;
	  dy_x[vndx] = vub[vndx] ; }
	else
	if (dy_cbar[vndx] >= 0)
	{ dy_status[vndx] = vstatNBLB ;
	  dy_x[vndx] = vlb[vndx] ; }
	else
	{ dy_status[vndx] = vstatNBUB ;
	  dy_x[vndx] = vub[vndx] ; } }
      else
      if (vlb[vndx] > -dy_tols->inf)
      { dy_status[vndx] = vstatNBLB ;
	dy_x[vndx] = vlb[vndx] ; }
      else
      if (vub[vndx] < dy_tols->inf)
      { dy_status[vndx] = vstatNBUB ;
	dy_x[vndx] = vub[vndx] ; }
      else
      { dy_status[vndx] = vstatNBFR ;
	dy_x[vndx] = 0 ; } }
#   ifndef DYLP_NDEBUG
    if (dy_opts->print.crash >= 4)
    { dyio_outfmt(dy_logchn,dy_gtxecho,"\n\t  %s (%d) %s",
		  consys_nme(dy_sys,'v',vndx,FALSE,NULL),vndx,
		  dy_prtvstat(dy_status[vndx])) ;
      if (flgon(dy_status[vndx],vstatNONBASIC|vstatNBFR))
	dyio_outfmt(dy_logchn,dy_gtxecho," with value %g.",dy_x[vndx]) ;
      else
	dyio_outchr(dy_logchn,dy_gtxecho,'.') ; }
#   endif
  }
/*
  Ok, status is set. Now it's time to calculate initial values for the basic
  variables and objective.  Once we have values for the basic variables, see
  how bad it looks, in terms of primal infeasibility, and adjust the status
  for any that are pinned against a bound or out of bounds.
*/
# ifndef DYLP_NDEBUG
  if (dy_opts->print.crash >= 2)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
		"\n\tcalculating basic variable values ...") ; }
# endif
  if (dy_calcprimals() == FALSE)
  { errmsg(316,rtnnme,dy_sys->nme) ;
    return (dyrFATAL) ; }
  dy_lp->z = dy_calcobj() ;
  dy_setfinalstatus() ;
/*
  Make the check for primal and/or dual feasibility, and set the initial
  simplex phase accordingly.
*/
  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 (retval) ; }
  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 ; }

# ifndef DYLP_NDEBUG
/*
  Some debug printing, to dump the initial basis & variables.
*/
  if (dy_opts->print.crash >= 3)
  { dyio_outfmt(dy_logchn,dy_gtxecho,
	        "\n  Pos'n Constraint\tDual\tPrimal\t\t   Status\tValue") ;
    for (cndx = 1 ; cndx <= dy_sys->concnt ; cndx++)
    { vndx = dy_basis[cndx] ;
      dyio_outfmt(dy_logchn,dy_gtxecho,"\n   %4d  %-13s%7g",cndx,
		  consys_nme(dy_sys,'c',cndx,FALSE,NULL),dy_y[cndx]) ;
      dyio_outfmt(dy_logchn,dy_gtxecho,"  (%4d) %-13s %-7s %7g",vndx,
		  consys_nme(dy_sys,'v',vndx,FALSE,NULL),
		  dy_prtvstat(dy_status[vndx]),dy_x[vndx]) ; } }

  if (dy_opts->print.crash >= 2)
  { dyio_outfmt(dy_logchn,dy_gtxecho,"\n\tinitial objective %g",dy_lp->z) ;
    if (dy_lp->infeascnt != 0)
    { dyio_outfmt(dy_logchn,dy_gtxecho,", %d infeasible vars, infeas. = %g",
		  dy_lp->infeascnt,dy_lp->infeas) ; }
    dyio_outfmt(dy_logchn,dy_gtxecho,", target simplex %s.",
	        dy_prtlpphase(dy_lp->simplex.next,FALSE)) ; }
# endif

  return (dyrOK) ; }
