if(quot < (*theta)) {
	  (*theta) = quot;
	  (*row_nr) = i;
	}
      }
    }

  if((*theta) < 0) {
    fprintf(stderr, "Warning: Numerical instability, qout = %g\n",
	    (double)(*theta));
    fprintf(stderr, "pcol[%d] = %18g, rhs[%d] = %18g , upbo = %g\n",
	    (*row_nr), (double)f, (*row_nr), (double)lp->rhs[(*row_nr)],
	    (double)lp->upbo[lp->bas[(*row_nr)]]);
  }
  if((*row_nr) == 0) {
    if(lp->upbo[colnr] == lp->infinite) {
      Doiter   = FALSE;
      DoInvert = FALSE;
      Status   = UNBOUNDED;
    }
    else {
      i = 1;
      while(pcol[i] >= 0 && i <= lp->rows)
	i++;
      if(i > lp->rows) { /* empty column with upperbound! */
	lp->lower[colnr] = FALSE;
	lp->rhs[0] += lp->upbo[colnr]*pcol[0];
	Doiter = FALSE;
	DoInvert = FALSE;
      }
      else if(pcol[i]<0) {
	(*row_nr) = i;
      }
    }
  }
  if((*row_nr) > 0)
    Doiter = TRUE;
  if(lp->trace)
    fprintf(stderr, "row_prim:%d, pivot element:%18g\n", (*row_nr),
	    (double)pcol[(*row_nr)]);

  return((*row_nr) > 0);
} /* rowprim */

static short rowdual(lprec *lp, int *row_nr)
{
  int   i;
  REAL  f, g, minrhs;
  short artifs;

  (*row_nr) = 0;
  minrhs = -lp->epsb;
  i = 0;
  artifs = FALSE;
  while(i < lp->rows && !artifs) {
    i++;
    f = lp->upbo[lp->bas[i]];
    if(f == 0 && (lp->rhs[i] != 0)) {
      artifs = TRUE;
      (*row_nr) = i;
    }
    else {
      if(lp->rhs[i] < f - lp->rhs[i])
	g = lp->rhs[i];
      else
	g = f - lp->rhs[i];
      if(g < minrhs) {
	minrhs = g;
	(*row_nr) = i;
      }
    }
  }

  if(lp->trace) {  
    if((*row_nr) > 0) { 
      fprintf(stderr,
	      "row_dual:%d, rhs of selected row:           %18g\n",
	      (*row_nr), (double)lp->rhs[(*row_nr)]);
      if(lp->upbo[lp->bas[(*row_nr)]] < lp->infinite)
	fprintf(stderr,
		"\t\tupper bound of basis variable:    %18g\n",
		(double)lp->upbo[lp->bas[(*row_nr)]]);
    }
    else
      fprintf(stderr, "row_dual: no infeasibilities found\n");
  }
    
  return((*row_nr) > 0);
} /* rowdual */

static short coldual(lprec *lp,
		     int row_nr,
		     int *colnr,
		     short minit,
		     REAL *prow,
		     REAL *drow)
{
  int  i, j, k, r, varnr, *rowp, row;
  REAL theta, quot, pivot, d, f, g, *valuep, value;
  
  Doiter = FALSE;
  if(!minit) {
    for(i = 0; i <= lp->rows; i++) {
      prow[i] = 0;
      drow[i] = 0;
    }

    drow[0] = 1;
    prow[row_nr] = 1;

    for(i = lp->eta_size; i >= 1; i--) {
      d = 0;
      f = 0;
      k = lp->eta_col_end[i] - 1;
      r = lp->eta_row_nr[k];
      j = lp->eta_col_end[i - 1];
      
      /* this is one of the loops where the program consumes a lot of CPU
	 time */
      /* let's help the compiler by doing some pointer arithmetic instead
	 of array indexing */
      for(rowp = lp->eta_row_nr + j, valuep = lp->eta_value + j;
	  j <= k;
	  j++, rowp++, valuep++) {
	f += prow[*rowp] * *valuep;
	d += drow[*rowp] * *valuep;
      }

      my_round(f, lp->epsel);
      prow[r] = f;
      my_round(d, lp->epsd);
      drow[r] = d;
    }

    for(i = 1; i <= lp->columns; i++) {
      varnr = lp->rows + i;
      if(!lp->basis[varnr]) {
	matrec *matentry;

	d = - Extrad * drow[0];
	f = 0;
	k = lp->col_end[i];
	j = lp->col_end[i - 1];

	/* this is one of the loops where the program consumes a lot
	   of cpu time */
	/* let's help the compiler with pointer arithmetic instead
	   of array indexing */
	for(matentry = lp->mat + j;
	    j < k;
	    j++, matentry++) {
	  row = (*matentry).row_nr;
	  value = (*matentry).value;
	  d += drow[row] * value;
	  f += prow[row] * value;
	}

	my_round(f, lp->epsel);
	prow[varnr] = f;
	my_round(d, lp->epsd);
	drow[varnr] = d;
      }
    }
  }

  if(lp->rhs[row_nr] > lp->upbo[lp->bas[row_nr]])
    g = -1;
  else
    g = 1;

  pivot = 0;
  (*colnr) = 0;
  theta = lp->infinite;

  for(i = 1; i <= lp->sum; i++) {
    if(lp->lower[i])
      d = prow[i] * g;
    else
      d = -prow[i] * g;
    
    if((d < 0) && (!lp->basis[i]) && (lp->upbo[i] > 0)) {
      if(lp->lower[i])
	quot = -drow[i] / (REAL) d;
      else
	quot = drow[i] / (REAL) d;
      if(quot < theta) {
	theta = quot;
	pivot = d;
	(*colnr) = i;
      }
      else if((quot == theta) && (my_abs(d) > my_abs(pivot))) {
	pivot = d;
	(*colnr) = i;
      }
    }
  }

  if(lp->trace)
    fprintf(stderr, "col_dual:%d, pivot element:  %18g\n", (*colnr),
	    (double)prow[(*colnr)]);

  if((*colnr) > 0)
    Doiter = TRUE;

  return((*colnr) > 0);
} /* coldual */

