{Priklad reseni linearni optimalizace simplexovou metodou
3 promenne, 3 podminky, po 1 od kazdeho typu
}

CONST
   nv = 3;
   mv = 5;
   np = 4;
   mp = 7;
TYPE
   RealArrayMPbyNP = ARRAY [1..mp,1..np] OF real;
   IntegerArrayN = ARRAY [1..nv] OF integer;
   IntegerArrayM = ARRAY [1..mv] OF integer;
   IntegerArrayNP = ARRAY [1..np] OF integer;

PROCEDURE simplx(VAR a: RealArrayMPbyNP;
          m,n,m1,m2,m3: integer;
             VAR icase: integer;
             VAR izrov: IntegerArrayN;
             VAR iposv: IntegerArrayM);
{Simplexova metoda pro reseni ulohy linearniho programovasni.
Vstupni parametry N je pocet promennych, M1 je pocet podminek
 <=, M2 je pocet podminek >=, M3 pocet podminek =, M je pocet
 vsech podminek, M = M1 + M2 + M3. Matice A[1..MP,1..NP] kde
 NP >= N+1 ,  MP >= m+2. V 1. radku zadavame ucinek, na 1. miste
 absolutni clen. Vdalsich radcich jsou postupne podminky v poradi
 typ 1(<=), 2(>=), 3(=).  V radku podminky na vstupu je vzdy nejdrive
 prava strana na 1 miste a pak nasleduje N koeficientu podminky,
 koeficienty se udavaji prevedene na pravou stranu, tj. vynasobene -1.
 Radek M+2 neni na vstupu zadan a je vyhrazen pro pomocnou ucinkovou fci.
 Na vystupu A[1,1] je maximum ucinku.
 Pokud IPOSV[J] (1..M) = I <= N, pak promenna X[I] je levostranna
 a jeji hodnota v optimu je A[J+1,1].
 Pokud IPOSV[J] (1..M) = K > N, tak rozdil v K-N te podmince je A[J+1,1].
 Pokud IZROV[J] (1..N) = I <= N, pak promenna X[I] = 0 (pravostranna).
 Pokud IZROV[J] (1..N) = K > N,  pak K-N ta podminka je splnena presne.
 reprezentovane nyni sloupcem J+1 matice A.
 ICASE je na vystupu 0, pokud je nalezeno konecne reseni,
 ICASE = 1 pokud je ucinek neomezeny,
 ICASE = -1 pokud nelze vyhovet podminkam.
 Pokud iterace nekonverguje, lze zvetsit pozadovanou presnost EPS. }

LABEL 1,2,3,4,5,99;
CONST
   eps = 1.0e-6;
VAR
   nl2,nl1,m12,kp,kh,k,is,ir,ip,i: integer;
   q1,bmax: real;
   l1: ^IntegerArrayNP;
   l2,l3: ^IntegerArrayM;

PROCEDURE simp1(VAR a: RealArrayMPbyNP;
                   mm: integer;
               VAR ll: IntegerArrayNP;
             nll,iabf: integer;
               VAR kp: integer;
             VAR bmax: real);
VAR
   k: integer;
   test: real;
BEGIN
   kp := ll[1];
   bmax := a[mm+1,kp+1];
   FOR k := 2 TO nll DO BEGIN
      IF iabf = 0 THEN
         test := a[mm+1,ll[k]+1]-bmax
      ELSE
         test := abs(a[mm+1,ll[k]+1])-abs(bmax);
      IF test > 0.0 THEN BEGIN
         bmax := a[mm+1,ll[k]+1];
         kp := ll[k]
      END
   END
END;

PROCEDURE simp2(VAR a: RealArrayMPbyNP;
                  m,n: integer;
               VAR l2: IntegerArrayM;
                  nl2: integer;
               VAR ip: integer;
                   kp: integer;
               VAR q1: real);
LABEL 1,2,99;
CONST
   eps = 1.0e-6;
VAR
   k,ii,i: integer;
   qp,q0,q: real;
