Here is a Pascal program to solve small problems using the simplex algorithm.
From: Stephen Gale 
Newsgroups: sci.op-research
Subject: Re: SIMPLEX code for PC
Date: Sun, 5 May 1996 07:21:50 -0600
Attached below is a fairly simple Simplex Solver (written for
Turbo Pascal 3.0).  Please let me know of how you use the
program, as well as any problems and any successes.
=====================================================================
Stephen F. Gale, B.Sc., CCP, I.S.P.
SFGale@FreeNet.Calgary.Ab.Ca
http://www.freenet.calgary.ab.ca/~sfgale/
=====================================================================
95/96 Program Director
Association for Systems Management - Calgary Chapter
http://www.freenet.calgary.ab.ca/populati/communit/asmcal/asmcal.html
=====================================================================
************************
*** SIMPLEX SOFTWARE ***
************************
program linearoptimization;
const rowmx    = 72;
      colmx    = 112;
      mxval    = 1.0E+35;
      zero     = 0.0;
      eqzero   = 0.00000001;
var   matrix   : array [0..rowmx, 0..colmx] of real;
      basis    : array [1..rowmx] of integer;
      basisp   : array [1..rowmx] of integer;
      minmax   : real;
      error    : integer;
      name     : string[70];
      filename : string[14];
      ncon     : integer; {number of constraints}
      nvar     : integer; {number of variables}
      nltcon   : integer; {number of less than constraints}
      neqcon   : integer; {number of equal to constraints}
      ngtcon   : integer; {number of greater than contraints}
      trows1   : integer;
      trows2   : integer;
      tcols1   : integer;
      tcols2   : integer;
      tcols3   : integer;
function getmatrix( row, col : integer ) : real;
begin
   getmatrix := matrix[ row, col ];
end;
procedure putmatrix( row, col : integer; value : real );
begin
   matrix[ row, col ] := value;
end;
procedure initmatrix;
var row, col : integer;
begin
   for row := 0 to rowmx do
      for col := 0 to colmx do
         putmatrix( row, col, zero );
end;
procedure price( var xcol : integer; trow : integer; var error : integer );
var quant, val : real;
    col        : integer;
begin
   quant := -eqzero;
   for col := 1 to tcols3 do
      begin
         val := getmatrix( trow, col );
         if ( val < quant ) then
            begin
               xcol := col;
               quant := val;
            end;
      end;
   error := 0;
   if ( quant = -eqzero ) then
      error := 1;
end;
procedure leave( var xrow : integer; xcol : integer; var error : integer );
var quant, val : real;
    row        : integer;
begin
   quant := mxval;
   for row := 1 to ncon do
      begin
         val := getmatrix( row, xcol );
         if ( val > eqzero ) then
            begin
               val := getmatrix( row, tcols2 ) / val;
               if ( val < quant ) then
                  begin
                     xrow := row;
                     quant := val;
                  end;
            end;
      end;
   error := 0;
   if ( quant = mxval ) then
      error := 2;
end;
procedure pivot( xrow, xcol : integer );
var value, val, vl : real;
    row, col       : integer;
begin
   value := getmatrix( xrow, xcol );
   for row := 1 to trows2 do
      if ( row <> xrow ) then
         begin
            vl := getmatrix( row, xcol );
            for col := 1 to tcols2 do
               if ( col <> xcol ) then
                  begin
                     val := getmatrix( row, col ) - vl * getmatrix( xrow, col ) / value;
                     if ( abs( val ) < eqzero ) then
                        val := zero;
                     putmatrix( row, col, val );
                  end;
         end;
   for col := 1 to tcols2 do
      putmatrix( xrow, col, getmatrix( xrow, col ) / value );
   for row := 1 to trows2 do
      putmatrix( row, xcol, zero );
   putmatrix( xrow, xcol, 1.0);
   basis[ xrow ] := xcol;
end;
procedure optimize( trow : integer; var error : integer);
var xrow, xcol, iterate : integer;
begin
   repeat
      price( xcol, trow, error );
      if ( error = 0 ) then
         leave( xrow, xcol, error );
      if ( error = 0 ) then
         pivot( xrow, xcol );
   until ( error <> 0 )
