/* blp.c  BLP program */

/* 10/2001: fixed marginal-value sign on = constraints */


#include <stdio.h>
#include <stdlib.h>
#include "blp.h"
#include <math.h>
/*
    8/26/91 Repaired range_cb()
*/
#define MAX_VARS 	50
#define MAX_ROWS 	30
#define NSTATES 	7+1
#define SIZE_LHS 	57
#define SIZE_NAME 	8
#define SIZE_LONG_NAME	SIZE_NAME+4
#define SIZE_LABEL	SIZE_NAME+1
#define SIZE_RHS 	10
#define DIM_ROWS 	(MAX_ROWS+2)
#define DIM_COLS 	(MAX_ROWS+MAX_VARS+1)
#define TAB 		9

#define BLANK   ' '
#define CR   	'\n'   /* RSB 5/16/06 
#define CR   	0x0d */
#define EOS  	'\0'
#define FALSE 	0
#define PERIOD 	'.'
#define TRUE 	!FALSE
#define OK	TRUE

typedef double REAL;

REAL 	a[MAX_VARS+1];
double 	a_mat[DIM_ROWS][DIM_COLS];
int 	basis_index[DIM_ROWS];           /* column number, 0=artif */
char    buf[100];                        /* Temporary sprint buffer */
char    chtype[256];
int 	col_status[DIM_COLS];            /* +=basic, 0= @lb, -= @ub */
double	cost[DIM_COLS];
char	isart[DIM_COLS];		/* TRUE=artificial column */
char	lhs[SIZE_LHS+1];
double	lowerb[DIM_COLS];
char	minmax[12];
char 	name[SIZE_NAME+1];
char	rel[1+1];
char	rhs[10+1];
char	rlhs[SIZE_LHS+1];
char	rname[DIM_ROWS][SIZE_LABEL];
double	rrhs[DIM_ROWS];
char	rrel[DIM_ROWS][3];
int 	row_type[DIM_ROWS];              /* 1=LE, 0=EQ, -1=GE */
REAL	save_rhs[DIM_ROWS];
char	titlecons[50];
char 	titlevar[50];
double 	upperb[DIM_COLS];
char    varname[MAX_VARS+1][SIZE_LONG_NAME+1];


double 	adj;
double	direction;			/* -1 = max, +1 = min */
double	element;
double	elt;
double 	eps_infe;
double	eps_pivot;
double	eps_price;
double	eps_print = 0.000001;		/* epsilon for printing model */
double	eps_ratio;
double	eps_ub;
double 	infinity;
double	min_cbar;
double	min_ratio;
REAL 	objective = 1.0;
double	piv_elem;
double	ratio;
double 	slk;
double	sum;
double	x;
double	z_1;
double	z_mult;


int 	artif;
int	best_col;
int	best_row;
int	bounded = FALSE;			/* Bounded variables used? */
int	check;
int	col;
int	debug 		= FALSE;		/* Display tableaus option? */
int	dispModel	= TRUE;			/* Display model option? */
int	echo		= TRUE;
int	end_stat;
int 	err_code;
int	errCol;					/* Input column of error */
int	fixed_col;
int	i_last;
int	infeasible;
int	inputMatrix = FALSE;			/* Input model as matrix? */
int	last_col;
int 	last_row;
int	last_tableau = FALSE;
int	leav_var;
int	max_cons;
int	max_vars;
int     nnames;
int	nrpt = -1;
int	num_artif;
int 	num_cons;
int	num_lbs;
int	num_vars;
int	obj_row;
int	optimal;
int	out_col;
int 	phase;
int	ph1_row;
int	ph2_row;
int	piv_col;
int	piv_row;
int	pivots;
int	ratio_row;
int 	rhs_col;
int	row;
int	show_exchanges = FALSE;
int	show_range = FALSE;
int	s_col;
int	tokLen;					/* # chars in input buffer */
int	unbounded;
int 	wbase;
int	wrpt[2];
int	wcons;
int	wvar;


FILE *fopen(),*luin,*luout;
double atof();

char lhstate[][NSTATES] =
    {
/*   Len Alph Dig Pd  +- Blk Rel   )                       */
    {  1,  1,  2,  3,  2,  0,  4,  5 },   /* State 0: Start  */
    {  SIZE_NAME,
	   1,  1,  1, -1, -1, -1,  5 },   /* State 1: Name   */
    { 99, -1,  2,  3, -1,  2, -1, -1 },   /* State 2: Integer */
    { 99, -1,  3, -2, -1,  3, -1, -1 },    /* State 3: Fraction */
    {  2, -1, -1, -1, -1,  4,  4, -1 },	  /* State 4: Relation */
    {  SIZE_LABEL,
	  -1, -1, -1, -1, -1, -1, -1 }	  /* State 5: Label */
    };

char tokState[][6] =
	{
/*		Name	Int	Fract	Relat	Label	*/
	{ 0,	1,	-2,	-2,	-2,	-2 },	/* 0. Start */
	{ 0,	-3,	-3,	-3,	-3,	 2 },	/* 1. Min or MAX */
	{ 0,	4,	3,	3,	-4,	-5 },	/* 2. Obj label */
	{ 0,	4,	-6,	-6,	-7,	-7 },	/* 3. Obj coeff */
	{ 0,	-8,	3,	3,	-4,	 5 },	/* 4. Obj variable */
	{ 0,	 7,	 6,	 6,	-9,	-10},	/* 5. Row label */
	{ 0,	 7,	-6,	-6,	-7,	-11},	/* 6. Coeff */
	{ 0,	-8,	 6,	 6,	 8,	-11},	/* 7. Variable */
	{ 0,	-12,	 9,	 9,	-13,	-14},	/* 8. Relation */
	{ 0,	-15,	-15,	-15,	-15,	 5 }	/* 9. RHS */
	};