static void iteration(lprec *lp,
		      int row_nr,
		      int varin,
		      REAL *theta,
		      REAL up,
		      short *minit,
		      short *low,
		      short primal)
{
  int  i, k, varout;
  REAL f;
  REAL pivot;
  
  lp->iter++;
  
  if(((*minit) = (*theta) > (up + lp->epsb))) {
    (*theta) = up;
    (*low) = !(*low);
  }

  k = lp->eta_col_end[lp->eta_size + 1];
  pivot = lp->eta_value[k - 1];

  for(i = lp->eta_col_end[lp->eta_size]; i < k; i++) {
    f = lp->rhs[lp->eta_row_nr[i]] - (*theta) * lp->eta_value[i];
    my_round(f, lp->epsb);
    lp->rhs[lp->eta_row_nr[i]] = f;
  }

  if(!(*minit)) {
    lp->rhs[row_nr] = (*theta);
    varout = lp->bas[row_nr];
    lp->bas[row_nr] = varin;
    lp->basis[varout] = FALSE;
    lp->basis[varin] = TRUE;

    if(primal && pivot < 0)
      lp->lower[varout] = FALSE;

    if(!(*low) && up < lp->infinite) {
      (*low) = TRUE;
      lp->rhs[row_nr] = up - lp->rhs[row_nr];
      for(i = lp->eta_col_end[lp->eta_size]; i < k; i++)
	lp->eta_value[i] = -lp->eta_value[i];
    }

    addetacol(lp);
    lp->num_inv++;
  }

  if(lp->trace) {
    fprintf(stderr, "Theta = %g ", (double)(*theta));
    if((*minit)) {
      if(!lp->lower[varin])
	fprintf(stderr,
		"Iteration: %d, variable %d changed from 0 to its upper bound of %g\n",
		lp->iter, varin, (double)lp->upbo[varin]);
      else
	fprintf(stderr,
		"Iteration: %d, variable %d changed its upper bound of %g to 0\n",
		lp->iter, varin, (double)lp->upbo[varin]);
    }
    else
      fprintf(stderr,
	      "Iteration: %d, variable %d entered basis at: %g\n",
	      lp->iter, varin, (double)lp->rhs[row_nr]);
    if(!primal) {
      f = 0;
      for(i = 1; i <= lp->rows; i++)
	if(lp->rhs[i] < 0)
	  f -= lp->rhs[i];
	else
	  if(lp->rhs[i] > lp->upbo[lp->bas[i]])
	    f += lp->rhs[i] - lp->upbo[lp->bas[i]];
      fprintf(stderr, "feasibility gap of this basis: %g\n",
	      (double)f);
    }
    else
      fprintf(stderr,
	      "objective function value of this feasible basis: %g\n",
	      (double)lp->rhs[0]);
  }
} /* iteration */