BEGIN
   ip := 0;
   FOR i := 1 TO nl2 DO
      IF a[l2[i]+1,kp+1] < -eps THEN GOTO 1;
   GOTO 99;
1: q1 := -a[l2[i]+1,1]/a[l2[i]+1,kp+1];
   ip := l2[i];
   FOR i := i+1 TO nl2 DO BEGIN
      ii := l2[i];
      IF a[ii+1,kp+1] < -eps THEN BEGIN
         q := -a[ii+1,1]/a[ii+1,kp+1];
         IF q < q1 THEN BEGIN
            ip := ii;
            q1 := q
         END
         ELSE IF q = q1 THEN BEGIN
            FOR k := 1 TO n DO BEGIN
               qp := -a[ip+1,k+1]/a[ip+1,kp+1];
               q0 := -a[ii+1,k+1]/a[ii+1,kp+1];
               IF q0 <> qp THEN GOTO 2
            END;
2:          IF q0 < qp THEN ip := ii
         END
      END
   END;
99:
END;

PROCEDURE simp3(VAR a: RealArrayMPbyNP;
          i1,k1,ip,kp: integer);
VAR
   kk,ii: integer;
   piv: real;
BEGIN
   piv := 1.0/a[ip+1,kp+1];
   FOR ii := 1 TO i1+1 DO BEGIN
      IF ii-1 <> ip THEN BEGIN
         a[ii,kp+1] := a[ii,kp+1]*piv;
         FOR kk := 1 TO k1+1 DO
            IF kk-1 <> kp THEN
               a[ii,kk] := a[ii,kk] -a[ip+1,kk]*a[ii,kp+1]
      END
   END;
   FOR kk := 1 TO k1+1 DO
      IF kk-1 <> kp THEN
         a[ip+1,kk] := -a[ip+1,kk]*piv;
   a[ip+1,kp+1] := piv
END;

BEGIN
   IF m <> m1+m2+m3 THEN BEGIN
      writeln('pause in routine SIMPLX');
      writeln('bad input constraint counts');
      readln
   END;
   new(l1);
   new(l2);
   new(l3);
   nl1 := n;
   FOR k := 1 TO n DO BEGIN
      l1^[k] := k;
      izrov[k] := k
   END;
   nl2 := m;
   FOR i := 1 TO m DO BEGIN
      IF a[i+1,1] < 0.0 THEN BEGIN
         writeln('pause in routine SIMPLX');
         writeln('bad input tableau');
         readln
      END;
      l2^[i] := i;
      iposv[i] := n+i
   END;
   FOR i := 1 TO m2 DO l3^[i] := 1;
   ir := 0;
   IF m2+m3 = 0 THEN
      GOTO 5;
   ir := 1;
   FOR k := 1 TO n+1 DO BEGIN
      q1 := 0.0;
      FOR i := m1+1 TO m DO
         q1 := q1+a[i+1,k];
      a[m+2,k] := -q1
   END;
3: simp1(a,m+1,l1^,nl1,0,kp,bmax);
   IF (bmax <= eps) AND (a[m+2,1] < -eps) THEN BEGIN
      icase := -1;
      GOTO 99
   END
   ELSE IF (bmax <= eps)
      AND (a[m+2,1] <= eps) THEN BEGIN
      m12 := m1+m2+1;
      IF m12 <= m THEN BEGIN
         FOR ip := m12 TO m DO BEGIN
            IF iposv[ip] = ip+n THEN BEGIN
               simp1(a,ip,l1^,nl1,1,kp,bmax);
               IF bmax > 0.0 THEN GOTO 1
            END
         END
      END;
      ir := 0;
      m12 := m12-1;
      IF m1+1 > m12 THEN GOTO 5;
      FOR i := m1+1 TO m12 DO
         IF l3^[i-m1] = 1 THEN
            FOR k := 1 TO n+1 DO a[i+1,k] := -a[i+1,k];
      GOTO 5
   END;
   simp2(a,m,n,l2^,nl2,ip,kp,q1);
   IF ip = 0 THEN BEGIN
      icase := -1;
      GOTO 99
   END;