char	*errorMsgs[] = {
	"Syntax error",
	"Syntax error on input",
	"Begin a problem with 'MIN' or 'MAX'",		/* 2 */
	"Objective function label missing",
	"Relation (<,=,>) not used in objective function",
	"Objective function missing",			/* 5 */
	"Two numbers may not be adjacent",
	"No variable for coefficient",
	"Two variable names may not be adjacent",
	"No left-hand-side for constraint",
	"Two labels may not be adjacent",		/* 10 */
	"Label unexpectedly encountered",
	"Righthand-side not a nonnegative constant",
	"Two relations (<,=,>) may not be adjacent",
	"Righthand side missing",
	"Constraint label missing",			/* 15 */
	"Maximum number of constraints exceeded",
	"Maximum number of variables exceeded",
	"Element too long.\n(Variable and row names limited to 8 characters.)",
	"Duplicate row name",
	"Label must begin with letter"			/* 20 */
	};

struct TERMS
    {
    double  tvalue;             /* linear multiplier for this term */
    int     tid;                /* Name id # */
    };

struct  TERMS term[MAX_VARS];
int     nterms;

main(argc,argv)
int argc; char *argv[];
	{
	if(getArg(argc,argv) != OK) return;
	initParam();			/* Prepare arrays for read */
	if(inputMatrix) generate2();
	else 		generate3();
	setprob();
	/*****************
	if (dispModel) showModel();
	 * Print tableau and stop
	 *
	simplex();
	if (end_stat == 1)	out_opt();
	else if (end_stat == 2)	out_infe();
	else			out_unbnd();
	 *****************/
	 tableau();
	}

parse()
    {
	generate3();
	setprob();
	showModel();
    }

generate3()
	{
	register int state, tokType, i;
	double	coeff, tokValue();

	chset();

	for (state=0, buf[0]=EOS; TRUE;)
		{
		if((tokType = getTok()) < 0) break;
		if((state = tokState[state][tokType]) < 0) errExit(-state);
		switch (state)
			{
			case 1:		/* Min/max */
				if (strnicmp("max",buf,3) == 0) z_mult = -1.0;
				else if (strnicmp("min",buf,3) == 0) z_mult = 1.0;
				else errExit(2);
				strncpy(minmax,strupr(buf),3);
				objective	= z_mult;
				break;
			case 2:		/* Obj label */
				num_cons	= -1;
			case 5:		/* Row label */
				if (++num_cons > MAX_ROWS) errExit(16);
				coeff		= 1.0;
				buf[--tokLen]	= EOS;
				if(tokLen == 0) errExit(20);
				strncpy(&rname[num_cons][0],buf,SIZE_LABEL);
				if (checkLabel() != OK) errExit(19);
				break;
			case 3:		/* obj coeff */
			case 6:		/* row coeff */
				coeff	= tokValue();
				break;
			case 4:		/* Obj variable */
			case 7:		/* Row variable */
				if ((i = nameNum(buf)) < 0) errExit(17);
				else	a_mat[num_cons][i] = coeff;
				break;
			case 8:		/* relation */
				row_type[num_cons] =
					(buf[0]=='<') ? 1 : ((buf[0]=='>')? -1 : 0);
				break;
			case 9:		/* RHS */
				a_mat[num_cons][0] = rrhs[num_cons] = tokValue();
				if (rrhs[num_cons] < 0) errExit(12);
				break;
			}
		}
	last_row	= num_cons;
	last_col	= num_vars	= nnames;
	}

getTok()
	/* Gets next token from luin.  Returns:
	 *	>0 = token type from lhstate[][] table or
	 *	-2 = input error at column errCol
	 *	-3 = if EOF
	 * Character version of token is in buf[], length tokLen, null terminated.
	 * Input is echoed to stdout.
	 */

	{
	int	inch, newstate, state;

	for(newstate = tokLen = 0; newstate >= 0;)
		{
		state	= newstate;
		if ((inch = getc(luin)) == EOF)	newstate = -3;
		else	newstate = lhstate[state][chtype[inch]];

		if (newstate>=0)
			{
			if (echo) putchar(inch);
			if (inch == CR) errCol = 0;
			else	{
				++errCol;      /* Don't save white space */
				if(chtype[inch] != 5 ) 
				  buf[tokLen++] = Toupper(inch);
				}
			}
		}
	switch (newstate)
		{
		case -1:	ungetc(inch,luin);
				buf[tokLen] = EOS; /* Check length */
				if (tokLen > lhstate[state][0]) errExit(18);
				return state;
		case -2:	printf("Error in column %d\n",errCol);
				errExit(1);
		case -3:	if (state == 0) return newstate;
				buf[tokLen] = EOS;
				return state;
		}
	}

checkLabel()
	{
	register int row;
	if(num_cons == 0) return OK;
	for (row=0; row<num_cons; ++row)
		if (strncmp(&rname[row][0],&rname[num_cons][0],SIZE_LABEL)
			== 0) return !OK;
	return OK;
	}

errExit(i)
int i;
	{
	register int k;

	printf("\n");
	for (k=1; k < errCol; ++k) putchar(' ');
	printf("^\nInput Error %d>> %s\nProcessing terminated\n",i,errorMsgs[i]);
	exit(-1);
	}

double tokValue()
	{
	switch (tokLen)
		{
		case 0:	return 1.0;
		case 1: switch (buf[0])
				{
				case '+':
				case '.':	return 1.0;
				case '-':	return -1.0;
				}
			break;
		case 2: if (!strncmp(buf, "+.", 2)) return 1.0;
			if (!strncmp(buf, "-.", 2)) return -1.0;
	}
    return atof(buf);
	}

