#include <stdio.h>
#include "opt.h"

#define	Debug	0

/*
	Simplex Method (linear programming)
	solve the canonical form of LP problem
	
	maximize
	    c0 x0 +     c1 x1 + ... +     cn-m+1 xn-m+1 + d0
	
	subject to
	
	  a0,0 x0 +   a0,1 x1 + ... +   a0,n-m+1 xn-m+1 + xn-m		= b0
	  a1,0 x0 +   a1,1 x1 + ... +   a1,n-m+1 xn-m+1     + xn-m+1	= b1
	   		......
	am-1,0 x0 + am-1,1 x1 + ... + am-1,n-m+1 xn-m+1          + xn-1	= bm-1
	
	inputs: m --- the number of equations
		n --- the number of variables
			m < n
		a --- coefficient matrix
			note that matrix a satisfies
				a[i][n-m+j] = 1 (when i=j)
				          = 0 (when i \neq j)
				for all i=0,1,...,m-1
					j=0,1,...,m-1
		b --- constant vector
			note that b[i] >= 0 for all i=0,1,...,m-1
		c --- coefficient vector
			note that c[n-m+j] = 0 for all j=0,1,...,m-1
		d --- constant
		eps --- small value

	outputs: ox --- optimal solution
		 ov --- optimal value

	return value:	success : find optimal solution
			failure : diverge (optimal value = infinity)
*/

extern status LPsimplex
(int m, int n, dbl a[], dbl b[], dbl c[], dbl *d, dbl eps, dbl ox[], dbl *ov)
{
	int	i, j;
	int	*base;
	status	st;
	
	base = allc(dbl, m);
	for (i=0; i<m; i++) {
		base[i] = n-m+i;
	}
	
	st = LPsimplextableau(m, n, a, b, c, d, eps, base);
	if (st == success) {
		LPsimplexsolution(m, n, b, d, base, ox, ov);
	}
	
	return st;
}

extern status LPsimplextableau
(int m, int n, dbl a[], dbl b[], dbl c[], dbl *d, dbl eps, int base[])
{
	int	i, j, r, s;
	dbl	v, min, mul;
	
	for (;;) {
		/*	column pivot	*/
		/* find c[j] > 0 */
		s = n;
		for (j=0; j<n; j++) {
			if (c[j] > eps) {
				s = j;
				break;
			}
		}
		if (s == n) {
			/* when c0 through cn-1 are all non-positive */
			return(success);
		}
#if Debug
		printf("col pivot = %d\n", s);
#endif
		
		/*	row pivot	*/
		r = m;
		for (i=0; i<m; i++) {
			if (a[n*i+s] <= eps) continue;
			v = b[i]/a[n*i+s];
			if (r == m || v < min) {
				r = i;
				min = v;
			}
		}
		if (r == m) {
			/* when objective func. diverges */
			return(failure);
		}
#if Debug
		printf("row pivot = %d\n", r);
#endif

		/*	pivot operation		*/
		base[r] = s;
		
		mul = 1/a[n*r+s];
		matrixrowscalar(m, n, a, r, mul);
		matrixrowscalar(m, 1, b, r, mul);
		
		for (i=0; i<m; i++) {
			if (i == r) continue;
			mul = - a[n*i+s];
			matrixrowadd(m, n, a, i, mul, r);
			matrixrowadd(m, 1, b, i, mul, r);
		}
		
		mul = c[s];
		for (j=0; j<n; j++) {
			c[j] -= mul * a[n*r+j];
		}
		*d += mul * b[r];
#if Debug
		matrixfprint(stdout, m, n, a);
		matrixfprint(stdout, m, 1, b);
		matrixfprint(stdout, 1, n, c);
#endif
	}
}

extern status TwoPhasemethodtableau
(int m, int n, dbl a[], dbl b[], dbl c[], dbl *d,
 dbl cadj[], dbl *dadj, dbl eps, int base[])
{
	int	i, j, r, s;
	dbl	v, min, mul;
	
	for (;;) {
		/*	column pivot	*/
		/* find c[j] > 0 */
		s = n;
		for (j=0; j<n; j++) {
			if (c[j] > eps) {
				s = j;
				break;
			}
		}
		if (s == n) {
			/* when c0 through cn-1 are all non-positive */
			return(success);
		}
#if Debug
		printf("col pivot = %d\n", s);
#endif
		
		/*	row pivot	*/
		r = m;
		for (i=0; i<m; i++) {
			if (a[n*i+s] <= eps) continue;
			v = b[i]/a[n*i+s];
			if (r == m || v < min) {
				r = i;
				min = v;
			}
		}
		if (r == m) {
			/* when objective func. diverges */
			return(failure);
		}
#if Debug
		printf("row pivot = %d\n", r);
#endif

		/*	pivot operation		*/
		base[r] = s;
		
		mul = 1/a[n*r+s];
		matrixrowscalar(m, n, a, r, mul);
		matrixrowscalar(m, 1, b, r, mul);
		
		for (i=0; i<m; i++) {
			if (i == r) continue;
			mul = - a[n*i+s];
			matrixrowadd(m, n, a, i, mul, r);
			matrixrowadd(m, 1, b, i, mul, r);
		}
		
		mul = c[s];
		for (j=0; j<n; j++) {
			c[j] -= mul * a[n*r+j];
		}
		*d += mul * b[r];
		
		mul = cadj[s];
		for (j=0; j<n; j++) {
			cadj[j] -= mul * a[n*r+j];
		}
		*dadj += mul * b[r];
#if Debug
		matrixfprint(stdout, m, n, a);
		matrixfprint(stdout, m, 1, b);
		matrixfprint(stdout, 1, n, c);
#endif
	}
}

extern void LPsimplexsolution(int m, int n, dbl b[], dbl *d, int base[],
			      dbl ox[], dbl *ov)
{
	int	i, j;
	
	for (j=0; j<n; j++) {
		ox[j] = 0.00;
	}
	for (i=0; i<m; i++) {
		ox[base[i]] = b[i];
	}
	*ov = *d;
}