From jrs296@psu.edu Thu Apr 20 20:20 EDT 2000
X-Mailer: Microsoft Outlook Express Macintosh Edition - 4.5 (0410)
Date: Thu, 20 Apr 2000 20:20:20 -0400
Subject: Code for performing Simplex Method
From: "Josh Sweetwood" <jrs296@psu.edu>
To: vstein@math.psu.edu
Mime-version: 1.0
X-Priority: 3
Content-Type: multipart/mixed; boundary="MS_Mac_OE_3039106821_258926_MIME_Part"
Content-Length: 15845

Dr. Vaserstein,

This is the C++ code for performing the Simplex Method.  It is
relatively simple to use.  You just enter a matrix in the standard tableau
for we use in class and it shows you the steps it goes through to get the
answer (max).  I can talk about it in class on Tuesday.  If you have a
problem with the format of the attachments just let me know.  Feel free to
e-mail me with any questions.  Thanks.

Josh Sweetwoo
---------
// Fraction Class by Josh Sweetwood

#include <iostream.h>

class Fraction
{
private:
int Numer;
int Denom;

const Fraction Reduce( );

public:
// Constructors
Fraction( int Num = 0, int Den = 1 ) :
Numer( Num ), Denom( Den ) { }
Fraction( const Fraction & Value );

// Destructor
~Fraction( ) { }

// Assignments
const Fraction & operator = ( const Fraction & Value );
const Fraction & operator += ( const Fraction & Value );
const Fraction & operator -= ( const Fraction & Value );
const Fraction & operator *= ( const Fraction & Value );
const Fraction & operator /= ( const Fraction & Value );

// Accessors
int Get_Numer( ) const { return Numer; }
int Get_Denom( ) const { return Denom; }
double Get_Decimal( )
{
double Decimal;
Decimal = 0;
if ( Denom != 0 )
Decimal = ( double(Numer) / double(Denom) );
return Decimal;
}

// Nonmember friends
private:
friend ostream & operator <<
( ostream & Out, const Fraction & Value );
friend Fraction operator *
( const Fraction & A, const Fraction & B );
friend Fraction operator /
( const Fraction & A, const Fraction & B );
friend Fraction operator +
( const Fraction & A, const Fraction & B );
friend Fraction operator -
( const Fraction & A, const Fraction & B );
friend Fraction operator -
( const Fraction & A );
friend int Gcd ( int M, int N );
friend int operator <
( const Fraction & A, const Fraction & B );
friend int operator >
( const Fraction & A, const Fraction & B );
friend int operator <=
( const Fraction & A, const Fraction & B );
friend int operator >=
( const Fraction & A, const Fraction & B );
friend int operator ==
( const Fraction & A, const Fraction & B );
friend int operator !=
( const Fraction & A, const Fraction & B );
};

//------------------------------------------------------//
// Implementation Starts Here                           //
//------------------------------------------------------//

Fraction::Fraction( const Fraction & Value )
{
Numer = Value.Numer;
Denom = Value.Denom;
}

int Gcd ( int M, int N )
{
int Rem;

if (M < 0 )
M *= -1;
if (N < 0 )
N *= -1;

while ( N > 0 )
{
Rem = M % N;
M = N;
N = Rem;
}
return M;
}

const Fraction Fraction::Reduce ( )
{
int ThisGcd;

ThisGcd = Gcd( Numer, Denom );
Numer /= ThisGcd;
Denom /= ThisGcd;

if ( (Numer < 0 && Denom < 0) || (Numer >= 0 && Denom < 0) )
{
Numer *= -1;
Denom *= -1;
}

return *this;
}

const Fraction & Fraction::operator = ( const Fraction & Value )
{
if ( this != &Value )
{
Numer = Value.Numer;
Denom = Value.Denom;
Reduce( );
}
return *this;
}

const Fraction & Fraction::operator += ( const Fraction & Value )
{
int LHNumer, RHNumer, LHDenom, RHDenom;
Fraction RValue;

Reduce( );
RValue = Value;
RValue.Reduce( );

LHNumer = Numer * RValue.Denom;
LHDenom = Denom * RValue.Denom;
RHNumer = RValue.Numer * Denom;
RHDenom = RValue.Denom * Denom;

Numer = LHNumer + RHNumer;
Denom = LHDenom;

Reduce( );

return *this;
}

const Fraction & Fraction::operator -= ( const Fraction & Value )
{
int LHNumer, RHNumer, LHDenom, RHDenom;
Fraction RValue;

Reduce( );
RValue = Value;
RValue.Reduce( );

LHNumer = Numer * RValue.Denom;
LHDenom = Denom * RValue.Denom;
RHNumer = RValue.Numer * Denom;
RHDenom = RValue.Denom * Denom;

Numer = LHNumer - RHNumer;
Denom = LHDenom;

Reduce( );

return *this;
}