nameNum()
    {
    register int i;
    if ((i = findname(buf)) < 0) i = putname(buf);
    return i;
    }

findname(name)
char *name;
    {
    register int i;
    for (i=1; i <= nnames; ++i)
	if (strncmp(name, &varname[i][0], SIZE_LONG_NAME+1) == 0) return i;

    return -1;
    }

putname(name)
char *name;
    {
    if (nnames >= MAX_VARS) return -1;
    else
	{
	strncpy(&varname[++nnames][0], name, SIZE_LONG_NAME+1);
	return nnames;
	}
    }

getArg(argc,argv)
int argc; char *argv[];
	{
	register int i, j;
	register char *ch; 
	char buf2[80];

	if (argc == 1) { help(); return !OK; }
	for(i=1; i<argc; ++i)
		{
		ch	= argv[i];
		if (*ch == '-')
			{
			while(*(++ch))
				switch (Tolower(*ch)) {
				case 'e': echo		= !echo;break;
				case 'i': inputMatrix	= TRUE; break;
				case 'f': last_tableau 	= TRUE; break;
				case 'm': dispModel	= !dispModel; break;
				case 'r': show_range 	= TRUE; break;
				case 't': debug 	= TRUE; break;
				case 'x': show_exchanges= TRUE; break;
				default:  printf("Unknown option '%c'\n",*ch);
				}
			}
		else strcpy(buf,argv[i]);
		}

	if ((luin = fopen(buf,"r")) == NULL)
		{
		printf("Data file %s not found\n",buf);
		return !OK;
		}

	ctitle("*** BLP - LINEAR PROGRAMMING ***",70);
	/*********
	 sprintf(buf2,"*** DATA FILE %s ***",strupr(buf));
	if (echo) ctitle(buf2, 70);
	**********/
	return OK;
	}

solve()
    {
    generate2();
    setprob();
    if(debug) tableau();
    simplex();
    if (end_stat == 1) out_opt();
    else if (end_stat == 2) out_infe();
    else out_unbnd();                               /*write outfile*/
    }

initParam()
	/*     i n i t i a l i z e    p a r a m e t e r s*/
    {
	register int i,j;

	max_vars 	= (MAX_VARS);
	max_cons 	= (MAX_ROWS);
	infinity	= 1E30;
	unbounded	= FALSE;
	artif 		= 0;                                /*artificials "column"*/
	rhs_col 	= 0;                                  /*amat column for rhs*/
	eps_infe 	= .00001;
	eps_ratio 	= .000000000001;
	eps_pivot 	= .00001;
	eps_price 	= .00001;
	eps_ub 		= .00001;
	fixed_col 	= 10000;                        /*fixed column status*/
	upperb[artif] 	= infinity;
	num_lbs 	= 0;

	for (j=0; j <= DIM_COLS; ++j)
		{
		cost[j] 	= 0.0;
		upperb[j] 	= infinity;
		lowerb[j] 	= 0.0;
		isart[j]	= FALSE;
		for (i=0; i <= num_cons; ++i) a_mat[i][j] = 0.0;
		}

    }

generate2()
    {
    register int i, j;
    int dummy, error, k;
    char ch;

    num_cons 	= last_row;
    nnames 	= 0;         /* Start new list of names */
    z_mult 	= objective;

    fscanf(luin,"%s %d %d",minmax, &num_vars, &num_cons);
    z_mult = objective = (strnicmp("max",minmax,3) == 0) ? -1.0 : 1.0;
    printf("Processing %d variables and %d constraints from file %s\n",
	num_vars,num_cons,buf);
    last_col	= num_vars;
    last_row	= num_cons;

    for (i=1; i<= num_vars; ++i)
	{
	fscanf(luin," %s ",&varname[i][0]);
	strupr(&varname[i][0]);
	}

    for (i=0; i <= num_cons; ++i)
	{
	fscanf(luin," %s ",&rname[i][0]);		/* Row name */
	strupr(&rname[i][0]);
	for (k=1; k<=num_vars; ++k)
	    fscanf(luin,"%lf", &a_mat[i][k]);	/* Coefficients */

	if (!i) continue;			/* objective function row 0 */
	fscanf(luin," %s ",buf);			/* Constraint type */
	ch 		= buf[0];
	row_type[i] 	= (ch == '<')? 1 : ((ch == '>')? -1 : 0);
	fscanf(luin, " %lf", &rrhs[i]);
	a_mat[i][0]	= rrhs[i];
	}
    nnames = num_vars;
    }

