/* $Id: lpsolve.c,v 1.2 1999/01/03 02:07:20 sybalsky Exp $ (C) Copyright Venue, All Rights Reserved
 */

/************************************************************************/
/*									*/
/*	(C) Copyright 1989-95 Venue. All Rights Reserved.		*/
/*	Manufactured in the United States of America.			*/
/*									*/
/************************************************************************/

#include "version.h"

#include <stdlib.h>
#include "lpkit.h"
#include "lpglob.h"

#include "lispemul.h"

#ifdef DOS
#include "devif.h"
#endif /* DOS */

extern int KBDEventFlg;
extern int *KEYBUFFERING68k;
#ifdef DOS
extern MouseInterface currentmouse;
#endif /* DOS */

/* Globals used by solver */
short JustInverted;
short Status;
short Doiter;
short DoInvert;
short Break_bb;

void set_globals(lprec *lp) {
  sstate *st;
  int i;

  if (Lp != NULL) Lp->active = FALSE;
  Lp = lp;
  Rows = lp->rows;
  Columns = lp->columns;
  Sum = Rows + Columns;
  Non_zeros = lp->non_zeros;
  Mat = lp->mat;
  Col_no = lp->col_no;
  Col_end = lp->col_end;
  Row_end = lp->row_end;
  Rh = lp->rh;
  Rhs = lp->rhs;
  Orig_rh = lp->orig_rh;
  Must_be_int = lp->must_be_int;
  Orig_upbo = lp->orig_upbo;
  Orig_lowbo = lp->orig_lowbo;
  Upbo = lp->upbo;
  Lowbo = lp->lowbo;
  Bas = lp->bas;
  Basis = lp->basis;
  Lower = lp->lower;
  Eta_alloc = lp->eta_alloc;
  Eta_size = lp->eta_size;
  Num_inv = lp->num_inv;
  Eta_value = lp->eta_value;
  Eta_row_nr = lp->eta_row_nr;
  Eta_col_end = lp->eta_col_end;
  Solution = lp->solution;
  Best_solution = lp->best_solution;
  Infinite = lp->infinite;
  Epsilon = lp->epsilon;
  Epsb = lp->epsb;
  Epsd = lp->epsd;
  Epsel = lp->epsel;

  /* ?? MB */
  TREJ = TREJ;
  TINV = TINV;

  Maximise = lp->maximise;
  Floor_first = lp->floor_first;
  lp->active = TRUE;
}

void ftran(int start, int end, REAL *pcol) {
  int i, j, k, r;
  REAL theta;

  for (i = start; i <= end; i++) {
    k = Eta_col_end[i] - 1;
    r = Eta_row_nr[k];
    theta = pcol[r];
    if (theta != 0)
      for (j = Eta_col_end[i - 1]; j < k; j++)
        pcol[Eta_row_nr[j]] += theta * Eta_value[j]; /* cpu expensive line */
    pcol[r] *= Eta_value[k];
  }
  /* round small values to zero */
  for (i = 0; i <= Rows; i++) my_round(pcol[i], Epsel);
} /* ftran */

void btran(REAL *row) {
  int i, j, k;
  REAL f;

  for (i = Eta_size; i >= 1; i--) {
    f = 0;
    k = Eta_col_end[i] - 1;
    for (j = Eta_col_end[i - 1]; j <= k; j++) f += row[Eta_row_nr[j]] * Eta_value[j];
    my_round(f, Epsel);
    row[Eta_row_nr[k]] = f;
  }
} /* btran */

short Isvalid(lprec *lp) {
  int i, j, *rownum, *colnum;
  int *num, row_nr;

  if (!lp->row_end_valid) {
    MALLOC(num, lp->rows + 1, int);
    MALLOC(rownum, lp->rows + 1, int);
    for (i = 0; i <= lp->rows; i++) {
      num[i] = 0;
      rownum[i] = 0;
    }
    for (i = 0; i < lp->non_zeros; i++) rownum[lp->mat[i].row_nr]++;
    lp->row_end[0] = 0;
    for (i = 1; i <= lp->rows; i++) lp->row_end[i] = lp->row_end[i - 1] + rownum[i];
    for (i = 1; i <= lp->columns; i++)
      for (j = lp->col_end[i - 1]; j < lp->col_end[i]; j++) {
        row_nr = lp->mat[j].row_nr;
        if (row_nr != 0) {
          num[row_nr]++;
          lp->col_no[lp->row_end[row_nr - 1] + num[row_nr]] = i;
        }
      }
    free(num);
    free(rownum);
    lp->row_end_valid = TRUE;
  }
  if (lp->valid) return (TRUE);
  CALLOC(rownum, lp->rows + 1, int);
  CALLOC(colnum, lp->columns + 1, int);
  for (i = 1; i <= lp->columns; i++)
    for (j = lp->col_end[i - 1]; j < lp->col_end[i]; j++) {
      colnum[i]++;
      rownum[lp->mat[j].row_nr]++;
    }
  for (i = 1; i <= lp->columns; i++)
    if (colnum[i] == 0)
      if (lp->names_used)
        fprintf(stderr, "Warning: Variable %s not used in any constraints\n", lp->col_name[i]);
      else
        fprintf(stderr, "Warning: Variable %d not used in any constaints\n", i);
  free(rownum);
  free(colnum);
  lp->valid = TRUE;
  return (TRUE);
}