static int solvelp(lprec *lp)
{
  int    i, j, varnr;
  REAL   f, theta;
  short  primal;
  REAL   *drow, *prow, *Pcol;
  short  minit;
  int    colnr, row_nr;
  short  *test; 

  if(lp->do_presolve)
    presolve(lp);

  CALLOC(drow, lp->sum + 1);
  CALLOC(prow, lp->sum + 1);
  CALLOC(Pcol, lp->rows + 1);
  CALLOC(test, lp->sum +1); 

  lp->iter = 0;
  minit = FALSE;
  Status = RUNNING;
  DoInvert = FALSE;
  Doiter = FALSE;

  for(i = 1, primal = TRUE; (i <= lp->rows) && primal; i++)
    primal = (lp->rhs[i] >= 0) && (lp->rhs[i] <= lp->upbo[lp->bas[i]]);

  if(lp->trace) {
    if(primal)
      fprintf(stderr, "Start at feasible basis\n");
    else
      fprintf(stderr, "Start at infeasible basis\n");
  }

  if(!primal) {
    drow[0] = 1;

    for(i = 1; i <= lp->rows; i++)
      drow[i] = 0;

	/* fix according to Joerg Herbers */
	btran(lp, drow);

    Extrad = 0;

    for(i = 1; i <= lp->columns; i++) {
      varnr = lp->rows + i;
      drow[varnr] = 0;

      for(j = lp->col_end[i - 1]; j < lp->col_end[i]; j++)
	if(drow[lp->mat[j].row_nr] != 0)
	  drow[varnr] += drow[lp->mat[j].row_nr] * lp->mat[j].value;

      if(drow[varnr] < Extrad)
	Extrad = drow[varnr];
    }
  }
  else
    Extrad = 0;

  if(lp->trace)
    fprintf(stderr, "Extrad = %g\n", (double)Extrad);

  minit = FALSE;

  while(Status == RUNNING) {
    Doiter = FALSE;
    DoInvert = FALSE;

    if(primal) {
      if(colprim(lp, &colnr, minit, drow)) {
	setpivcol(lp, lp->lower[colnr], colnr, Pcol);
	
	if(rowprim(lp, colnr, &row_nr, &theta, Pcol))
	  condensecol(lp, row_nr, Pcol);
      }
    }
    else /* not primal */ {
      if(!minit)
	rowdual(lp, &row_nr);

      if(row_nr > 0 ) {
	if(coldual(lp, row_nr, &colnr, minit, prow, drow)) {
	  setpivcol(lp, lp->lower[colnr], colnr, Pcol);

	  /* getting div by zero here. Catch it and try to recover */
	  if(Pcol[row_nr] == 0) {
	    fprintf(stderr,
		    "An attempt was made to divide by zero (Pcol[%d])\n",
		    row_nr);
	    fprintf(stderr,
		    "This indicates numerical instability\n");
	    Doiter = FALSE;
	    if(!JustInverted) {
	      fprintf(stderr,
		      "Trying to recover. Reinverting Eta\n");
	      DoInvert = TRUE;
	    }
	    else {
	      fprintf(stderr, "Can't reinvert, failure\n");
	      Status = FAILURE;
	    }
	  }
	  else {
	    condensecol(lp, row_nr, Pcol);
	    f = lp->rhs[row_nr] - lp->upbo[lp->bas[row_nr]];

	    if(f > 0) {
	      theta = f / (REAL) Pcol[row_nr];
	      if(theta <= lp->upbo[colnr])
		lp->lower[lp->bas[row_nr]] = !lp->lower[lp->bas[row_nr]];
	    }
	    else /* f <= 0 */
	      theta = lp->rhs[row_nr] / (REAL) Pcol[row_nr];
	  }
	}
	else
	  Status = INFEASIBLE;
      }
      else {
	primal   = TRUE;
	Doiter   = FALSE;
	Extrad   = 0;
	DoInvert = TRUE;
      }	  
    }

    if(Doiter)
      iteration(lp, row_nr, colnr, &theta, lp->upbo[colnr], &minit,
		&lp->lower[colnr], primal);
    
    if(lp->num_inv >= lp->max_num_inv)
      DoInvert = TRUE;

    if(DoInvert) {
      if(lp->print_at_invert)
	fprintf(stderr, "Inverting: Primal = %d\n", primal);
      invert(lp);
    }
  } 

  lp->total_iter += lp->iter;
 
  free(drow);
  free(prow);
  free(Pcol);
  free(test);

  return(Status);
} /* solvelp */


static short is_int(lprec *lp, int i)
{
  REAL   value, error;

  value = lp->solution[i];
  error = value - (REAL)floor((double)value);

  if(error < lp->epsilon)
    return(TRUE);

  if(error > (1 - lp->epsilon))
    return(TRUE);

  return(FALSE);
} /* is_int */