setprob()
    {
	/*   Set up problem  */

	ph2_row	= 0;                              /*phase 2 row index*/
	last_col= num_vars;                      /*last column to pivot*/
	last_row= num_cons;                    /*last row to pivot*/
	ph1_row	= last_row + 1;

					/* set status of dec'n vars & OF*/
	for (col=0;col <= num_vars; ++col)
		{
		col_status[col]	= 0;
		cost[col] 	= a_mat[0][col];
		a_mat[0][col] 	*= z_mult;
		}

					/* get starting basis*/
	num_artif=0;
	for (row=1;row <= num_cons;++row)
	    {                               /*for each constraint,*/
	    if (row_type[row] > 0)  {
				/* a less-than constraint*/
				/* add a slack to the starting basis*/
				/* nnames should == last_col */
			if (addslk("SLK:",&rname[row][0]) < 0)
		    	{
		    	hit_max(2);
		    	exit(2);
		    	}
			++last_col;                             /*add new column*/
			a_mat[row][last_col]	= 1.0;          /*add slack*/
			upperb[last_col]	= infinity;     /*ub = infinity*/
			row_type[row]		= last_col;     /*set row pointer*/
			col_status[last_col]	= row;          /*set variable status*/
			basis_index[row]	= last_col;	/*set basis pointer*/
	   		}
	  else if (row_type[row] == 0) {

				/* an equality constraint*/
				/* add artificial to basis and tableau
				   for ranging reports, hide from basis */
			if (addslk("ART:",&rname[row][0]) < 0)
		    	{
		    	hit_max(2);
		    	exit(2);
		    	}
			++last_col;                             /*add new column*/
			isart[last_col]		= TRUE;		/*flag as artificial */
			a_mat[row][last_col]	= 1.0;          /*add art */
			upperb[last_col]	= infinity;     /*ub = infinity*/
			col_status[last_col]	= row;            /*set variable status*/
			basis_index[row] 	= last_col;    	/*artificial pointer*/
			++num_artif;
	  		}
	  else {
				/* a greater-than constraint*/
			/* add an artificial variable to the starting basis*/
			if (addslk("SRP:", &rname[row][0]) < 0)
		    	{
		    	hit_max(2);
		    	exit(2);
		    	}
			++last_col;                             /*add new column*/
			a_mat[row][last_col] 	= -1;		/*add surplus*/
			row_type[row]		= - last_col; 	/*set row pointer*/
			upperb[last_col] 	= infinity;	/*ub = infinity*/
			col_status[last_col] 	= 0;		/*set variable status*/

			if (addslk("ART:", &rname[row][0]) < 0)
		    	{
		    	hit_max(2);
		    	exit(2);
		    	}
			++last_col;                             /*add new artificial*/
			isart[last_col]		= TRUE;		/*flag as artificial */
			a_mat[ph1_row][last_col]= 1;		/*phase 1 obj fcn*/
			a_mat[row][last_col]	= 1.0;          /*add art */
			upperb[last_col]	= infinity;     /*ub = infinity*/
			col_status[last_col]	= row;            /*set variable status*/
			basis_index[row] 	= last_col;    	/*artificial pointer*/
			++num_artif;
	  		}
	}
				/* adjust rhs for lower bounds*/
	if (num_lbs > 0) adj_lbs();

	phase = (num_artif>0) ? 1 : 2;    /*set initial phase*/
				/* if phase 1 needed, set objective function*/
	if (phase == 2)  return;

	++last_row;                      /*add one more row to pivot*/
	min_cbar = 0;
	for (col=1,sum=0; col <= num_vars;sum=0,++col) {
	    for (row=1; row <= num_cons; ++row) {
		if (basis_index[basis_index[row]])   sum -= a_mat[row][col];
		}
		a_mat[ph1_row][col] = sum;
		if (sum <= min_cbar) {           /*check for best piv col*/
		    min_cbar = sum;
		    piv_col = col;
		}
	}

				/*  calculate phase 1 solution value*/
	z_1 = 0.0;
	for (row=1; row <= num_cons; ++row)
	    {
	    if (row_type[row] <= 0)  z_1 += a_mat[row][rhs_col]; 
	    }
	a_mat[ph1_row][rhs_col] = z_1;            /*set phase 1 value*/
}

simplex()
	/*    OPTIMIZATION DRIVER - RETURNS end_stat */
    {
    pivots=0;
    infeasible=0;
    unbounded = 0;

    if (phase == 1)
	{                                  /* perform phase 1*/
		if (show_exchanges) printf("BEGIN PHASE 1\n");
	obj_row = ph1_row;
	opt();                                  /*optimize*/
	if (end_stat == 1) infeas();            /*infeas check*/
	if (end_stat == 1) phase = 2;           /*skip ph 2?*/
	}

    if (phase == 2)
	{                                 /* perform phase 2*/
		if (show_exchanges) printf("BEGIN PHASE 2\n");
	obj_row = ph2_row;
	// last_row = num_cons;                /* Do not reset */
	opt();                                  /*optimize*/
	}

    if(last_tableau) tableau();
    }


opt()                   /*  o p t i m i z e    w. r. t.    o b j . r o w*/
    {
    end_stat = 0;
    if (debug) tableau();
    while (end_stat == 0) {
	if (debug) prt_stat();                             /*print current status*/
	incoming();                             /*pivot selection*/
	if (optimal)
	    {                                   /* stop if optimal */
	    end_stat=1;
	    break;
	    }
	outgoing();                             /*set incoming stat*/
	++pivots;
	ratio_test();                           /*ratio test*/
	if (unbounded)
	    {                       /*stop if unbounded*/
	    end_stat=3;
	    break;
	    }
	pivot();                                /*pivot*/
	if (debug) tableau();                   /*print tableau*/
	}
}


	/**********************************************************/
	/*  t e r m i n a t i n g    s t a t u s   r o u t i n e s*/
	/**********************************************************/

infeas()                    /* infeasibility check*/
    {
    for (row=1; row <= num_cons; ++row)
	{
	if (basis_index[row] == 0 && a_mat[row][rhs_col] > eps_infe)
	    {
	    infeasible = row;
	    end_stat=2;
	    return;
	    }
	}
    return;
    }

incoming()           /* Pivot selection:  pivot column */
    {
    min_cbar = -eps_price;
    best_col = 0;

    for (col=1; col <= last_col; ++col)
	{
	if (col_status[col] > 0 || isart[col]) continue; /* skip if basic/art */
	else if (col_status[col] == 0)
	    {                       /* currently at lower bound*/
	    if (a_mat[obj_row][col] <= min_cbar)
		{
		min_cbar = a_mat[obj_row][col];
		best_col = col;
		 }
	    }
	else
	    {                      /* currently at upper bound*/
	    if (-a_mat[obj_row][col] <= min_cbar)
		{
		min_cbar = -a_mat[obj_row][col];
		best_col = -col;
		}
	    }
	}
    if (best_col == 0)  optimal=1;
    else
	{
	optimal=0;
	piv_col=abs(best_col);
	}
    }

