mlyacc/examples/pascal/test/t1.p

   1 junk simplex(input, output);
   2
   3 { two-phase simplex algorithm: version Feb. 24, 1988 }
   4
   5 { copyright K. Steiglitz }
   6 { Computer Science Dept. }
   7 { Princeton University 08544 }
   8 { ken@princeton.edu }
   9
  10 var
  11     maxpivots = 1000;   { maximum no. of pivots }
  12     large = 1.0e+31;    { large number used in search for minimum cost column }
  13     lowlim = -1.0e+31;  { large negative number to test for unboundedness }
  14     mmax = 32;          { max. no. of rows }
  15     ncolmax = 50;       { max. no. of columns allowed in tableau }
  16     eps = 1.0e-8;       { for testing for zero }
  17
  18 const
  19     done, unbounded, optimal: boolean;          { flags for simplex }
  20     result: (toomanycols, unbound, infeas, toomanypivots, opt);
  21     m: 1..mmax;         { no. of rows - 1, the rows are numbered 0..m }
  22     numpivots: integer; { pivot count }
  23     pivotcol, pivotrow: integer;        { pivot column and row }
  24     pivotel: real;                      { pivot element }
  25     cbar: real;         { price when searching for entering column }
  26     carry: array[-1..mmax, -1..mmax] of real;   { inverse-basis matrix of the
  27                                                   revised simplex method }
  28     phase: 1..2;        { phase }
  29     price: array[0..mmax] of real;      { shadow prices = row -1 of carry =
  30                                           -dual variables }
  31     basis: array[0..mmax] of integer;  { basis columns, negative integers
  32                                           artificial }
  33     ncol: 1..ncolmax;                   { number of columns }
  34     tab: array[0..mmax, 1..ncolmax] of real;    { tableau }
  35     lhs: array[0..mmax] of real;        { left-hand-side }
  36     d: array[1..ncolmax] of real;       { current cost vector }
  37     c: array[1..ncolmax] of real;       { cost vector in original problem }
  38     curcol: array[-1..mmax] of real;    { current column }
  39     curcost: real;                      { current cost }
  40     i, col, row: integer;               { miscellaneous variables }
  41
  42 procedure columnsearch;
  43 { looks for favorable column to enter basis.
  44   returns lowest cost and its column number, or turns on the flag optimal }
  45
  46 var
  47   i , col : integer;
  48   tempcost: real;           { minimum cost, temporary cost of column }
  49
  50   begin  { columnsearch }
  51     for i:= 0 to m do price[i]:= -carry[-1, i];  { set up price vector }
  52     optimal:= false;
  53     cbar:= large;
  54     pivotcol:= 0;
  55     for col:= 1 to ncol do
  56       begin
  57         tempcost:= d[col];
  58         for i:= 0 to m do tempcost:= tempcost - price[i]*tab[i, col];
  59         if( cbar > tempcost ) then
  60           begin
  61            cbar:= tempcost;
  62            pivotcol:= col
  63           end
  64       end;  { for col }
  65     if ( cbar > -eps ) then optimal:= true
  66   end;   { columnsearch }
  67
  68
  69 procedure rowsearch;
  70 {  looks for pivot row. returns pivot row number,
  71    or turns on the flag unbounded }
  72
  73 var
  74   i, j: integer;
  75   ratio, minratio: real;        { ratio and minimum ratio for ratio test }
  76
  77   begin  { rowsearch }
  78     for i:= 0 to m do           { generate column }
  79       begin
  80         curcol[i]:= 0.0;        { current column = B inverse * original col. }
  81         for  j:= 0 to m do curcol[i]:=
  82                            curcol[i] + carry[i, j]*tab[j, pivotcol]
  83       end;
  84   curcol[-1]:= cbar;            { first element in current column }
  85   pivotrow:= -1;
  86   minratio:= large;
  87   for i:= 0 to m do                             { ratio test }
  88     begin
  89       if( curcol[i] > eps ) then
  90         begin
  91           ratio:= carry[i, -1]/curcol[i];
  92             if( minratio > ratio ) then         { favorable row }
  93               begin
  94                 minratio:= ratio;
  95                 pivotrow:= i;
  96                 pivotel:= curcol[i]
  97               end
  98             else { break tie with max pivot }
  99               if ( (minratio = ratio) and (pivotel < curcol[i]) ) then
 100                   begin
 101                     pivotrow:= i;
 102                     pivotel:= curcol[i]
 103                   end
 104         end  { curcol > eps }
 105       end;  { for i }
 106     if ( pivotrow = -1 ) then  unbounded:= true         { nothing found }
 107                          else  unbounded:= false
 108   end;  { rowsearch }
 109
 110
 111 procedure pivot;
 112 { pivots }
 113
 114   var
 115     i, j: integer;
 116
 117   begin { pivot }
 118     basis[pivotrow]:= pivotcol;
 119     for j:= -1 to m do carry[pivotrow, j]:= carry[pivotrow, j]/pivotel;
 120     for i:= -1 to m do
 121       if( i<> pivotrow ) then
 122         for j:= -1 to m do
 123           carry[i, j]:= carry[i, j] - carry[pivotrow, j]*curcol[i];
 124     curcost:= -carry[-1, -1]
 125   end;  { pivot }
 126
 127
 128 procedure changephase;
 129 { changes phase from 1 to 2, by switching to original cost vector }
 130
 131   var
 132     i, j, b: integer;
 133
 134   begin  { changephase }
 135     phase:= 2;
 136     for i:= 0 to m do if( basis[i] <= 0 ) then
 137       writeln( '...artificial basis element ', basis[i]:5,
 138                ' remains in basis after phase 1');
 139     for j:= 1 to ncol do d[j]:= c[j];   { switch to original cost vector }
 140     for j:= -1 to m do
 141       begin
 142         carry[-1, j]:= 0.0;
 143         for i:= 0 to m do
 144           begin
 145             b:= basis[i];       { ignore artificial basis elements that are }
 146             if( b >= 1 ) then   { still in basis }
 147               carry[-1, j]:= carry[-1, j] - c[b]*carry[i,j];
 148           end  { for i }
 149       end;  { for j }
 150     curcost:= -carry[-1, -1]
 151   end;   { changephase }
 152
 153 procedure setup;
 154 { sets up test problem, lhs = tab*x, x >= 0, min c*x }
 155 { nrow = number of rows; ncol = number of cols }
 156 { tab = tableau; lhs = constants }
 157
 158 var
 159   i, j, nrow: integer;
 160
 161 begin { setup }
 162   readln(nrow);         { read number of rows }
 163   readln(ncol);         { read number of columns }
 164   m:= nrow - 1;         { rows are numbered 0..m }
 165   for j:= 1 to ncol do
 166    read(c[j]);          { cost vector }
 167   for i:= 0 to m do
 168    begin
 169     read(lhs[i]);       { left-hand-side }
 170     for j:= 1 to ncol do
 171      read(tab[i, j]);    { tableau }
 172    end;
 173
 174   done:= false;                      { initialize carry matrix, etc. }
 175   phase:= 1;
 176   for i:= -1 to m do for j:= -1 to mmax do carry[i, j]:= 0.0;
 177   for i:= 0 to m do carry[i, i]:= 1.0;            { artificial basis }
 178   for i:= 0 to m do
 179     begin
 180       carry[i, -1]:= lhs[i];      { -1 col of carry = left-hand-side }
 181       carry[-1, -1]:= carry[-1, -1] - lhs[i]        { - initial cost }
 182     end;
 183   curcost:= -carry[-1, -1];
 184   for i:= 0 to m do basis[i]:= -i;       { initial, artificial basis }
 185   if( ncol <= ncolmax ) then               { check number of columns }
 186     for col:= 1 to ncol do      { initialize cost vector for phase 1 }
 187       begin
 188         d[col]:= 0.0;
 189         for row:= 0 to m do d[col]:= d[col] - tab[row, col]
 190       end
 191   else
 192     begin
 193       writeln('...termination: too many columns for storage');
 194       done:= true;
 195       result:= toomanycols
 196     end;
 197   numpivots:= 0;
 198 end; { setup }
 199
 200
 201 begin  { simplex }
 202   setup;
 203   while( (numpivots < maxpivots) and (not done) and
 204          ( (curcost > lowlim) or (phase = 1) ) ) do
 205     begin
 206       columnsearch;
 207       if( not optimal ) then
 208         begin                         { not optimal }
 209           rowsearch;
 210             if( unbounded ) then
 211               begin
 212                 done:= true;
 213                 result:= unbound;
 214                 writeln('problem is unbounded')
 215               end
 216             else
 217               begin
 218                 pivot;
 219                 numpivots:= numpivots + 1;
 220                 if ( (numpivots = 1 ) or ( numpivots mod 10 = 0 ) ) then
 221                       writeln('pivot ', numpivots:4, ' cost= ', curcost:12)
 222               end
 223         end  { not optimal }
 224       else                            { optimal }
 225           if( phase = 1 ) then
 226             begin
 227               if( curcost > eps ) then
 228                 begin
 229                   done:= true;
 230                   result:= infeas;
 231                   writeln('problem is infeasible')
 232                 end
 233               else
 234                 begin
 235                   if ( (numpivots <> 1 ) and ( numpivots mod 10 <> 0 ) ) then
 236                     writeln('pivot ', numpivots:4, ' cost= ', curcost:12);
 237                   writeln('phase 1 successfully completed');
 238                   changephase
 239                 end
 240             end  { if phase = 1 }
 241           else
 242             begin
 243               if ( (numpivots <> 1 ) and ( numpivots mod 10 <> 0 ) ) then
 244                 writeln('pivot ', numpivots:4, ' cost= ', curcost:12);
 245               writeln('phase 2 successfully completed');
 246               done:= true;
 247               result:= opt
 248             end
 249     end;  { while }
 250   if(((curcost <= lowlim) and (phase = 2) ) then
 251     begin
 252       if ( (numpivots <> 1 ) and ( numpivots mod 10 <> 0 ) ) then
 253         writeln('pivot ', numpivots:4, ' cost= ', curcost:12);
 254       result:= unbound;
 255       writeln('problem is unbounded')
 256     end;
 257   if ( numpivots >= maxpivots ) then
 258     begin
 259       writeln('...termination: maximum number of pivots exceeded');
 260       result:= toomanypivots,
 261     end;
 262
 263   if result = opt then
 264     begin
 265       writeln('optimal solution reached');
 266       writeln('cost    =', -carry[-1,-1]:10:6);
 267       for i:= 0 to m do
 268         writeln('x(', basis[i]:4, ')= ', carry[i,-1]:10:6)
 269
 270 end.