static void resize_eta(void) {
  Eta_alloc *= 2;
  REALLOC(Eta_value, Eta_alloc, REAL);
  Lp->eta_value = Eta_value;
  REALLOC(Eta_row_nr, Eta_alloc, int);
  Lp->eta_row_nr = Eta_row_nr;
} /* resize_eta */

static void condensecol(int row_nr, REAL *pcol) {
  int i, elnr;

  elnr = Eta_col_end[Eta_size];

  if (elnr + Rows + 2 > Eta_alloc) /* maximum local growth of Eta */
    resize_eta();

  for (i = 0; i <= Rows; i++)
    if (i != row_nr && pcol[i] != 0) {
      Eta_row_nr[elnr] = i;
      Eta_value[elnr] = pcol[i];
      elnr++;
    }
  Eta_row_nr[elnr] = row_nr;
  Eta_value[elnr] = pcol[row_nr];
  elnr++;
  Eta_col_end[Eta_size + 1] = elnr;
} /* condensecol */

static void addetacol(void) {
  int i, j, k;
  REAL theta;

  j = Eta_col_end[Eta_size];
  Eta_size++;
  k = Eta_col_end[Eta_size];
  theta = 1 / (REAL)Eta_value[k - 1];
  Eta_value[k - 1] = theta;
  for (i = j; i < k - 1; i++) Eta_value[i] *= -theta;
  JustInverted = FALSE;
} /* addetacol */

static void setpivcol(short lower, int varin, REAL *pcol) {
  int i, colnr;
  REAL f;

  if (lower)
    f = 1;
  else
    f = -1;
  for (i = 0; i <= Rows; i++) pcol[i] = 0;
  if (varin > Rows) {
    colnr = varin - Rows;
    for (i = Col_end[colnr - 1]; i < Col_end[colnr]; i++) pcol[Mat[i].row_nr] = Mat[i].value * f;
    pcol[0] -= Extrad * f;
  } else if (lower)
    pcol[varin] = 1;
  else
    pcol[varin] = -1;
  ftran(1, Eta_size, pcol);
} /* setpivcol */

static void minoriteration(int colnr, int row_nr) {
  int i, j, k, wk, varin, varout, elnr;
  REAL piv, theta;

  varin = colnr + Rows;
  elnr = Eta_col_end[Eta_size];
  wk = elnr;
  Eta_size++;
  if (Extrad != 0) {
    Eta_row_nr[elnr] = 0;
    Eta_value[elnr] = -Extrad;
    elnr++;
  }
  for (j = Col_end[colnr - 1]; j < Col_end[colnr]; j++) {
    k = Mat[j].row_nr;
    if (k == 0 && Extrad != 0)
      Eta_value[Eta_col_end[Eta_size - 1]] += Mat[j].value;
    else if (k != row_nr) {
      Eta_row_nr[elnr] = k;
      Eta_value[elnr] = Mat[j].value;
      elnr++;
    } else
      piv = Mat[j].value;
  }
  Eta_row_nr[elnr] = row_nr;
  Eta_value[elnr] = 1 / (REAL)piv;
  elnr++;
  theta = Rhs[row_nr] / (REAL)piv;
  Rhs[row_nr] = theta;
  for (i = wk; i < elnr - 1; i++) Rhs[Eta_row_nr[i]] -= theta * Eta_value[i];
  varout = Bas[row_nr];
  Bas[row_nr] = varin;
  Basis[varout] = FALSE;
  Basis[varin] = TRUE;
  for (i = wk; i < elnr - 1; i++) Eta_value[i] /= -(REAL)piv;
  Eta_col_end[Eta_size] = elnr;
} /* minoriteration */

static void rhsmincol(REAL theta, int row_nr, int varin) {
  int i, j, k, varout;
  REAL f;

  if (row_nr > Rows + 1) {
    fprintf(stderr, "Error: rhsmincol called with row_nr: %d, rows: %d\n", row_nr, Rows);
    fprintf(stderr, "This indicates numerical instability\n");
    /* exit(1); */ ERROR(ERR_NUM);
  }
  j = Eta_col_end[Eta_size];
  k = Eta_col_end[Eta_size + 1];
  for (i = j; i < k; i++) {
    f = Rhs[Eta_row_nr[i]] - theta * Eta_value[i];
    my_round(f, Epsb);
    Rhs[Eta_row_nr[i]] = f;
  }
  Rhs[row_nr] = theta;
  varout = Bas[row_nr];
  Bas[row_nr] = varin;
  Basis[varout] = FALSE;
  Basis[varin] = TRUE;
} /* rhsmincol */