outgoing()            /*  Set entering variable status */
    {
    if (best_col < 0)  direction= -1.0;               /* now @ upper?? */
    else if (best_col > 0) direction= 1.0;            /* now @ lower?? */
    else err_exit(1);   /*error check*/
    }


ratio_test()
    {
    min_ratio = upperb[piv_col];
    best_row = 0;
					/* scan the pivot column*/
    for (row=1; row <= num_cons; ++row)
	{
	element = a_mat[row][piv_col] * direction;
	if (fabs(element)<eps_ratio) continue;     /*0 elem?*/
	if (element > 0)
	    {                      /* ratio 1:  basic to LB */
	    ratio = a_mat[row][rhs_col]/element;
	    ratio_row = row;
	 }
	else
	    {                     /* ratio 2:  basic to UB */
	    ratio=(upperb[basis_index[row]]-a_mat[row][rhs_col])/(-element);
	    ratio_row = -row;
	    }
			/* test for minimum ratio*/
	if (ratio < min_ratio ||
	   (ratio == min_ratio && basis_index[row]==artif))
	       {
	       min_ratio = ratio;
	       best_row = ratio_row;
	       }
	}
					/* >>> ONLY FOR UNBOUNDED CODE */
    if (best_row < 0)
	{
	unbounded = TRUE;
	end_stat  = 3;
	}

    piv_row = abs(best_row);
				      /* set leaving variable status */
    if (piv_row != 0)	/* <<< ERROR HERE, UNBOUNDED VAR can leave at UB */
	{                             /* min found?*/
	leav_var = basis_index[piv_row];
	col_status[leav_var] = (best_row > 0) ? 0 : -1;  /*to LB,UB*/
	}
		    /* unbounded or entering variable switches bounds*/
     else if (upperb[piv_col] == infinity)
	{
	unbounded = TRUE;                       /*problem unbounded*/
	end_stat=3;
	}
    }


pivot()
    {
    if (best_row != 0)
	{                     /*  normal basis exchange*/
	if (best_col<0) adjrhsin();             /*incoming from ub?*/
	if (show_exchanges)
		{
		int	ii;
		char	*outname = "ARTIFICIAL";
		if ((ii=basis_index[piv_row])!=0) outname = &varname[ii][0];
		printf("%s leaves and %s enters solution\n",
			outname, &varname[piv_col][0]);
		}
	piv_tabl();                             /*pivot tableau*/
	if (best_row<0) adjrhsout();            /*leaving at ub?*/
	basis_index[piv_row]=piv_col;           /*new basic variable*/
	col_status[piv_col]=piv_row;            /*new incoming status*/
	}

    else
	{                    /* bounds switch*/
	adjrhsin();                                /*adjust rhs for nb*/
	if      (col_status[piv_col] == 0)   col_status[piv_col]= -1;
	else if (col_status[piv_col] == -1)  col_status[piv_col]=0;
	else    err_exit(2);                       /*error in col_status*/
	}
    }


adjrhsin()                 /*      adjust rhs for incoming nonbasic*/
    {
    adj = upperb[piv_col]*direction;
    for (row=0; row <= last_row; ++row)
	a_mat[row][rhs_col] -= adj*a_mat[row][piv_col];
    }

adjrhsout()               /* adjust rhs for basic leaving at upper bound*/
    {
    out_col = basis_index[piv_row];
    adj     = upperb[out_col];
    for (row=0; row <= last_row; ++row)
	a_mat[row][rhs_col] -= adj*a_mat[row][out_col];
    }

piv_tabl()                            /* pivot tableau */
    {
    piv_elem = a_mat[piv_row][piv_col];

				/* update pivot row*/
    for (col=0; col <= last_col;++col)
	    a_mat[piv_row][col]=a_mat[piv_row][col]/piv_elem;

				/* update non-pivot rows*/
    for (row=0; row <= last_row; ++row)
	{
	if (row == piv_row)  continue;           /*skip pivot row*/
	elt = a_mat[row][piv_col];
	for (col=0; col <= last_col; ++ col)
	    {
	    if (col != piv_col)
		{
		a_mat[row][col] -= elt*a_mat[piv_row][col];
		if (fabs(a_mat[row][col]) < eps_pivot)  a_mat[row][col]=0.0;
		}
	    }
	}

				/* zero out pivot column*/
    for (row=0; row <= last_row; ++row) a_mat[row][piv_col]=0.0;
    a_mat[piv_row][piv_col] = 1.0;
    }


prt_stat()                      /* print current status*/
    {
    printf("Iteration %d  Objective = %11.3lf \n",pivots,-z_mult*a_mat[0][0]);
    }

out_opt()
    {
    var_summ();                     /* print variable summary */
    if(show_range) ranger();
    }

out_infe()
     /*      write infeasible problem */
   {
   printf("PROBLEM INFEASIBLE\n");
   infe_summ();
   }

infe_summ()
    {
    register int row;

    strcpy(titlevar," INFEASIBILITY REPORT");
    nrpt = 0;

    printf("The problem is infeasible -- contradictory constraints have been\n");
    printf("posed, such as:   X<10 and X>20.     Specifically, the following\n");
    printf("constraints could not be satisfied:\n\n");
    printf("Constraint        RHS Violation\n");
    printf("----------        -------------\n");
    for(row=1; row <= num_cons; ++row)
	{
	if (basis_index[row]==0 && a_mat[row][rhs_col]>eps_infe)
	    {
	    printf("%2d. %8s       %12.3lf\n",
		row, &rname[row][0],rrhs[row]);
	    }
	}

    titlecons[0]=EOS;
}


