github.com/Interlisp.maiko
1
0
Fork 0
mirror of https://github.com/Interlisp/maiko.git synced 2026-01-15 15:57:13 +00:00
Code Issues Projects Releases Wiki Activity
Interlisp.maiko/src/lpsolve.c
Nick Briggs 6fb4b9a189
Remove sccs id lines (#98) 
* Remove static char *id = from all source files.

The same information is included in a comment in each source file.

* Remove unused template file 'id'
2020-12-19 19:08:52 -08:00

1248 lines
34 KiB
C
Raw Blame History

/* $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 * * * * */
Reference in New Issue View Git Blame Copy Permalink