end;
procedure simplex( var error : integer );
var val      : real;
    row, col : integer;
    flag     : boolean;
label 1000;
begin
   if ( ncon <> nltcon ) then
      begin
         optimize( trows1, error );
         if ( error > 1 ) then exit;
         error := 3;
         for row := 1 to ncon do
            if ( basis[ row ] > tcols3 ) then
               begin
                  if ( getmatrix( row, tcols2 ) > eqzero ) then exit;
                  flag := false;
                  col := 1;
                  repeat
                     if ( abs( getmatrix( row, col ) ) >= eqzero ) then
                        begin
                           pivot( row, col );
                           flag := true;
                        end;
                     col := col + 1;
                  until ( (flag) or (col > tcols3) );
               end;
      end;
   error := 0;
   optimize( trows2, error );
end;
procedure reader( var error : integer );
var row, col, column : integer;
    value, amt       : real;
    filevar          : text;
begin
   error := 0;
   writeln(con,'Problem file should be in the following format:');
   writeln(con,' Line 1 : Up to 70 character problem description');
   writeln(con,' Line 2 : +1 (for max) or -1 (for min); # of constraints; # of variables');
   writeln(con,' Line 3 : # of <= constraints; # of = constraints; # of >= constraints');
   writeln(con,' Next   : Constraints coefficients and RHS value for each constraint');
   writeln(con,' Last   : Objective function coefficients');
   writeln(con);
   write  (con,'Enter the filename containing the problem: ');
   readln (con,filename);
   assign(filevar,filename);
   reset(filevar);
   { read the problem description }
   readln( filevar, name );
   { read the minmax, number of constraints, number of variables }
   readln( filevar, minmax, ncon, nvar );
   minmax := -minmax;
   { read the number of less than, equal to, greater than contraints}
   readln( filevar, nltcon, neqcon, ngtcon );
   if ( ncon <> nltcon + neqcon + ngtcon ) then error := -1;
   trows1 := ncon + 1;
   trows2 := ncon + 2;
   tcols1 := nvar + ncon + ngtcon;
   tcols2 := tcols1 + 1;
   tcols3 := nvar + nltcon + ngtcon;
   {prepare matrix and basis}
   for row := 1 to trows2 do
      for col := 1 to tcols2 do
         putmatrix( row, col, zero );
   for row := 1 to ncon do
      basis[ row ] := 0;
   {prepare artificial and surplus variables}
   for row := 1 to ncon do
      if ( row <= nltcon ) then
         begin
            column := nvar + row;
            basis[ row ] := column;
            putmatrix( row, column, +1.0 );
         end
      else
         begin
            column := nvar + ngtcon + row;
            basis[ row ] := column;
            putmatrix( row, column, +1.0 );
            if ( row > nltcon + neqcon ) then
               begin
                  column := nvar - neqcon + row;
                  putmatrix( row, column, -1.0 );
                  putmatrix( trows1, column, +1.0 );
               end
         end;
   {read matrix and right hand side}
   for row := 1 to ncon do
      begin
         for col := 1 to nvar do
            begin
               read( filevar, value );
               putmatrix( row, col, value );
            end;
         read( filevar, value );
         putmatrix( row,      0, value );
         putmatrix( row, tcols2, value );
      readln( filevar );
      end;
   { read the coefficients of the objective function }
   for col := 1 to nvar do
      begin
         read( filevar, value);
         putmatrix(      0, col, value * minmax );
         putmatrix( trows2, col, value * minmax );
      end;
   readln( filevar );
   { calculate artifical variables }
   for col := 1 to nvar do
      begin
         value := zero;
         for row := nltcon+1 to ncon do
            value := value - getmatrix( row, col );
         putmatrix( trows1, col, value );
      end;
   close(filevar);
end;
procedure stats;
begin
   writeln;
   writeln(' * Your Variables : 1 through ', nvar);
   if ( nltcon > 0 ) then
      writeln(' * Slack Variables : ',nvar+1,' through ',nvar+nltcon);
   if ( ngtcon > 0 ) then
      writeln(' * Surplus Variables : ',nvar+nltcon+1,' through ',tcols3);
   if ( nltcon <> ncon ) then
      writeln(' * Artificial Variables : ',tcols3+1,' through ',tcols1);