void invert(void) {
  int i, j, v, wk, numit, varnr, row_nr, colnr, varin;
  REAL f, theta;
  REAL *pcol;
  short *frow;
  short *fcol;
  int *rownum, *col, *row;
  int *colnum;

  if (Lp->print_at_invert)
    fprintf(stderr, "Start Invert iter %7d eta_size %4d rhs[0] %16.4f \n", Lp->iter, Eta_size,
            -Rhs[0]);

  CALLOC(rownum, Rows + 1, int);
  CALLOC(col, Rows + 1, int);
  CALLOC(row, Rows + 1, int);
  CALLOC(pcol, Rows + 1, REAL);
  CALLOC(frow, Rows + 1, short);
  CALLOC(fcol, Columns + 1, short);
  CALLOC(colnum, Columns + 1, int);

  for (i = 0; i <= Rows; i++) frow[i] = TRUE;

  for (i = 0; i < Columns; i++) fcol[i] = FALSE;

  for (i = 0; i < Rows; i++) rownum[i] = 0;

  for (i = 0; i <= Columns; i++) colnum[i] = 0;

  for (i = 0; i <= Rows; i++)
    if (Bas[i] > Rows)
      fcol[Bas[i] - Rows - 1] = TRUE;
    else
      frow[Bas[i]] = FALSE;

  for (i = 1; i <= Rows; i++)
    if (frow[i])
      for (j = Row_end[i - 1] + 1; j <= Row_end[i]; j++) {
        wk = Col_no[j];
        if (fcol[wk - 1]) {
          colnum[wk]++;
          rownum[i - 1]++;
        }
      }
  for (i = 1; i <= Rows; i++) Bas[i] = i;
  for (i = 1; i <= Rows; i++) Basis[i] = TRUE;
  for (i = 1; i <= Columns; i++) Basis[i + Rows] = FALSE;

  for (i = 0; i <= Rows; i++) Rhs[i] = Rh[i];

  for (i = 1; i <= Columns; i++) {
    varnr = Rows + i;
    if (!Lower[varnr]) {
      theta = Upbo[varnr];
      for (j = Col_end[i - 1]; j < Col_end[i]; j++) Rhs[Mat[j].row_nr] -= theta * Mat[j].value;
    }
  }
  for (i = 1; i <= Rows; i++)
    if (!Lower[i]) Rhs[i] -= Upbo[i];
  Eta_size = 0;
  v = 0;
  row_nr = 0;
  Num_inv = 0;
  numit = 0;
  while (v < Rows) {
    row_nr++;
    if (row_nr > Rows) row_nr = 1;
    v++;
    if (rownum[row_nr - 1] == 1)
      if (frow[row_nr]) {
        v = 0;
        j = Row_end[row_nr - 1] + 1;
        while (!(fcol[Col_no[j] - 1])) j++;
        colnr = Col_no[j];
        fcol[colnr - 1] = FALSE;
        colnum[colnr] = 0;
        for (j = Col_end[colnr - 1]; j < Col_end[colnr]; j++)
          if (frow[Mat[j].row_nr]) rownum[Mat[j].row_nr - 1]--;
        frow[row_nr] = FALSE;
        minoriteration(colnr, row_nr);
      }
  }
  v = 0;
  colnr = 0;
  while (v < Columns) {
    colnr++;
    if (colnr > Columns) colnr = 1;
    v++;
    if (colnum[colnr] == 1)
      if (fcol[colnr - 1]) {
        v = 0;
        j = Col_end[colnr - 1] + 1;
        while (!(frow[Mat[j - 1].row_nr])) j++;
        row_nr = Mat[j - 1].row_nr;
        frow[row_nr] = FALSE;
        rownum[row_nr - 1] = 0;
        for (j = Row_end[row_nr - 1] + 1; j <= Row_end[row_nr]; j++)
          if (fcol[Col_no[j] - 1]) colnum[Col_no[j]]--;
        fcol[colnr - 1] = FALSE;
        numit++;
        col[numit - 1] = colnr;
        row[numit - 1] = row_nr;
      }
  }
  for (j = 1; j <= Columns; j++)
    if (fcol[j - 1]) {
      fcol[j - 1] = FALSE;
      setpivcol(Lower[Rows + j], j + Rows, pcol);
      row_nr = 1;
      while ((!(frow[row_nr] && pcol[row_nr])) && row_nr <= Rows)
        row_nr++; /* this sometimes sets row_nr to Rows + 1 and makes
                     rhsmincol crash. Solved in 2.0? MB */
      if (row_nr == Rows + 1) {
        printf("Inverting failed");
        ERROR(ERR_NUM);
      }
      frow[row_nr] = FALSE;
      condensecol(row_nr, pcol);
      theta = Rhs[row_nr] / (REAL)pcol[row_nr];
      rhsmincol(theta, row_nr, Rows + j);
      addetacol();
    }
  for (i = numit - 1; i >= 0; i--) {
    colnr = col[i];
    row_nr = row[i];
    varin = colnr + Rows;
    for (j = 0; j <= Rows; j++) pcol[j] = 0;
    for (j = Col_end[colnr - 1]; j < Col_end[colnr]; j++) pcol[Mat[j].row_nr] = Mat[j].value;
    pcol[0] -= Extrad;
    condensecol(row_nr, pcol);
    theta = Rhs[row_nr] / (REAL)pcol[row_nr];
    rhsmincol(theta, row_nr, varin);
    addetacol();
  }
  for (i = 1; i <= Rows; i++) my_round(Rhs[i], Epsb);
  if (Lp->print_at_invert)
    fprintf(stderr, "End Invert                eta_size %4d rhs[0] %16.4f\n", Eta_size, -Rhs[0]);

  JustInverted = TRUE;
  DoInvert = FALSE;
  free(rownum);
  free(col);
  free(row);
  free(pcol);
  free(frow);
  free(fcol);
  free(colnum);
} /* invert */