out_unbnd()                      /* write unbounded output file*/
    {
    printf("UNBOUNDED PROBLEM\n");
    printf("Variable %s can be increased without limit.\n",
	&varname[piv_col][0]);
    printf("Please reformulate the model.\n");
    nrpt = -1;

    }


string(buf,len,ch)      /* Sets buffer to string of 'len' chars */
char buf[],ch;
int len;
    	{
	setmem(buf,len,ch);
   	buf[len] = '\0';
    	}

addslk(prefix, suffix)
char *prefix, *suffix;
    {
    extern char buf[];

    strcpy(buf, prefix);
    strcat(buf, suffix);
    buf[SIZE_LONG_NAME] = EOS;
    return nameNum(buf);
    }

chset()
    {
    extern char chtype[];
    register int i;

    for (i=0; i < 256;)   chtype[i++] = -2;         /* Default = error */
    for (i='A'; i <= 'Z';) chtype[i++] = 1;         /* 1. Upper case alpha */
    for (i='a'; i <= 'z';) chtype[i++] = 1;         /* 1. Lower case alpha */
    for (i='0'; i <= '9';) chtype[i++] = 2;          /* 2. Digit */
    chtype['.']     = 3;                            /* 3. Period */
    chtype['+']     = 4;                            /* 4. Signs */
    chtype['-']     = 4;
    chtype[EOS]     = 5;                            /* 5. End chars */
    chtype['\n']    = 5;
    chtype[BLANK]   = 5;
    chtype['\t']    = 5;
    chtype['<']     = 6;                             /* 6. Relation */
    chtype['=']     = 6;
    chtype['>']     = 6;
    chtype[')']	    = 7;
    }


tableau()                             /* Print out Tableau */
    {
    int col1,col2;
    char *dashes;
    FILE *lu;
    dashes="---------------------------------------";
	static char* dashes7  = "-------";
	static char* dashes14 = "--------------";
	static int printcols = 50;

    lu = stdout;
    for (col1=0; col1<=last_col;col1 += printcols)
		{
		col2 = (last_col-col1 < printcols-1) ? last_col : col1+printcols;

		fprintf(lu,"\n  Name ");
		for (col=col1; col <= col2; ++col)
	    	fprintf(lu,"  %12s", &varname[col][0]);

		fprintf(lu,"\n   Cost");
		for (col=col1; col <= col2; ++col)
			fprintf(lu,"  %12.3lf",cost[col]);
		fprintf(lu,"\n Status");
		for (col=col1; col <= col2; ++col)
			fprintf(lu," %13d",col_status[col]);

		//fprintf(lu,"\n%s%s\n",dashes,dashes);
		fprintf(lu,"\n%s",dashes7);
		for (col=col1;col <= col2;++col) fprintf(lu,"%s",dashes14);
		fprintf(lu,"\n");

		fprintf(lu,"Row Bas");
		for (col=col1;col <= col2;++col)
		    fprintf(lu,"     Column%3d",col);
		//fprintf(lu,"\n%s%s\n",dashes,dashes);
		fprintf(lu,"\n%s",dashes7);
		for (col=col1;col <= col2;++col) fprintf(lu,"%s",dashes14);
		fprintf(lu,"\n");

		for (row=0;row <= last_row;++row)
	    	{
	    	fprintf(lu,"%2d  %3d",row,basis_index[row]);
	    	for (col=col1;col <= col2;++col)
		    	fprintf(lu,"  %12.3lf",a_mat[row][col]);
	    	fprintf(lu,"\n");
	    	}
		fprintf(lu,"\n\n");
		}
    }

adj_lbs() { /* Dummy function */ }

hit_max(i)
int i;
    {
    printf("MAX DIMENSIONS EXCEEDED  Error code %d\n",i);
    }

var_summ()                   /* Decision Variable and Constraint Summary  */
    {
    char *type,*slk_id;
    int c_status;
    double marginal,slack;

    nrpt = 1;

    printf("Solution OPTIMAL in %d iterations:    %3.3s %s = %lg\n",
	 pivots, minmax, &rname[0][0], -z_mult*a_mat[0][0]);
    ctitle("*** DECISION VARIABLE SUMMARY ***",60);

    sprintf(buf,"Unit %s",&rname[0][0]);
    printf("Variable            Opt Value   %13.13s    Marg Value ",buf);
    if (bounded) printf("Upper Bound   Lower Bound\n");
    else printf("\n");
    printf("------------      -----------    ------------    ---------- ");
    if (bounded) printf("-----------   -----------\n");
    else printf("\n");

    for (col=1; col <= num_vars; ++col)
	{
	printf("%2d. %-8.8s",col,&varname[col][0]);

	if (col_status[col] < 0) x=upperb[col]+lowerb[col];
	else if (col_status[col] == 0) x=lowerb[col];
	else x=a_mat[col_status[col]][rhs_col]+lowerb[col];

	if (x >= eps_infe) printf(" %16.4lf %15.3lf",x,cost[col]);
	else  printf("            .     %15.3lf",cost[col]);

	if (fabs(a_mat[0][col]) < eps_price) printf("          .   ");
	else printf(" %13.3lf",a_mat[0][col]*z_mult);

	if (!bounded) { printf("\n"); continue; }

	if (upperb[col] == infinity) printf("    +Infinity");
	else printf("%13.3lf",upperb[col]+lowerb[col]);

	if (lowerb[col] > 0.0) printf(				"%13.3lf\n",lowerb[col]);
	else  printf("        .   \n");
	}

	ctitle("*** CONSTRAINT SUMMARY ***",72);
	printf(
	"Constraint   Rel         RHS           Slack/Surplus    Marginal Value\n");
	printf(
	"------------  --    ------------    ------------------  --------------\n");

	for (row=1; row <= num_cons; ++row)
	    {
	    s_col=row_type[row];
	    if (s_col > 0)
		{                              /* LE constraint */
		type="<=";
		slk_id="SLK";
		c_status=col_status[s_col];
		slack = (c_status == 0) ? 0.0 : a_mat[c_status][rhs_col];
		marginal=a_mat[0][s_col]*(-z_mult);
		}
	    else if (s_col == 0)
		{                             /* EQ constraint */
		type=" =";
		marginal=a_mat[0][row+num_vars]*(-z_mult);
		if (basis_index[row] != artif)
		    {
		    slk_id="   ";
		    slack=0.0;
		    }
		else
		    {
		    slk_id="ART";
		    slack=a_mat[row][rhs_col];
		    }
		}
	    else
		{                            /* GE constraint */
		s_col= -s_col;
		type = ">=";
		slk_id="SRP";
		marginal=a_mat[0][s_col]*z_mult;
		c_status=col_status[s_col];
		slack = (c_status == 0)? 0.0 : a_mat[c_status][rhs_col];
		}

	    printf("%2d. %-8.8s  %2.2s %15.3lf",
		row, &rname[row][0], type, rrhs[row]);

	    if (slack == 0.0)  printf("            .         ");
	    else printf(" %15.3lf (%.3s)",slack,slk_id);

	    if (fabs(marginal) < .0001) printf("            .\n");
	    else printf(" %15.3lf\n",marginal);
	}
    }