const Fraction & Fraction::operator *= ( const Fraction & Value )
{
Fraction RValue;

Reduce( );
RValue = Value;
RValue.Reduce( );

Numer *= RValue.Numer;
Denom *= RValue.Denom;

Reduce( );

return *this;
}

const Fraction & Fraction::operator /= ( const Fraction & Value )
{
Fraction RValue;

Reduce( );
RValue = Value;
RValue.Reduce( );

Numer *= RValue.Denom;
Denom *= RValue.Numer;

Reduce( );

return *this;
}

ostream & operator << ( ostream & Out, const Fraction & Value )
{
Fraction RValue;

RValue = Value;
RValue.Reduce( );

if ( RValue.Denom == 0 )
Out << "INF";
else if ( RValue.Denom == 1 || RValue.Numer == 0 )
Out << RValue.Numer;
else
Out << RValue.Numer << "/" << RValue.Denom;

return Out;
}

Fraction operator + ( const Fraction & A, const Fraction & B )
{
int LHNumer, RHNumer, LHDenom, RHDenom;

RA = A;
RB = B;

RA.Reduce( );
RB.Reduce( );

LHNumer = RA.Numer * RB.Denom;
LHDenom = RA.Denom * RB.Denom;
RHNumer = RB.Numer * RA.Denom;
RHDenom = RB.Denom * RA.Denom;

}

Fraction operator - ( const Fraction & A, const Fraction & B )
{
int LHNumer, RHNumer, LHDenom, RHDenom;

RA = A;
RB = B;

RA.Reduce( );
RB.Reduce( );

LHNumer = RA.Numer * RB.Denom;
LHDenom = RA.Denom * RB.Denom;
RHNumer = RB.Numer * RA.Denom;
RHDenom = RB.Denom * RA.Denom;

}

Fraction operator - ( const Fraction & A )
{
Fraction RA;

RA = A;

RA.Reduce( );

return Fraction ( -RA.Numer, RA.Denom );
}

Fraction operator * ( const Fraction & A, const Fraction & B )
{

RA = A;
RB = B;

RA.Reduce( );
RB.Reduce( );

Answer = Fraction ( RA.Numer * RB.Numer, RA.Denom * RB.Denom );

}

Fraction operator / ( const Fraction & A, const Fraction & B )
{

RA = A;
RB = B;

RA.Reduce( );
RB.Reduce( );

Answer = Fraction ( RA.Numer * RB.Denom, RA.Denom * RB.Numer );

}

int operator < ( const Fraction & A, const Fraction & B )
{
int LHCross, RHCross;
Fraction RA, RB;

RA = A;
RB = B;

RA.Reduce( );
RB.Reduce( );

LHCross = B.Denom * A.Numer;
RHCross = A.Denom * B.Numer;

if ( LHCross < RHCross )
return 1;
else
return 0;
}

int operator > ( const Fraction & A, const Fraction & B )
{
int LHCross, RHCross;
Fraction RA, RB;

RA = A;
RB = B;

RA.Reduce( );
RB.Reduce( );

LHCross = B.Denom * A.Numer;
RHCross = A.Denom * B.Numer;

if ( LHCross > RHCross )
return 1;
else
return 0;
}

int operator <= ( const Fraction & A, const Fraction & B )
{
int LHCross, RHCross;
Fraction RA, RB;

RA = A;
RB = B;

RA.Reduce( );
RB.Reduce( );

LHCross = B.Denom * A.Numer;
RHCross = A.Denom * B.Numer;

if ( LHCross <= RHCross )
return 1;
else
return 0;
}

int operator >= ( const Fraction & A, const Fraction & B )
{
int LHCross, RHCross;
Fraction RA, RB;

RA = A;
RB = B;

RA.Reduce( );
RB.Reduce( );

LHCross = B.Denom * A.Numer;
RHCross = A.Denom * B.Numer;

if ( LHCross >= RHCross )
return 1;
else
return 0;
}

int operator == ( const Fraction & A, const Fraction & B )
{
Fraction RA, RB;

RA = A;
RB = B;

RA.Reduce( );
RB.Reduce( );

if ( RA.Numer == RB.Numer && RA.Denom == RB.Denom )
return 1;
else
return 0;
}

