This computes 

min  c x
       A x = b
       x >=0
       
through the Simplex method

How to use this notebook

First you need input the problem data in the section Problem data. You are required to provide the matrix A, c (costs), b (rhs). Other quantities will be computed automatically, like m (# rows) and n (# columns).
Alternatively you may use the data already present in the notebook.

Then you are ready to solve Phase One. To this aim you have to set up the data for Phase One in the section Setting Phase One.
In this section you needn't provide data. Everything is computed automatically.
If you set manual to True the choice of the column to enter the base has to be done manually, otherwise is done automatically.

Now you may iterate the simplex steps to find an initial feasible solution. You do this in the section Iteration. (see instructions within the section).

Once a feasible solution is found the data must be reset to solve Phase Two.
This is done in Resetting for Phase Two. Follow carefully instructions within the section.

At this point you go back to Iteration and carry out the simplex steps until you find a solution or discover that the problem is unbounded.

Input Data

Manual instance

A = {{1,0,0,1,0, 1,0},
     {0,1,0,0,1, 2,1},
     {1,1,2,1,1,1,0},
     {0,1,3,1,0, 0,0}};

m=Length[A];n=Length[A[[1]]];

b = {4,4,2,5};

If[Length[b]!=m, Print["Check b !"]];
If[Length[c]!=n, Print["Check c !"]];

Assignment

This is a random nn x nn weighted assignment instance

nn=4;

A=Table[0,{i,2 nn},{j,nn^2}];

Do[A[[i,(j-1)nn+i]]=1;A[[nn+j,(j-1)nn+i]]=1,
{i,nn},{j,nn}];
A=Drop[A,-1]

b=Table[1,{2 nn-1}]

c=Table[Random[Integer,{1,10}],{nn^2}]

m=Length[A];n=Length[A[[1]]];

Setting Phase I

Off[General::spell1];

AA=A;
cc=c;
nn=n;
A=Transpose[Join[Transpose[A],IdentityMatrix[m]]];
c=Table[If[i<=n,0,1],{i,n+m}];
Do[If[b[[i]]<0,A[[i,n+i]]=-1],{i,m}];
n=n+m;
base = Table[i,{i,nn+1,n}];
manual=False;
resume={};
MatrixForm[A]

Resetting  for Phase II

first compute this

base

v

artinbase=Select[base,#>nn &]

If different from {} and value=0 then there are degenerate artificial vars in base.
Do the following for each one of them

artinbase

artificial index to be pushed out

indexout=artinbase[[1]];

artinbase=Drop[artinbase,1];

baseout=Position[base,indexout][[1,1]]

index to be brought in (new matrix must be non singular and only indices with non zero values in the following vector must be considered)

A

WAi=(W.A)[[baseout]]

notzero=Select[WAi, #!=0 &]

indexin=Position[WAi,notzero[[1]]][[1,1]]

base updating

base[[baseout]]=indexin;

base

Now do the following and go to Iteration

A=AA;

c=cc;
n=nn;

Iteration

Current base

base

B is the basis matrix

B = Transpose[Transpose[A][[ base ]]];
MatrixForm[B]

W = Inverse[B];
MatrixForm[W]

current x in base

x = W. b

cb = c[[ base ]]

dual variables

u = cb . W

current value

v=u . b

resume=AppendTo[resume,{base,v}];

reduced costs

cr = c - u . A

enteringset= Select[
             Transpose[{N[cr,2],Range[n]}],#[[1]] < 0 &];

variables eligible to enter the basis:
first row    -> reduced costs
second row -> variable indices

If[enteringset=={},
   Print["optimal base !"];p=0;,p=1;
   TableForm[Transpose[enteringset]]]

If[p>0,
    If[manual,p=Input["value for p"],
              p=enteringset[[1,2]]
      ]]

If[p>0,y = W . Transpose[A][[ p ]];]

x

y

Selection of variable to leave the basis

If[p>0,
   Print["{x[[i]],y[[i]],i} for y[[i]]>0 "];Print[""];
   xy=Select[Transpose[{x,y,Range[m]}], #[[2]]>0 &];
   If[xy=={},
      Print["Problem  unbounded  !"];
      ,
      Print[xy];Print[""];
      xdivy=Table[{xy[[i,1]]/xy[[i,2]],xy[[i,3]]},{i,Length[xy]}];
      Print["{x[[i]]/y[[i]],i} for y[[i]]>0 "];Print[""];
      Print[xdivy];Print[""];
      xdivyT=Transpose[xdivy];
      xdivyT[[1]]-=Min[xdivyT[[1]]];
      xdivy=Transpose[xdivyT];
      q=Select[xdivy, #[[1]]==0 &][[1,2]];
      Print["q = ",q];
   ]]

If[p>0,base[[q]] = p]

Final solution

solution=Table[0,{n}];
Do[solution[[base[[i]]]]=x[[i]],{i,m}];
solution

value=solution.c

MatrixForm[resume]