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;
  Fraction RA, RB, Answer;
 
  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;
 
  Answer.Numer = LHNumer + RHNumer;
  Answer.Denom = LHDenom;
 
  Answer.Reduce( );
 
  return Answer;
 }
 
 Fraction operator - ( const Fraction & A, const Fraction & B )
 {
  int LHNumer, RHNumer, LHDenom, RHDenom;
  Fraction RA, RB, Answer;
 
  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;
 
  Answer.Numer = LHNumer - RHNumer;
  Answer.Denom = LHDenom;
 
  Answer.Reduce( );
 
  return Answer;
 }
 
 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 )
 {
  Fraction RA, RB, Answer;
 
  RA = A;
  RB = B;
 
  RA.Reduce( );
  RB.Reduce( );
 
  Answer = Fraction ( RA.Numer * RB.Numer, RA.Denom * RB.Denom );
 
  Answer.Reduce( );
 
  return Answer;
 }
 
 Fraction operator / ( const Fraction & A, const Fraction & B )
 {
  Fraction RA, RB, Answer;
 
  RA = A;
  RB = B;
 
  RA.Reduce( );
  RB.Reduce( );
 
  Answer = Fraction ( RA.Numer * RB.Denom, RA.Denom * RB.Numer );
 
  Answer.Reduce( );
 
  return Answer;
 }
 
 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);
int Bad_Row(int rows, int cols);
int Bad_Column(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 Bad_Row(int rows, int cols)
{
 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) )
   is_bad_row = 1;
 }
 
 return is_bad_row;
}

int Bad_Column(int rows, int cols)
{
 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 )
   is_bad_column = 1;
 }
 
 return is_bad_column;
}

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;
 
  if ( Bad_Row(rows, cols) )
  {
   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;
 
  if ( Bad_Column(rows, cols) )
  {
   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;
}