PROCEDURE gaussj(VAR a: matice;VAR b: vektor2);
VAR
   big,dum,pivinv: double;
   i,icol,irow,j,k,l,ll: integer;
   indxc,indxr,ipiv:vektor3;
BEGIN
   FOR j := 1 TO st DO ipiv[j] := 0;
   FOR i := 1 TO st DO BEGIN
      big := 0.0;
      FOR j := 1 TO st DO
         IF ipiv[j] <> 1 THEN
            FOR k := 1 TO st DO
               IF ipiv[k] = 0 THEN
                  IF abs(a[j,k]) >= big THEN BEGIN
                     big := abs(a[j,k]);
                     irow := j;
                     icol := k;
                  END
               ELSE IF ipiv[k] > 1 THEN BEGIN
                  writeln('pause 1 in GAUSSJ - singular matrix');
                  readln
               END;
      ipiv[icol] := ipiv[icol]+1;

      { Odkomentovani vyradi pivoting }
      icol:=i;
      irow:=i;

      IF irow <> icol THEN BEGIN
         FOR l := 1 TO st DO BEGIN
            dum := a[irow,l];
            a[irow,l] := a[icol,l];
            a[icol,l] := dum;
         END;
         dum := b[irow];
         b[irow] := b[icol];
         b[icol] := dum;
      END;
      indxr[i] := irow;
      indxc[i] := icol;
      IF a[icol,icol] = 0.0 THEN BEGIN
         writeln('pause 2 in GAUSSJ - singular matrix');
         readln;
      END;
      pivinv := 1.0/a[icol,icol];
      a[icol,icol] := 1.0;
      FOR l := 1 TO st DO
         a[icol,l] := a[icol,l]*pivinv;
         b[icol] := b[icol]*pivinv;
      FOR ll := 1 TO st DO
         IF ll <> icol THEN BEGIN
            dum := a[ll,icol];
            a[ll,icol] := 0.0;
            FOR l := 1 TO st DO
               a[ll,l] := a[ll,l]-a[icol,l]*dum;
               b[ll] := b[ll]-b[icol]*dum;
         END
   END;
   FOR l := st DOWNTO 1 DO
      IF indxr[l] <> indxc[l] THEN
         FOR k := 1 TO st DO BEGIN
            dum := a[k,indxr[l]];
            a[k,indxr[l]] := a[k,indxc[l]];
            a[k,indxc[l]] := dum;
         END;
END;