1: simp3(a,m+1,n,ip,kp);
   IF iposv[ip] >= n+m1+m2+1 THEN BEGIN
      FOR k := 1 TO nl1 DO
         IF l1^[k] = kp THEN GOTO 2;
2:    nl1 := nl1-1;
      FOR is := k TO nl1 DO l1^[is] := l1^[is+1]
   END
   ELSE BEGIN
      IF iposv[ip] < n+m1+1 THEN GOTO 4;
      kh := iposv[ip]-m1-n;
      IF l3^[kh] = 0 THEN GOTO 4;
      l3^[kh] := 0
   END;
   a[m+2,kp+1] := a[m+2,kp+1]+1.0;
   FOR i := 1 TO m+2 DO a[i,kp+1] := -a[i,kp+1];
4: is := izrov[kp];
   izrov[kp] := iposv[ip];
   iposv[ip] := is;
   IF ir <> 0 THEN
      GOTO 3;
5: simp1(a,0,l1^,nl1,0,kp,bmax);
   IF bmax <= 0.0 THEN BEGIN
      icase := 0;
      GOTO 99
   END;
   simp2(a,m,n,l2^,nl2,ip,kp,q1);
   IF ip = 0 THEN BEGIN
      icase := 1;
      GOTO 99
   END;
   simp3(a,m,n,ip,kp);
   GOTO 4;
99:
   dispose(l3);
   dispose(l2);
   dispose(l1)
END;

VAR a: RealArrayMPbyNP;
    m,n,m1,m2,m3: integer;
    icase,i,j: integer;
    izrov: IntegerArrayN;
    iposv: IntegerArrayM;

BEGIN
{zadani ulohy}
    m1:=1; m2:=1;m3:=1;
    m:=m1+m2+m3;
    n:=3;
    a[1,2] :=1; a[1,3] :=1; a[1,4] :=1; a[1,1] :=0;
    a[2,2] :=-3; a[2,3] :=-1; a[2,4] :=-1; a[2,1] :=3;
    a[3,2] :=-1; a[3,3] :=-3; a[3,4] :=1; a[3,1] :=1;
    a[4,2] :=-2; a[4,3] :=0; a[4,4] :=-2; a[4,1] :=1;
{KONEC ZADANI}
    simplx(a,m,n,m1,m2,m3,icase,izrov,iposv);
    Writeln('case =', icase);
    Writeln('vystupni matice a');
    FOR i:=1 TO m+2 DO BEGIN
          FOR j:=1 TO n+1 DO  Write(a[i,j]:10:3,' ');
          Writeln;
          END;
    Writeln('Pisu izrov');
    FOR j:=1 TO n DO Write(izrov[j]:5);
    Writeln;
    Writeln('Pisu iposv');
    FOR i:=1 TO m DO Write(iposv[i]:5);
    Writeln;
    Readln;
{Vypis  vysledku}
    Writeln('Vypis vysledku');
    Writeln('Maximalni hodnota ucinkove funkce = ',a[1,1]:10:3);
    Writeln(' -------------------------------------------- ');
    FOR j:=1 TO n DO
         IF izrov[j] <= n THEN Writeln('x[',izrov[j]:0,'] = 0');
    FOR j:=1 TO n DO
         IF iposv[j] <= n THEN Writeln('x[',iposv[j]:0,'] = ',a[j+1,1]:10:3);
    Writeln(' -------------------------------------------- ');
    Writeln('Vypis rozdilu v jednotlivych podminkach');
    FOR j:=1 TO n DO
         IF izrov[j] > n THEN Writeln('rozdil v podm.',izrov[j]-n:0,' = 0');
    FOR j:=1 TO n DO
         IF iposv[j] > n THEN Writeln('rozdil v podm.',iposv[j]-n:0,
                      ' = ',a[j+1,1]:10:3);
    Readln;

END.