static short colprim(int *colnr, short minit, REAL *drow) {
  int varnr, i, j;
  REAL f, dpiv;

  dpiv = -Epsd;
  (*colnr) = 0;
  if (!minit) {
    for (i = 1; i <= Sum; i++) drow[i] = 0;
    drow[0] = 1;
    btran(drow);
    for (i = 1; i <= Columns; i++) {
      varnr = Rows + i;
      if (!Basis[varnr])
        if (Upbo[varnr] > 0) {
          f = 0;
          for (j = Col_end[i - 1]; j < Col_end[i]; j++) f += drow[Mat[j].row_nr] * Mat[j].value;
          drow[varnr] = f;
        }
    }
    for (i = 1; i <= Sum; i++) my_round(drow[i], Epsd);
  }
  for (i = 1; i <= Sum; i++)
    if (!Basis[i])
      if (Upbo[i] > 0) {
        if (Lower[i])
          f = drow[i];
        else
          f = -drow[i];
        if (f < dpiv) {
          dpiv = f;
          (*colnr) = i;
        }
      }
  if (Lp->trace)
    if ((*colnr) > 0)
      fprintf(stderr, "col_prim:%7d, reduced cost: % 18.10f\n", (*colnr), dpiv);
    else
      fprintf(stderr, "col_prim: no negative reduced costs found, optimality!\n");
  if (*colnr == 0) {
    Doiter = FALSE;
    DoInvert = FALSE;
    Status = OPTIMAL;
  }
  return ((*colnr) > 0);
} /* colprim */