end;
procedure setbasis;
var row, col : integer;
    flag     : boolean;
begin
   for col := 1 to nvar+ncon do
      begin
         flag := false;
         row := 1;
         repeat
            if ( basis[ row ] = col ) then
               flag := true
            else
               row := row + 1;
         until ( (flag) or (row > ncon) );
         if (flag) then
            basisp[ col ] := row
         else
            basisp[ col ] := 0;
      end;
end;
procedure problem;
var row, col : integer;
begin
   {filename and problem description}
   writeln('Filename: ',filename);
   writeln(' Problem: ',name);
   writeln;
   {objective function}
   if minmax < 0 then
      writeln('Maximize:')
   else
      writeln('Minimize:');
   for col := 1 to nvar do
      begin
         write(minmax*getmatrix( trows2, col ):18:8,' * Var#',col:3);
         if col <> nvar then
            write(' + ');
         writeln;
      end;
   writeln;
   {constraints}
   writeln('Subject to:');
   for row := 1 to ncon do
      begin
         for col := 1 to nvar do
            begin
               if (col = 1) then
                  write('  Constraint #',row:3,'...')
               else
                  write(' ':20);
               write(getmatrix( row, col):18:8,' * Var#',col:3);
               if col <> nvar then
                  writeln(' + ');
            end;
         if (row <= nltcon) then
            write(' <= ')
         else if (row > nltcon+neqcon) then
            write(' >= ')
         else
            write(' =  ');
         writeln(getmatrix( row, tcols2 ):18:8);
      end;
end;
procedure answer;
var smallpos         : real;
    smallptr         : integer;
    largeneg         : real;
    largeptr         : integer;
    value            : real;
    row, col, row1   : integer;
begin
      setbasis;
      stats;
      {objective value}
      writeln;
      writeln('* Value of Objective Function: ',
               -minmax*getmatrix( trows2, tcols2 ):18:8);
      writeln;
      writeln('* Variable Analysis *');
      writeln('Var':3,'Value':18,'Reduced Cost':18);
      writeln('---':3,'-----':18,'------------':18);
      for col := 1 to nvar do
         begin
            write(col:3);
            if (basisp[ col ] <> 0) then
               write(getmatrix( basisp[ col ], tcols2 ):18:8)
            else
               write(0.0:18:8);
            write(-minmax * getmatrix( trows2, col ):18:8);
            writeln;
         end;
      writeln;
      writeln('* Constraint Analysis *');
      writeln('Row':3,'RHS Value':18,'Slack/Surplus':18,'Shadow Price':18);
      writeln('---':3,'---------':18,'-------------':18,'------------':18);
      for row := 1 to ncon do
         begin
            if (row <= nltcon) then
               col := nvar + row
            else if (row > nltcon+neqcon) then
               col := nvar + row - neqcon
            else
               col := nvar + ngtcon + row;
            write(row:3);
            write(getmatrix( row, 0 ):18:8);
            if (basisp[ col ] <> 0) then
               write(getmatrix( basisp[ col ], tcols2 ):18:8)
            else
               write(0.0:18:8);
            write(-minmax * getmatrix( trows2, col ):18:8);
            writeln;
         end;
      writeln(' ');
      writeln('* Sensitivity Analysis - Right Hand Side Ranging *');
      writeln('Row':3,'Lower':18,'Current':18,'Upper':18,'OutLo':6,'OutUp':6);
      writeln('---':3,'-----':18,'-------':18,'-----':18,'-----':6,'-----':6);
      for row1 := 1 to ncon do
         begin
            if (row1 <= nltcon) then
               col := nvar + row1
            else if (row1 > nltcon+neqcon) then
               col := nvar + row1 - neqcon
            else
               col := nvar + ngtcon + row1;
            smallpos := +mxval;
            smallptr := 0;
            largeneg := -mxval;
            largeptr := 0;
            for row := 1 to ncon do
               begin
                  value := getmatrix( row, col );
                  if ( value <> zero ) then
                     begin
                        value := getmatrix( row, tcols2 ) / value;
                        if (value>zero) and (value