Program DemGaussj;

{Demonstruje reseni systemu linearnich rovnic a vypocet inverzni
matice pomoci Gauss-Jordanovy eliminace s
uplnym hledanim hlavniho prvku}

USES Crt;

CONST
   np = 20;
   mp = 20;
TYPE
   RealArrayNP = ARRAY [1..np] OF real;
   IntegerArrayNP = ARRAY [1..np] OF integer;
   RealArrayNPbyNP = ARRAY [1..np,1..np] OF real;
   RealArrayNPbyMP = ARRAY [1..np,1..mp] OF real;

PROCEDURE gaussj(VAR a: RealArrayNPbyNP;
                     n: integer;
                 VAR b: RealArrayNPbyMP;
                     m: integer);
{
 Resi system linearnich rovnic s mnoha pravymi stranami
 Gauss-Jordanovou eliminaci s uplnym hledanim hlavniho prvku.
 Matice a[n,n], m pravych stran v matici b[n,m]
 na vystupu je v a inverzni matice a^-1
 a v b jsou reseni rovnic
}
VAR
   big,dum,pivinv: real;
   i,icol,irow,j,k,l,ll: integer;
   indxc,indxr,ipiv: ^IntegerArrayNP;
BEGIN
   new(indxc);
   new(indxr);
   new(ipiv);
   FOR j := 1 TO n DO ipiv^[j] := 0;
   FOR i := 1 TO n DO BEGIN
      big := 0.0;
      FOR j := 1 TO n DO
         IF ipiv^[j] <> 1 THEN
            FOR k := 1 TO n 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;
      IF irow <> icol THEN BEGIN
         FOR l := 1 TO n DO BEGIN
            dum := a[irow,l];
            a[irow,l] := a[icol,l];
            a[icol,l] := dum
         END;
         FOR l := 1 TO m DO BEGIN
            dum := b[irow,l];
            b[irow,l] := b[icol,l];
            b[icol,l] := dum
         END
      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 n DO
         a[icol,l] := a[icol,l]*pivinv;
      FOR l := 1 TO m DO
         b[icol,l] := b[icol,l]*pivinv;
      FOR ll := 1 TO n DO
         IF ll <> icol THEN BEGIN
            dum := a[ll,icol];
            a[ll,icol] := 0.0;
            FOR l := 1 TO n DO
               a[ll,l] := a[ll,l]-a[icol,l]*dum;
            FOR l := 1 TO m DO
               b[ll,l] := b[ll,l]-b[icol,l]*dum
         END
   END;
   FOR l := n DOWNTO 1 DO
      IF indxr^[l] <> indxc^[l] THEN
         FOR k := 1 TO n DO BEGIN
            dum := a[k,indxr^[l]];
            a[k,indxr^[l]] := a[k,indxc^[l]];
            a[k,indxc^[l]] := dum
         END;
   dispose(ipiv);
   dispose(indxr);
   dispose(indxc)
END;

VAR a,apom,ainv,ajed: RealArrayNPbyNP;
    n,i,j,k,m: integer;
    indx: IntegerArrayNP;
    b,bpom,bvys: RealArrayNPbyMP;
    detzn,dethod,sum,sumi,msum,msumi,normb,normx:Real;
    fvst,fvys:Text;
    vstnam:String;
Begin
    Writeln('Zadej jmeno vstupniho souboru');
    Readln(vstnam);
    ASSIGN(fvst,vstnam);
    RESET(fvst);
    Readln(fvst,n);
    FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO Read(fvst,a[i,j]); Readln(fvst)
        END;
    FOR i:=1 TO n DO
        FOR j:=1 TO n DO
                   apom[i,j] := a[i,j];
    m :=1;
    FOR J:=1 to m DO
    FOR i:=1 TO n DO BEGIN Read(fvst,b[i,j]); bpom[i,j]:=b[i,j]; END;
    Readln(fvst);
    Writeln; Writeln('      ZADANI Ulohy');
    Writeln;
    FOR i:=1 TO n DO BEGIN
        Write('|   ');
        FOR j:=1 TO n DO Write(apom[i,j]:8:3);
        Writeln('    |   | x(',i:0,') |       | ',b[i,1]:8:3,' |');
        END;
    Writeln; Readln;
    gaussj(a,n,b,m);
    Writeln(' Vysledny vektor x ');
    FOR j:=1 TO m DO BEGIN
    FOR i:=1 TO n DO Write(b[i,j]:9:3); Writeln; END;
    Writeln; Writeln('         KONTROLA Spravnosti');
    FOR k:=1 TO m DO BEGIN
    FOR i:=1 TO n DO BEGIN bvys[i,k]:=0.;
        FOR j:=1 TO n DO bvys[i,k]:=bvys[i,k]+apom[i,j]*b[j,k];
        Writeln(' Zadana   ',bpom[i,k],'  Kontrola  ',bvys[i,k]);
        End;
        WRITELN; END;
    Readln;
    Writeln;
    normb:=0; normx:=0;
    FOR i:= 1 TO n DO BEGIN
            IF abs(bpom[i,1]) > normb THEN normb := abs(bpom[i,1]);
            IF abs(b[i,1]) > normx THEN normx := abs(b[i,1]);
    END;
    FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO Write(a[i,j]:8:3); Writeln
        END;
    Writeln;
    FOR i:=1 TO n DO
        FOR j:=1 TO n DO BEGIN
               ajed[i,j] := 0.;
               FOR k:=1 TO n DO ajed[i,j] := ajed[i,j] + apom[i,k]*a[k,j]
        END;
    Writeln;
    FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO Write(ajed[i,j],'   '); Writeln
        END;
    Readln;
    Writeln;
    msum := 0.; msumi := 0.;
    FOR i:= 1 TO n DO BEGIN
             sum:=0; sumi:=0;
             FOR j := 1 TO n DO BEGIN
                      sum := sum + abs(apom[i,j]);
                      sumi := sumi + abs(a[i,j])
             END;
             IF sum > msum THEN msum := sum;
             IF sumi > msumi THEN msumi := sumi;
    END;
    Writeln(' Radkova norma || A || =  ',msum);
    Writeln(' Radkova norma || A^-1 || =  ',msumi);
    Writeln(' Podminenost matice je || A ||* || A^-1 ||=  ',msum*msumi);
    Writeln;
    Writeln(' Radkova norma matice je souhlasna s maximovou normou vektoru');
    Writeln(' Maximova norma ||x|| = ',normx);
    Writeln(' Maximova norma ||b|| = ',normb);
    Writeln(' || A.x || = || b || =',normb,' <= || A || * ||x|| = ',
                                                    msum * normx);
    Readln
END.