err_exit(i)
int i;
    {
    printf(">>>ERROR EXIT %d <<<\n",i);
	exit(-1);
    }

ranger()
{
	range_cost();
	range_rhs();
}

range_rhs()
/* ****************************************************************** */
/*             ranging the right-hand-side */
/* ****************************************************************** */
{
	int min_row,max_row;
	register int cons, slk_col,i ;
	double 	direction, min_ratio, max_ratio;
	char	*fd1 = "%13.2lf %-12.12s ";
	char	*fd2 = "%13.2lf %-12.12s\n";
	char	*finf = "    +Infinity              ";
	char	*fminf= "    -Infinity\n";

	ctitle("*** RANGING THE RIGHT-HAND SIDE ***",80);
	printf("Constraint         RHS   ");
	printf("   Upper Limit  (Becomes 0)");
	printf("   Lower Limit  (Becomes 0)\n");
	printf("--------------------------  ------------------------");
	printf("   ------------------------\n");

	slk_col	= num_vars+1;
	for (cons=1; cons <= num_cons; ++cons, ++slk_col )
		{
		printf("%2d. %-8s %13.2lf",cons,&rname[cons][0],rrhs[cons]);

					/* Case A:  Slk/Srp basic? */
		if ((i=col_status[slk_col]) > 0)
			{
			if (row_type[cons]>0)		/* LE constraint? */
				{
				printf(finf);
				printf(fd2,rrhs[cons]-a_mat[i][0],&varname[slk_col][0]);
				}
			else	{			/* GE constraint */
				printf(fd1,rrhs[cons]+a_mat[i][0],&varname[slk_col][0]);
				printf(fminf);
				}
			continue;
			}

					/* Case B:  Slk/Srp nonbasic */
		direction	= (row_type[cons]>=0) ? 1.0 : -1.0;
		min_ratio	= infinity;
		max_ratio	= -infinity;
		for (row=1; row <= num_cons; ++row)
			{
			if (a_mat[row][slk_col] == 0.0) continue;
			ratio	= direction*a_mat[row][rhs_col] / a_mat[row][slk_col];
						/* min of + ratios */
			if (ratio > 0.0 ||
			   (ratio == 0.0 && a_mat[row][slk_col] > 0.0))
				{
				if (ratio <= min_ratio)
					{
					min_ratio 	= ratio;
					min_row		= row;
					}
				}
						/* max of - ratios */
			else 	{
				if (ratio >= max_ratio)
					{
					max_ratio	= ratio;
					max_row		= row;
					}
				}
			}

		if (max_ratio > -infinity) {
			printf(fd1,rrhs[cons]-max_ratio,
				&varname[basis_index[max_row]][0]);
		}
		else printf(finf);

		if (min_ratio < infinity) {
			printf(fd2,rrhs[cons]-min_ratio,
				&varname[basis_index[min_row]][0]);
		}
		else printf(fminf);
	}
}

range_cost()
{
	int var;

	sprintf(buf,"*** SENSITIVITY OF THE %s VALUES ***",	&rname[0][0]);
	ctitle(buf,80);
	printf("Variable        Current OF    Upper Limit (Becomes+)");
	printf("    Lower Limit (Becomes+)\n");
	printf("--------------------------    ----------------------");
	printf("    -----------------------\n");

	for (var=1; var <= num_vars; ++var) {
		printf("%2d. %-8.8s %13.3lf ",var,varname[var],cost[var]);
		if (col_status[var] > 0) range_cb(var);
		else range_cnb(var);
	}
}

range_cb(var)
int var;
/*
	Range cost on basic variable
*/
{
    double element,max_neg_ratio,min_pos_ratio,sign;
    int neg_col,pos_col;

    max_neg_ratio = -infinity;
    min_pos_ratio =  infinity;
	for (col=1;col <= last_col; ++col) {
		if (col_status[col] > 0 || isart[col] ||
		    fabs((element=a_mat[col_status[var]][col])) <= eps_pivot)
		    continue;

		sign = (col_status[col] < 0)? -1.0 : 1.0;
        ratio=a_mat[0][col]/element*z_mult;
		if (sign*z_mult*element >= 0) {
            if (ratio > 0.0 && ratio < min_pos_ratio) {
                min_pos_ratio=ratio;
                pos_col=col;
			}
		}
		else {
            if (ratio < 0.0 && ratio > max_neg_ratio)  {
                max_neg_ratio=ratio;
                neg_col=col;
			}
		}
	}

    if (min_pos_ratio < infinity) {
        printf("%13.3lf ",(cost[var]+min_pos_ratio));
        printf("%-12.12s",&varname[pos_col][0]);
	}
	else printf("    +Infinity             ");

    if (max_neg_ratio > -infinity) {
        printf("%13.3lf ",(cost[var]+max_neg_ratio));
        printf("%-12.12s\n",&varname[neg_col][0]);
	}
	else printf("    -Infinity             \n");
}