static short rowprim(int colnr, int *row_nr, REAL *theta, REAL *pcol) {
  int i;
  REAL f, quot;

  (*row_nr) = 0;
  (*theta) = Infinite;
  for (i = 1; i <= Rows; i++) {
    f = pcol[i];
    if (my_abs(f) < TREJ) f = 0;
    if (f != 0) {
      quot = 2 * Infinite;
      if (f > 0)
        quot = Rhs[i] / (REAL)f;
      else if (Upbo[Bas[i]] < Infinite)
        quot = (Rhs[i] - Upbo[Bas[i]]) / (REAL)f;
      my_round(quot, Epsel);
      if (quot < (*theta)) {
        (*theta) = quot;
        (*row_nr) = i;
      }
    }
  }
  if ((*row_nr) == 0)
    for (i = 1; i <= Rows; i++) {
      f = pcol[i];
      if (f != 0) {
        quot = 2 * Infinite;
        if (f > 0)
          quot = Rhs[i] / (REAL)f;
        else if (Upbo[Bas[i]] < Infinite)
          quot = (Rhs[i] - Upbo[Bas[i]]) / (REAL)f;
        my_round(quot, Epsel);
        if (quot < (*theta)) {
          (*theta) = quot;
          (*row_nr) = i;
        }
      }
    }

  if ((*theta) < 0) {
    fprintf(stderr, "Warning: Numerical instability, qout = %f\n", (*theta));
    fprintf(stderr, "pcol[%d] = % 18.10f, rhs[%d] = % 18.8f , upbo = % f\n", (*row_nr), f,
            (*row_nr), Rhs[(*row_nr)], Upbo[Bas[(*row_nr)]]);
  }
  if ((*row_nr) == 0) {
    if (Upbo[colnr] == Infinite) {
      Doiter = FALSE;
      DoInvert = FALSE;
      Status = UNBOUNDED;
    } else {
      i = 1;
      while (pcol[i] >= 0 && i <= Rows) i++;
      if (i > Rows) /* empty column with upperbound! */
      {
        Lower[colnr] = FALSE;
        Rhs[0] += 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:%7d, pivot element:% 18.10f\n", (*row_nr), pcol[(*row_nr)]);

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

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

  (*row_nr) = 0;
  minrhs = -Epsb;
  i = 0;
  artifs = FALSE;
  while (i < Rows && !artifs) {
    i++;
    f = Upbo[Bas[i]];
    if (f == 0 && (Rhs[i] != 0)) {
      artifs = TRUE;
      (*row_nr) = i;
    } else {
      if (Rhs[i] < f - Rhs[i])
        g = Rhs[i];
      else
        g = f - Rhs[i];
      if (g < minrhs) {
        minrhs = g;
        (*row_nr) = i;
      }
    }
  }

  if (Lp->trace) {
    if ((*row_nr) > 0) {
      fprintf(stderr, "row_dual:%7d, rhs of selected row:           % 18.10f\n", (*row_nr),
              Rhs[(*row_nr)]);
      if (Upbo[Bas[(*row_nr)]] < Infinite)
        fprintf(stderr, "               upper bound of basis variable:    % 18.10f\n",
                Upbo[Bas[(*row_nr)]]);
    } else
      fprintf(stderr, "row_dual: no infeasibilities found\n");
  }

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

static short coldual(int row_nr, int *colnr, short minit, REAL *prow, REAL *drow) {
  int i, j, r, varnr;
  REAL theta, quot, pivot, d, f, g;

  Doiter = FALSE;
  if (!minit) {
    for (i = 0; i <= Rows; i++) {
      prow[i] = 0;
      drow[i] = 0;
    }
    drow[0] = 1;
    prow[row_nr] = 1;
    for (i = Eta_size; i >= 1; i--) {
      d = 0;
      f = 0;
      r = Eta_row_nr[Eta_col_end[i] - 1];
      for (j = Eta_col_end[i - 1]; j < Eta_col_end[i]; j++) {
        /* this is where the program consumes most cpu time */
        f += prow[Eta_row_nr[j]] * Eta_value[j];
        d += drow[Eta_row_nr[j]] * Eta_value[j];
      }
      my_round(f, Epsel);
      prow[r] = f;
      my_round(d, Epsel);
      drow[r] = d;
    }
    for (i = 1; i <= Columns; i++) {
      varnr = Rows + i;
      if (!Basis[varnr]) {
        d = -Extrad * drow[0];
        f = 0;
        for (j = Col_end[i - 1]; j < Col_end[i]; j++) {
          d = d + drow[Mat[j].row_nr] * Mat[j].value;
          f = f + prow[Mat[j].row_nr] * Mat[j].value;
        }
        drow[varnr] = d;
        prow[varnr] = f;
      }
    }
    for (i = 0; i <= Sum; i++) {
      my_round(prow[i], Epsel);
      my_round(drow[i], Epsd);
    }
  }
  if (Rhs[row_nr] > Upbo[Bas[row_nr]])
    g = -1;
  else
    g = 1;
  pivot = 0;
  (*colnr) = 0;
  theta = Infinite;
  for (i = 1; i <= Sum; i++) {
    if (Lower[i])
      d = prow[i] * g;
    else
      d = -prow[i] * g;
    if ((d < 0) && (!Basis[i]) && (Upbo[i] > 0)) {
      if (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:%7d, pivot element:  % 18.10f\n", (*colnr), prow[(*colnr)]);

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

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

static void iteration(int row_nr, int varin, REAL *theta, REAL up, short *minit, short *low,
                      short primal, REAL *pcol) {
  int i, k, varout;
  REAL f;
  REAL pivot;

  Lp->iter++;
  (*minit) = (*theta) > (up + Epsb);
  if ((*minit)) {
    (*theta) = up;
    (*low) = !(*low);
  }
  k = Eta_col_end[Eta_size + 1];
  pivot = Eta_value[k - 1];
  for (i = Eta_col_end[Eta_size]; i < k; i++) {
    f = Rhs[Eta_row_nr[i]] - (*theta) * Eta_value[i];
    my_round(f, Epsb);
    Rhs[Eta_row_nr[i]] = f;
  }
  if (!(*minit)) {
    Rhs[row_nr] = (*theta);
    varout = Bas[row_nr];
    Bas[row_nr] = varin;
    Basis[varout] = FALSE;
    Basis[varin] = TRUE;
    if (primal && pivot < 0) Lower[varout] = FALSE;
    if (!(*low) && up < Infinite) {
      (*low) = TRUE;
      Rhs[row_nr] = up - Rhs[row_nr];
      for (i = Eta_col_end[Eta_size]; i < k; i++) Eta_value[i] = -Eta_value[i];
    }
    addetacol();
    Num_inv++;
  }
  if (Lp->trace) {
    fprintf(stderr, "Theta = %16.4g ", (*theta));
    if ((*minit)) {
      if (!Lower[varin])
        fprintf(stderr, "Iteration:%6d, variable%5d changed from 0 to its upper bound of %12f\n",
                Lp->iter, varin, Upbo[varin]);
      else
        fprintf(stderr, "Iteration:%6d, variable%5d changed its upper bound of %12f to 0\n",
                Lp->iter, varin, Upbo[varin]);
    } else
      fprintf(stderr, "Iteration:%6d, variable%5d entered basis at:% 18.10f\n", Lp->iter, varin,
              Rhs[row_nr]);
    if (!primal) {
      f = 0;
      for (i = 1; i <= Rows; i++)
        if (Rhs[i] < 0)
          f -= Rhs[i];
        else if (Rhs[i] > Upbo[Bas[i]])
          f += Rhs[i] - Upbo[Bas[i]];
      fprintf(stderr, "feasibility gap of this basis:% 18.10f\n", (double)f);
    } else
      fprintf(stderr, "objective function value of this feasible basis: % 18.10f\n",
              (double)Rhs[0]);
  }
} /* iteration */

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

  CALLOC(drow, Sum + 1, REAL);
  CALLOC(prow, Sum + 1, REAL);
  CALLOC(Pcol, Rows + 1, REAL);
  CALLOC(test, Sum + 1, short);

  Lp->iter = 0;
  minit = FALSE;
  Status = RUNNING;
  DoInvert = FALSE;
  Doiter = FALSE;
  i = 0;
  primal = TRUE;
  while (i != Rows && primal) {
    i++;
    primal = Rhs[i] >= 0 && Rhs[i] <= Upbo[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 <= Rows; i++) drow[i] = 0;
    Extrad = 0;
    for (i = 1; i <= Columns; i++) {
      varnr = Rows + i;
      drow[varnr] = 0;
      for (j = Col_end[i - 1]; j < Col_end[i]; j++)
        if (drow[Mat[j].row_nr] != 0) drow[varnr] += drow[Mat[j].row_nr] * Mat[j].value;
      if (drow[varnr] < Extrad) Extrad = drow[varnr];
    }
  } else
    Extrad = 0;
  if (Lp->trace) fprintf(stderr, "Extrad = %f\n", Extrad);
  minit = FALSE;

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

    if (primal) {
      if (colprim(&colnr, minit, drow)) {
        setpivcol(Lower[colnr], colnr, Pcol);
        if (rowprim(colnr, &row_nr, &theta, Pcol)) condensecol(row_nr, Pcol);
      }
    } else /* not primal */
    {
      if (!minit) rowdual(&row_nr);
      if (row_nr > 0) {
        if (coldual(row_nr, &colnr, minit, prow, drow)) {
          setpivcol(Lower[colnr], colnr, Pcol);
          /* getting div by zero here ... MB */
          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, "Reinverting Eta\n");
              DoInvert = TRUE;
            } else {
              fprintf(stderr, "Can't reinvert, failure\n");
              Status = FAILURE;
            }
          } else {
            condensecol(row_nr, Pcol);
            f = Rhs[row_nr] - Upbo[Bas[row_nr]];
            if (f > 0) {
              theta = f / (REAL)Pcol[row_nr];
              if (theta <= Upbo[colnr]) Lower[Bas[row_nr]] = !Lower[Bas[row_nr]];
            } else /* f <= 0 */
              theta = Rhs[row_nr] / (REAL)Pcol[row_nr];
          }
        } else
          Status = INFEASIBLE;
      } else {
        primal = TRUE;
        Doiter = FALSE;
        Extrad = 0;
        DoInvert = TRUE;
      }
    }
    if (Doiter) iteration(row_nr, colnr, &theta, Upbo[colnr], &minit, &Lower[colnr], primal, Pcol);
    if (Num_inv >= Lp->max_num_inv) DoInvert = TRUE;
    if (DoInvert) {
      if (Lp->print_at_invert) fprintf(stderr, "Inverting: Primal = %d\n", primal);
      invert();
    }
  }

  Lp->total_iter += Lp->iter;

  free(drow);
  free(prow);
  free(Pcol);
  free(test);

  return (Status);
} /* solvelp */

static short is_int(REAL value) {
  REAL tmp;

  tmp = value - (REAL)floor((double)value);
  if (tmp < Epsilon) return (TRUE);
  if (tmp > (1 - Epsilon)) return (TRUE);
  return (FALSE);
} /* is_int */

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

  /* zero all results of rows */
  memset(sol, '\0', (Rows + 1) * sizeof(REAL));

  if (Lp->scaling_used) {
    for (i = Rows + 1; i <= Sum; i++) sol[i] = Lowbo[i] * Lp->scale[i];
    for (i = 1; i <= Rows; i++) {
      basi = Bas[i];
      if (basi > Rows) sol[basi] += Rhs[i] * Lp->scale[basi];
    }
    for (i = Rows + 1; i <= Sum; i++)
      if (!Basis[i] && !Lower[i]) sol[i] += Upbo[i] * Lp->scale[i];

    for (j = 1; j <= Columns; j++) {
      f = sol[Rows + j];
      if (f != 0)
        for (i = Col_end[j - 1]; i < Col_end[j]; i++)
          sol[Mat[i].row_nr] +=
              (f / Lp->scale[Rows + j]) * (Mat[i].value / Lp->scale[Mat[i].row_nr]);
    }

    for (i = 0; i <= Rows; i++) {
      if (my_abs(sol[i]) < Epsb)
        sol[i] = 0;
      else if (Lp->ch_sign[i])
        sol[i] = -sol[i];
    }
  } else {
    for (i = Rows + 1; i <= Sum; i++) sol[i] = Lowbo[i];
    for (i = 1; i <= Rows; i++) {
      basi = Bas[i];
      if (basi > Rows) sol[basi] += Rhs[i];
    }
    for (i = Rows + 1; i <= Sum; i++)
      if (!Basis[i] && !Lower[i]) sol[i] += Upbo[i];
    for (j = 1; j <= Columns; j++) {
      f = sol[Rows + j];
      if (f != 0)
        for (i = Col_end[j - 1]; i < Col_end[j]; i++) sol[Mat[i].row_nr] += f * Mat[i].value;
    }

    for (i = 0; i <= Rows; i++) {
      if (my_abs(sol[i]) < Epsb)
        sol[i] = 0;
      else if (Lp->ch_sign[i])
        sol[i] = -sol[i];
    }
  }
} /* construct_solution */

static void calculate_duals(void) {
  int i;

  /* initialise */
  for (i = 1; i <= Rows; i++) Lp->duals[i] = 0;
  Lp->duals[0] = 1;
  btran(Lp->duals);
  if (Lp->scaling_used)
    for (i = 1; i <= 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. Can this happen? */
  for (i = 1; i <= Rows; i++) {
    if (Lp->basis[i])
      Lp->duals[i] = 0;
    else if (Lp->ch_sign[0] == Lp->ch_sign[i])
      Lp->duals[i] = -Lp->duals[i];
  }
}

int milpsolve(sstate *st, REAL *upbo, REAL *lowbo, short *sbasis, short *slower, int *sbas) {
  int i, j, failure, notint, is_worse;
  REAL theta, tmpreal;

  /* First, check for "time-out", time to give control back to lisp */
  if (st->next)
    st = st->next;
  else
    for (i = 0; i < 20; i++) { /* Need more state-saving space, so create 20 more. */
      sstate *newst;
      newst = (sstate *)malloc(sizeof(sstate));
      if (!st) ERROR(ERR_ST);
      newst->next = st->next;
      st->next = newst;
      newst->saved = newst->notint = 0;
    }

  if (st->saved) {
    if (st->saved == ST_SOLN) {
      st->saved = 0;
      return (OPTIMAL);
    }
  } else if (SolveCount++ > 100)
    return (TIMEOUT); /* Time out every 100 LP solves */
  else if ((KBDEventFlg > 0) && *KEYBUFFERING68k == ATOM_T)
    return (TIMEOUT); /* Time out on key/mouse clicks */
#ifdef DOS
  else if (currentmouse->Cursor.Moved)
    return (TIMEOUT); /* Time out if mouse moves in DOS */
#endif                /* DOS */

  if (Break_bb) return (BREAK_BB);
  Level++;
  Lp->total_nodes++;
  if (Level > Lp->max_level) Lp->max_level = Level;

  if (!st->saved) { /* We're coming into this fresh, rather than returnin from a TIMEOUT. */
    /* make fresh copies of upbo, lowbo, rh as solving changes them */
    memcpy(Upbo, upbo, (Sum + 1) * sizeof(REAL));
    memcpy(Lowbo, lowbo, (Sum + 1) * sizeof(REAL));
    memcpy(Basis, sbasis, (Sum + 1) * sizeof(short));
    memcpy(Lower, slower, (Sum + 1) * sizeof(short));
    memcpy(Bas, sbas, (Rows + 1) * sizeof(int));
    memcpy(Rh, Orig_rh, (Rows + 1) * sizeof(REAL));

    if (Lp->anti_degen) {
      for (i = 1; i <= Columns; i++) {
        tmpreal = (REAL)(rand() % 100 * 0.00001);
        if (tmpreal > Epsb) Lowbo[i + Rows] -= tmpreal;
        tmpreal = (REAL)(rand() % 100 * 0.00001);
        if (tmpreal > Epsb) Upbo[i + Rows] += tmpreal;
      }
      Lp->eta_valid = FALSE;
    }

    if (!Lp->eta_valid) {
      for (i = 1; i <= Columns; i++)
        if (Lowbo[Rows + i] != 0) {
          theta = Lowbo[Rows + i];
          if (Upbo[Rows + i] < Infinite) Upbo[Rows + i] -= theta;
          for (j = Col_end[i - 1]; j < Col_end[i]; j++) Rh[Mat[j].row_nr] -= theta * Mat[j].value;
        }
      invert();
      Lp->eta_valid = TRUE;
    }

    failure = solvelp();

    if (Lp->anti_degen) {
      memcpy(Upbo, upbo, (Sum + 1) * sizeof(REAL));
      memcpy(Lowbo, lowbo, (Sum + 1) * sizeof(REAL));
      memcpy(Rh, Orig_rh, (Rows + 1) * sizeof(REAL));

      for (i = 1; i <= Columns; i++)
        if (Lowbo[Rows + i] != 0) {
          theta = Lowbo[Rows + i];
          if (Upbo[Rows + i] < Infinite) Upbo[Rows + i] -= theta;
          for (j = Col_end[i - 1]; j < Col_end[i]; j++) Rh[Mat[j].row_nr] -= theta * Mat[j].value;
        }
      invert();
      Lp->eta_valid = TRUE;
      failure = solvelp();
    }

    if (failure == INFEASIBLE && Lp->verbose) fprintf(stderr, "level%4d INF\n", Level);
  } else
    failure =
        OPTIMAL; /* Coming back thru after a timeout; we got OPTIMAL last time, so do it again. */

  if (failure == OPTIMAL) /* there is a solution */
  {
    if (!st->saved) {
      construct_solution(Solution);

      /* if this solution is worse than the best sofar, this branch must die */
      if (Maximise)
        is_worse = Solution[0] <= Best_solution[0];
      else /* minimising! */
        is_worse = Solution[0] >= Best_solution[0];

      if (is_worse) {
        if (Lp->verbose)
          fprintf(stderr, "level%4d OPT NOB value %f bound %f\n", Level, Solution[0],
                  Best_solution[0]);
        Level--;
        return (MILP_FAIL);
      }

      /* check if solution contains enough ints */
      st->notint = notint = 0;
      if (Lp->bb_rule == FIRST_NI) {
        notint = 0;
        i = Rows + 1;
        while (i <= Sum && notint == 0) {
          if (Must_be_int[i] && !is_int(Solution[i]))
            if (lowbo[i] == upbo[i]) /* this var is already fixed */
            {
              fprintf(
                  stderr,
                  "Warning: integer var %d is already fixed at %d, but has non-integer value %g\n",
                  i - Rows, (int)lowbo[i], Solution[i]);
              fprintf(stderr, "Perhaps the -e option should be used\n");
            } else
              st->notint = notint = i;
          i++;
        }
      }
      if (Lp->bb_rule == RAND_NI) {
        int nr_not_int, select_not_int;
        nr_not_int = 0;
        for (i = Rows + 1; i <= Sum; i++)
          if (Must_be_int[i] && !is_int(Solution[i])) nr_not_int++;
        if (nr_not_int == 0)
          notint = 0;
        else {
          select_not_int = (rand() % nr_not_int) + 1;
          i = Rows + 1;
          while (select_not_int > 0) {
            if (Must_be_int[i] && !is_int(Solution[i])) select_not_int--;
            i++;
          }
          st->notint = notint = i - 1;
        }
      }
    } else
      notint = st->notint; /* Coming back in, use old value. */

    if (Lp->verbose == TRUE)
      if (notint)
        fprintf(stderr, "level %3d OPT     value %f\n", Level, Solution[0]);
      else
        fprintf(stderr, "level %3d OPT INT value %f\n", Level, Solution[0]);

    if (notint) /* there is at least one value not yet int */
    {
      /* set up two new problems */
      REAL *new_upbo, *new_lowbo;
      REAL new_bound;
      short *new_lower, *new_basis;
      int *new_bas;
      int resone, restwo;

      /* allocate room for them */
      MALLOC(new_upbo, Sum + 1, REAL);
      MALLOC(new_lowbo, Sum + 1, REAL);
      MALLOC(new_lower, Sum + 1, short);
      MALLOC(new_basis, Sum + 1, short);
      MALLOC(new_bas, Rows + 1, int);
      memcpy(new_upbo, upbo, (Sum + 1) * sizeof(REAL));
      memcpy(new_lowbo, lowbo, (Sum + 1) * sizeof(REAL));
      memcpy(new_lower, Lower, (Sum + 1) * sizeof(short));
      memcpy(new_basis, Basis, (Sum + 1) * sizeof(short));
      memcpy(new_bas, Bas, (Rows + 1) * sizeof(int));

      if (Floor_first) {
        if (st->saved)
          new_bound = st->bound;
        else
          st->bound = new_bound = ceil(Solution[notint]) - 1;
        /* this bound might conflict */
        if (st->saved >= ST_HI)
          resone = st->res1; /* We got the upper bound earlier; skip lower */
        else if (new_bound < lowbo[notint]) {
          resone = MILP_FAIL;
        } else /* bound feasible */
        {
          new_upbo[notint] = new_bound;
          Lp->eta_valid = FALSE;
          st->saved = ST_LO;
          resone = milpsolve(st, new_upbo, lowbo, new_basis, new_lower, new_bas);
          Lp->eta_valid = FALSE;
          st->res1 = resone;
          if ((resone == INT_SOLN) || (resone == TIMEOUT)) {
            Level--;
            free(new_upbo);
            free(new_lowbo);
            free(new_lower);
            free(new_basis);
            free(new_bas);
            return resone;
          }
        }
        new_bound += 1;
        if (new_bound > upbo[notint]) {
          restwo = MILP_FAIL;
        } else /* bound feasible */
        {
          new_lowbo[notint] = new_bound;
          st->saved = ST_HI;
          Lp->eta_valid = FALSE;
          restwo = milpsolve(st, upbo, new_lowbo, new_basis, new_lower, new_bas);
          Lp->eta_valid = FALSE;
          st->res2 = restwo;
          if ((restwo == INT_SOLN) || (restwo == TIMEOUT)) {
            Level--;
            free(new_upbo);
            free(new_lowbo);
            free(new_lower);
            free(new_basis);
            free(new_bas);
            return restwo;
          }
        }
      } else /* take ceiling first */
      {
        new_bound = ceil(Solution[notint]);
        /* this bound might conflict */
        if (new_bound > upbo[notint]) {
          resone = MILP_FAIL;
        } else /* bound feasible */
        {
          new_lowbo[notint] = new_bound;
          Lp->eta_valid = FALSE;
          resone = milpsolve(st, upbo, new_lowbo, new_basis, new_lower, new_bas);
          Lp->eta_valid = FALSE;
        }
        new_bound -= 1;
        if (new_bound < lowbo[notint]) {
          restwo = MILP_FAIL;
        } else /* bound feasible */
        {
          new_upbo[notint] = new_bound;
          Lp->eta_valid = FALSE;
          restwo = milpsolve(st, new_upbo, lowbo, new_basis, new_lower, new_bas);
          Lp->eta_valid = FALSE;
        }
      }
      if (resone && restwo) /* both failed and must have been infeasible */
        failure = INFEASIBLE;
      else
        failure = OPTIMAL;

      free(new_upbo);
      free(new_lowbo);
      free(new_basis);
      free(new_lower);
      free(new_bas);
    } else /* all required values are int */
    {
      if (Maximise)
        is_worse = Solution[0] < Best_solution[0];
      else
        is_worse = Solution[0] > Best_solution[0];

      if (!is_worse) /* Current solution better */
      {
        if (Lp->debug || (Lp->verbose && !Lp->print_sol))
          fprintf(stderr, "*** new best solution: old: %g, new: %g ***\n", (double)Best_solution[0],
                  (double)Solution[0]);
        memcpy(Best_solution, Solution, (Sum + 1) * sizeof(REAL));
        calculate_duals();
        if (Lp->print_sol) print_solution(Lp);
        if (Lp->break_at_int) {
          if (Maximise && (Best_solution[0] > Lp->break_value)) Break_bb = TRUE;
          if (!Maximise && (Best_solution[0] < Lp->break_value)) Break_bb = TRUE;
        }
        st->saved = INT_SOLN; /* Tell caller we found -a- solution */
        failure = INT_SOLN;   /* & remember that fact for next time. */
      }
    }
  }

  /* failure can have the values OPTIMAL, UNBOUNDED and INFEASIBLE. */
  /*                             INT_SOLN, TIMEOUT, or MILP_FAIL */
  Level--;
  st->saved = 0; /* We're done at this level, so mark the ST empty. */
  return (failure);
} /* milpsolve */

int solve(lprec *lp) {
  int result, i;

  if (!lp->active) set_globals(lp);

  lp->total_iter = 0;
  lp->max_level = 1;
  lp->total_nodes = 0;

  if (Isvalid(lp)) {
    if (Maximise && lp->obj_bound == Infinite)
      Best_solution[0] = -Infinite;
    else if (!Maximise && lp->obj_bound == -Infinite)
      Best_solution[0] = Infinite;
    else
      Best_solution[0] = lp->obj_bound;

    Level = 0;

    if (!lp->basis_valid) {
      for (i = 0; i <= lp->rows; i++) {
        lp->basis[i] = TRUE;
        lp->bas[i] = i;
      }
      for (i = lp->rows + 1; i <= lp->sum; i++) lp->basis[i] = FALSE;
      for (i = 0; i <= lp->sum; i++) lp->lower[i] = TRUE;
      lp->basis_valid = TRUE;
    }

    lp->eta_valid = FALSE;
    Break_bb = FALSE;
    result = milpsolve(lp->solve_states, Orig_upbo, Orig_lowbo, Basis, Lower, Bas);
    lp->eta_size = Eta_size;
    lp->eta_alloc = Eta_alloc;
    lp->num_inv = Num_inv;
  }
  return (result);
}

/* * * * lag_solve used to be here * * * * */
