Naive simplex

simplex.nb

simplex[a_, c_, bIdx1_, ab_, itr_] := 
 Module[{cols, colB, b, bInv, bIdx = bIdx1, rc, theta, j, aj, s, cb, 
   u, l, fc, thetas, n = 0},
  cols = Transpose[a];
  While[n++ < itr,
   colB = cols[[#]] & /@ bIdx;
   Print["basic indices=", bIdx];
   b = Transpose[colB];
   Print["basic matrix=", MatrixForm[b]];
   bInv = Inverse[b];
   Print["basic matrix invert=", MatrixForm[bInv]];
   s = LinearSolve[b, ab];
   Print["basic values=", MatrixForm[s]];
   cb = c[[#]] & /@ bIdx;
   fc = cb . s;
   Print["basic cost=", fc];
   rc = MapIndexed[#1 - cb . bInv . cols[[First[#2]]] &, c];
   Print["reduced cost vector=", MatrixForm[rc]];
   If[AllTrue[rc, # >= 0 &], Break[]];
   (* Use smallest indices rule. 
   Here the size of rc is always m so we can always choose the first \
position *)
   j = First[FirstPosition[rc, x_ /; x < 0]];
   Print["choose non basic index j=", j];
   aj = cols[[j]];
   u = bInv . aj;
   Print["u=", MatrixForm[u]];
   If[AllTrue[u, # <= 0 &], fc = -Infinity; Break[]];
   thetas = 
    MapIndexed[{If[#1 > 0, s[[First[#2]]]/#1, Infinity], First[#2]} &,
      u];
   Print["thetas=", MatrixForm[First /@ thetas]];
   (* Use smallest indices rule. 
   If more than one thetas are the smallest, 
   we take the one with minimum basic index, 
   so we need to take a look at bIdx. *)
   {theta, l} = 
    First[MinimalBy[MinimalBy[thetas, First[#] &], 
      bIdx[[#[[2]]]] &]];
   Print["choose theta=", theta];
   Print["choose l=", l];
   Print["replace ", bIdx[[l]], " with ", j];
   bIdx[[l]] = j;
   ];
  {bIdx, s, fc}
  ]