int operator != ( const Fraction & A, const Fraction & B )
{
Fraction RA, RB;

RA = A;
RB = B;

RA.Reduce( );
RB.Reduce( );

if ( RA.Numer != RB.Numer || RA.Denom != RB.Denom )
return 1;
else
return 0;
}
------------
#include "Fraction.h"

const int Maxcols = 10;
const int Maxrows = 10;

Fraction Matrix[Maxrows][Maxcols];

// Function Prototypes because my compiler sucks
void Enter_Dimensions(int & x, int & y);
void Print_Matrix(int rows, int cols);
void Enter_Matrix(int rows, int cols);
void Pivot(int rows, int cols, int prow, int pcol);
int Feasible(int rows, int cols);
int Optimal(int rows, int cols);
void Phase1_Pick_Pivot(int rows, int cols, int & prow, int & pcol);
void Phase2_Pick_Pivot(int rows, int cols, int & prow, int & pcol);

void Enter_Dimensions(int & x, int & y)
{
while ( x <= 0 || x > Maxrows )
{
cout << "Enter the number of    Rows (max = " << Maxrows << "): ";
cin >> x;
}
while ( y <= 0 || y > Maxcols )
{
cout << "Enter the number of Columns (max = " << Maxcols << "): ";
cin >> y;
}
}

void Print_Matrix(int rows, int cols)
{
for (int i = 0; i < rows; i++)
{
cout << "{";
for (int j = 0; j < cols; j++)
{
cout << Matrix[i][j];
if ( (j + 1) < cols )
cout << ",";
}  cout << "}" << endl;
}
}

void Enter_Matrix(int rows, int cols)
{
int entry;

for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
{
cout << "Enter integer for row " << (i+1)
<< ", column " << (j+1) << ": ";
cin >> entry;
Matrix[i][j] = entry;
}
}
}

void Pivot(int rows, int cols, int prow, int pcol)
{
Fraction Orig_Matrix[Maxrows][Maxcols];

// Save a copy of Matrix before pivoting
for (int i = 0; i < rows; i++)
{
for (int j = 0; j < cols; j++)
{
Orig_Matrix[i][j] = Matrix[i][j];
}
}

// Calculate "A"
Matrix[prow][pcol] = 1 / Orig_Matrix[prow][pcol];

// Calculate "B"
for (int j = 0; j < pcol; j++)
Matrix[prow][j] = Orig_Matrix[prow][j] / Orig_Matrix[prow][pcol];

for (int j = (pcol+1); j < cols; j++)
Matrix[prow][j] = Orig_Matrix[prow][j] / Orig_Matrix[prow][pcol];

// Calculate "C"
for (int i = 0; i < prow; i++)
Matrix[i][pcol] = -(Orig_Matrix[i][pcol] / Orig_Matrix[prow][pcol]);

for (int i = (prow+1); i < rows; i++)
Matrix[i][pcol] = -(Orig_Matrix[i][pcol] / Orig_Matrix[prow][pcol]);

// Calculate "D"
for (int i = 0; i < prow; i++)
{
for (int j = 0; j < pcol; j++)
{
Matrix[i][j] = (Orig_Matrix[i][j] -
((Orig_Matrix[prow][j] * Orig_Matrix[i][pcol]) /
Orig_Matrix[prow][pcol]));
}
}

for (int i = (prow+1); i < rows; i++)
{
for (int j = 0; j < pcol; j++)
{
Matrix[i][j] = (Orig_Matrix[i][j] -
((Orig_Matrix[prow][j] * Orig_Matrix[i][pcol]) /
Orig_Matrix[prow][pcol]));
}
}

for (int i = 0; i < prow; i++)
{
for (int j = (pcol+1); j < cols; j++)
{
Matrix[i][j] = (Orig_Matrix[i][j] -
((Orig_Matrix[prow][j] * Orig_Matrix[i][pcol]) /
Orig_Matrix[prow][pcol]));
}
}

for (int i = (prow+1); i < rows; i++)
{
for (int j = (pcol+1); j < cols; j++)
{
Matrix[i][j] = (Orig_Matrix[i][j] -
((Orig_Matrix[prow][j] * Orig_Matrix[i][pcol]) /
Orig_Matrix[prow][pcol]));
}
}

}

int Feasible(int rows, int cols)
{
int is_feasible = 1;

for (int i = 0; i < (rows-1); i++)
{
if ( Matrix[i][cols-1] < 0 )
is_feasible = 0;
}

return is_feasible;
}

int Optimal(int rows, int cols)
{
int is_optimal = 1;

for (int j = 0; j < (cols-1); j++)
{
if ( Matrix[rows-1][j] > 0 )
is_optimal = 0;
}

if ( !Feasible(rows, cols) )
is_optimal = 0;

return is_optimal;
}