static void construct_solution(lprec *lp)
{
  int    i, j, basi;
  REAL   f;

  /* zero all results of rows */
  memset(lp->solution, '\0', (lp->rows + 1) * sizeof(REAL));

  lp->solution[0] = -lp->orig_rh[0];

  if(lp->scaling_used) {
    lp->solution[0] /= lp->scale[0];

    for(i = lp->rows + 1; i <= lp->sum; i++)
      lp->solution[i] = lp->lowbo[i] * lp->scale[i];

    for(i = 1; i <= lp->rows; i++) {
      basi = lp->bas[i];
      if(basi > lp->rows)
	lp->solution[basi] += lp->rhs[i] * lp->scale[basi];
    }
    for(i = lp->rows + 1; i <= lp->sum; i++)
      if(!lp->basis[i] && !lp->lower[i])
	lp->solution[i] += lp->upbo[i] * lp->scale[i];

    for(j = 1; j <= lp->columns; j++) {
      f = lp->solution[lp->rows + j];
      if(f != 0)
	for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
	  lp->solution[lp->mat[i].row_nr] += (f / lp->scale[lp->rows+j])
	    * (lp->mat[i].value / lp->scale[lp->mat[i].row_nr]);
    }
  
    for(i = 0; i <= lp->rows; i++) {
      if(my_abs(lp->solution[i]) < lp->epsb)
	lp->solution[i] = 0;
      else if(lp->ch_sign[i])
	lp->solution[i] = -lp->solution[i];
    }
  }
  else { /* no scaling */
    for(i = lp->rows + 1; i <= lp->sum; i++)
      lp->solution[i] = lp->lowbo[i];

    for(i = 1; i <= lp->rows; i++) {
      basi = lp->bas[i];
      if(basi > lp->rows)
	lp->solution[basi] += lp->rhs[i];
    }

    for(i = lp->rows + 1; i <= lp->sum; i++)
      if(!lp->basis[i] && !lp->lower[i])
	lp->solution[i] += lp->upbo[i];

    for(j = 1; j <= lp->columns; j++) {
      f = lp->solution[lp->rows + j];
      if(f != 0)
	for(i = lp->col_end[j - 1]; i < lp->col_end[j]; i++)
	  lp->solution[lp->mat[i].row_nr] += f * lp->mat[i].value;
    }
  
    for(i = 0; i <= lp->rows; i++) {
      if(my_abs(lp->solution[i]) < lp->epsb)
	lp->solution[i] = 0;
      else if(lp->ch_sign[i])
	lp->solution[i] = -lp->solution[i];
    }
  }
} /* construct_solution */

static void calculate_duals(lprec *lp)
{
  int i;

  /* initialize */
  lp->duals[0] = 1;
  for(i = 1; i <= lp->rows; i++)
    lp->duals[i] = 0;

  btran(lp, lp->duals);

  if(lp->scaling_used)
    for(i = 1; i <= lp->rows; i++)
      lp->duals[i] *= lp->scale[i] / lp->scale[0];

  /* the dual values are the reduced costs of the slacks */
  /* When the slack is at its upper bound, change the sign. */
  for(i = 1; i <= lp->rows; i++) {
    if(lp->basis[i])
      lp->duals[i] = 0;
    /* added a test if variable is different from 0 because sometime you get
       -0 and this is different from 0 on for example INTEL processors (ie 0
       != -0 on INTEL !) PN */
    else if((lp->ch_sign[0] == lp->ch_sign[i]) && lp->duals[i])
      lp->duals[i] = - lp->duals[i];
  }
} /* calculate_duals */

static void check_if_less(REAL x,
			  REAL y,
			  REAL value)
{
  if(x >= y) {
    fprintf(stderr,
	    "Error: new upper or lower bound is not more restrictive\n");
    fprintf(stderr, "bound 1: %g, bound 2: %g, value: %g\n",
	    (double)x, (double)y, (double)value);
    /* exit(EXIT_FAILURE); */
  }
}

static void check_solution(lprec *lp,
			   REAL *upbo,
			   REAL *lowbo)
{
  int i;

  /* check if all solution values are within the bounds, but allow some margin
     for numerical errors */

#define CHECK_EPS 1e-2

  if(lp->columns_scaled)
    for(i = lp->rows + 1; i <= lp->sum; i++) {
      if(lp->solution[i] < lowbo[i] * lp->scale[i] - CHECK_EPS) {
	fprintf(stderr,
		"Error: variable %d (%s) has a solution (%g) smaller than its lower bound (%g)\n",
		i - lp->rows, lp->col_name[i - lp->rows],
		(double)lp->solution[i], (double)lowbo[i] * lp->scale[i]);
	/* abort(); */
      }