range_cnb(var)
int var;
/*
	Range cost on Nonbasic variable
*/
{
	double sign;

	sign = (col_status[var] == 0)? 1.0 : -1.0;
	if (z_mult*sign < 0) {
		printf("%13.3lf ",cost[var]-z_mult*a_mat[0][var]);
		printf("%-12.12s",&varname[var][0]);
		printf("    -Infinity\n");
	}
	else {
		printf("    +Infinity             ");
		printf("%13.3lf ",cost[var]-z_mult*a_mat[0][var]);
		printf("%-12.12s\n",&varname[var][0]);
	}
}

ctitle(title,width)
char *title; int width;
	{
	register int i, indent;
	printf("\n");
	indent	= (width - strlen(title))/2;
	for(i=1; i<= indent; ++i) printf(" ");
	printf("%s\n\n",title);
	}

help()
	{
	char	*msg =
"           B L P -- Linear Programming, (C) 2001, Richard S. Barr.\n\
Program usage: 		blp filename\n\
where 'filename' is a file containing the following LP model description:\n\
   1.  'MIN' or 'MAX' for type of optimization.\n\
   2.  A row label for the objective function.\n\
   3.  A linear expression for the objective function\n\
   4.  Up to %d constraints of the form:\n\
	  Row label   Linear expression    Relation   Non-negative constant\n\
where: Label or variable name p 1-%d letters and digits, beginning with letter\n\
       Label p a name, followed by a ')'  [e.g.  ROW1) ]\n\
       Linear expression p linear sum of variables multiplied by constants\n\
	  [e.g., 2X - 13.5ALPHA+ZETA].  The model may use up to %d variables.\n\
       Relation p the characters '<=', '>=', or '='\n\
Model description is free-form.  Names may not contain blanks.\n\
\n\
User options can be invoked using the form:	blp -options filename\n\
where 'options' is a string of one or more of the following characters\n\
	r = display Ranging reports\n\
	f = display Final simplex tableau\n\
	e = do not Echo input data\n\
	i = Input model in matrix format, rather than above algebraic format\n\
	m = do not display LP Model\n\
	x = display variables eXchanged at each iteration\n";

	printf(msg,MAX_ROWS,SIZE_NAME,MAX_VARS);
	}


#define LINELEN 79
outTerm(s,flag)
char *s; int flag;
	{
	register int sLen;
    static char outBuf[LINELEN+2], outLen, *indent = "        ";

	switch (flag)
		{
		case 1:			/* initialize and process s */
			outBuf[0]	= EOS;
			outLen		= 0;
		case 0:			/* Process string s */
			sLen		= strlen(s);
			if(sLen+outLen >= LINELEN)
				{
				printf("%s\n",outBuf);
				strcpy(outBuf,indent);
				}
			strcat(outBuf,s);
            outLen  = strlen(outBuf);
            break;
		case -1:		/* Flush */
			printf("%s\n",outBuf);
			outBuf[0]	= EOS;
			outLen		= 0;
		}
	}

makeTerm(val,first,name,s)		/* Convert nonzero term to string */
double val; int first; char *name, *s;
	{
	double 		absVal;

	absVal	= fabs(val);
	if (fabs(val-1.0)<eps_print)
		if (first) sprintf(s," %s",name);
		else sprintf(s," + %s",name);
	else if (fabs(val+1.0) < eps_print)
		sprintf(s," - %s",name);
	else 	if (first)
			if (val>=0.0) sprintf(s," %g %s", (float)absVal, name);
			else        sprintf(s," - %g %s", (float)absVal, name);
		else
			if (val>=0.0) sprintf(s," + %g %s", (float)absVal, name);
			else	    sprintf(s," - %g %s", (float)absVal, name);
	}

showModel()
	{
	register int 	i,j, firstTerm;
	char	 	s[80];

	ctitle("*** LP MODEL ***",70);
	if (z_mult < 1.0) outTerm("Maximize ",1);
	else	outTerm("Minimize ",1);
	outLhs(0,&cost[0],s);
	outTerm(s,-1);

	outTerm("subject to:",1);
	outTerm(s,-1);

	for (i=1; i<=num_cons; ++i)
		{
		outLhs(i,&a_mat[i][0],s);
		outRhs(i,s);
		}
	printf("\n");
	}

outLhs(i,rowVal,s)
int i; double rowVal[]; char *s;
	{
	register int j, firstTerm;
	double	absVal, val;
	char temp[15];

	sprintf(temp,"%s)",&rname[i][0]);
	sprintf(s,"%-9.9s",temp);
	outTerm(s,0);
	for (j=1, firstTerm=TRUE; j<=num_vars; ++j)
		{
		absVal	= fabs( val = rowVal[j] );
		if (absVal > eps_print)
			{
			makeTerm(val, firstTerm, &varname[j][0], s);
			outTerm(s,0);
			firstTerm = FALSE;
			}
		}
	}

outRhs(i,s)
int i; char *s;
	{
	char rel;

	if (row_type[i]<0) 	rel = '>';
	else if (row_type[i]>0) rel = '<';
	else 			rel = ' ';
	sprintf(s," %c= %g",rel, (float) a_mat[i][0]);
	outTerm(s,0);
	outTerm(s,-1);
	}