{
int row_greaterthan_zero = 0, is_bad_row = 0;

for (int i = 0; i < (rows-1); i++)
{
row_greaterthan_zero = 1;
for (int j = 0; j < (cols-1); j++)
{
if ( Matrix[i][j] < 0 )
row_greaterthan_zero = 0;
}
if ( row_greaterthan_zero && (Matrix[i][cols-1] < 0) )
}

}

{
int col_lessthan_zero = 0, is_bad_column = 0;

for (int j = 0; j < (cols-1); j++)
{
col_lessthan_zero = 1;
for (int i = 0; i < (rows-1); i++)
{
if ( Matrix[i][j] > 0 )
col_lessthan_zero = 0;
}
if ( col_lessthan_zero && Matrix[rows-1][j] > 0 )
}

}

void Phase1_Pick_Pivot(int rows, int cols, int & prow, int & pcol)
{
int init_pcol, init_prow;
int done = 0;
Fraction Curr_Min_Ratio;
int Curr_Min_Row;

for (int i = 0; i < (rows-1); i++)
{
if ( !done && Matrix[i][cols-1] < 0 )
{
init_prow = i;
for (int j = 0; j < (cols-1); j++)
{
if ( Matrix[init_prow][j] < 0 )
{
init_pcol = j;
done = 1;
}
}
}
}

pcol = init_pcol;  // found pivot column
Curr_Min_Ratio = Matrix[init_prow][cols-1] / Matrix[init_prow][init_pcol];
Curr_Min_Row = init_prow;

for (int i = 0; i < init_prow; i++)
{
if ( Matrix[i][cols-1] >= 0 && Matrix[i][init_pcol] > 0 )
{
if ( (Matrix[i][cols-1] / Matrix[i][init_pcol]) < Curr_Min_Ratio )
{
Curr_Min_Ratio = Matrix[i][cols-1] / Matrix[i][init_pcol];
Curr_Min_Row = i;
}
}
}

prow = Curr_Min_Row;
}

void Phase2_Pick_Pivot(int rows, int cols, int & prow, int & pcol)
{
int start = 1;
Fraction Curr_Min_Ratio;
int Curr_Min_Row;

for (int j = 0; j < (cols-1); j++)
{
if ( Matrix[rows-1][j] > 0 )
pcol = j;
}

for (int i = 0; i < (rows-1); i++)
{
if ( Matrix[i][pcol] > 0 )
{
if (start)
{
Curr_Min_Ratio = Matrix[i][cols-1] / Matrix[i][pcol];
Curr_Min_Row = i;
start = 0;
}
else
{
if ( (Matrix[i][cols-1] / Matrix[i][pcol]) < Curr_Min_Ratio )
{
Curr_Min_Ratio = Matrix[i][cols-1] / Matrix[i][pcol];
Curr_Min_Row = i;
}
}
}
}

prow = Curr_Min_Row;
}

// Main Program Starts Here
int main()
{
int rows = 0, cols = 0;

Enter_Dimensions(rows, cols);
cout << endl;
Enter_Matrix(rows, cols);
cout << endl;
cout << "Starting Matrix: " << endl;
Print_Matrix(rows, cols);
cout << endl;

// Perform Phase One
cout << "Phase 1: " << endl;
cout << "-------- " << endl;
cout << endl;
while ( !Feasible(rows, cols) )
{
int pivot_row, pivot_col;

{
cout << "Problem is infeasible." << endl;
return 1;
}

Phase1_Pick_Pivot(rows, cols, pivot_row, pivot_col);

Pivot(rows, cols, pivot_row, pivot_col);

cout << "Pivoted on ROW = " << (pivot_row+1)
<< ", COLUMN = " << (pivot_col+1) << endl;
Print_Matrix(rows, cols);
cout << endl;
}

// Perform Phase Two
cout << "Phase 2: " << endl;
cout << "-------- " << endl;
cout << endl;
while ( !Optimal(rows, cols) )
{
int pivot_row, pivot_col;

{
cout << "Problem is unbounded." << endl;
return 1;
}

Phase2_Pick_Pivot(rows, cols, pivot_row, pivot_col);

Pivot(rows, cols, pivot_row, pivot_col);

cout << "Pivoted on ROW = " << (pivot_row+1)
<< ", COLUMN = " << (pivot_col+1) << endl;
Print_Matrix(rows, cols);
cout << endl;
}

cout << "---------------------------" << endl;
cout << "Solution (max) = " << -(Matrix[rows-1][cols-1]) << endl;
cout << "---------------------------" << endl;

return 0;
}