Program DemLU;

USES Crt;

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

PROCEDURE lubksb(VAR a: RealArrayNPbyNP;
                     n: integer;
              VAR indx: IntegerArrayNP;
                 VAR b: RealArrayNP);
{
 Resi system rovnic s pravou stranou B kdyz vstupni matice A
 byla uz prevedena do LU tvaru procedurou LUDCMP,
 pole INDX obsahuje permutace radku provedene v LUDCMP.
 Pri vystupu je ve vektoru B reseni rovnic.
 Bere v potaz prave strany s mnoha 0 na pocatku a tak je efektivni
 i pri vypoctu inverzni matice.
}
VAR
   j,ip,ii,i: integer;
   sum: real;
   bpom: RealArrayNP;

BEGIN
   ii := 0;
   FOR i := 1 TO n DO BEGIN
      ip := indx[i];
      sum := b[ip];
      b[ip] := b[i];
      bpom[ip] := b[i];
      bpom[i] := sum;
      IF ii <> 0 THEN
         FOR j := ii TO i-1 DO
            sum := sum-a[i,j]*b[j]
      ELSE IF sum <> 0.0 THEN
         ii := i;
      b[i] := sum
   END;
   Writeln(' Tisk uvnitr reseni');
   Writeln(' Prehozeny vektor b');
   FOR i:= 1 TO n DO WRITE(bpom[i]:8:3); Writeln;
   Writeln(' Tisk reseni y rovnice Ly = b');
   FOR i:= 1 TO n DO WRITE(b[i]:8:3); Writeln;
   FOR i := n DOWNTO 1 DO BEGIN
      sum := b[i];
      FOR j := i+1 TO n DO
         sum := sum-a[i,j]*b[j];
      b[i] := sum/a[i,i]
   END;
   Writeln(' Tisk reseni x rovnice Ux = y');
   FOR i:= 1 TO n DO WRITE(b[i]:8:3); Writeln;
   Writeln(' Konec zpetneho reseni'); Writeln;
END;

PROCEDURE ludcmp(VAR a: RealArrayNPbyNP;
                     n: integer;
              VAR indx: IntegerArrayNP;
                 VAR d: real);
{
 Provadi LU rozklad matice a[n,n]
 na vystupu je matice A prevedena do tvaru
               na diagonale a nad diagonalou jsou prvky beta matice U
               pod diagonalou jsou prvky alpha matice L
 v poli indx[n] jsou zaznamenany radkove permutace A provadene pri
               castecnem hledani hlaniho prvku
 d = -+ 1 podle licheho ci sudeho poctu permutaci radku
}

CONST
   tiny = 1.0e-20;
VAR
   k,j,imax,i: integer;
   sum,dum,big: real;
   vv: ^RealArrayNP;
   mode:Word;
BEGIN
   new(vv);
   d := 1.0;
   FOR i := 1 TO n DO BEGIN
      big := 0.0;
      FOR j := 1 TO n DO
         IF abs(a[i,j]) > big THEN big := abs(a[i,j]);
      IF big = 0.0 THEN BEGIN
         writeln('pause in LUDCMP - singular matrix');
         readln
      END;
      vv^[i] := 1.0/big
   END;
   mode := LastMode;
   TEXTMODE(CO80);
   TextAttr := Black + LightGray*16;
   WINDOW(5,5,5+9*n,5+n);
   Gotoxy(1,1);
   FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO Write(a[i,j]:8:3,' '); Writeln
        END;
   Window(47,5,56,5+n);
   Gotoxy(1,1);
   FOR i:=1 TO n DO Writeln(vv^[i]:8:3);
   FOR j := 1 TO n DO BEGIN
      Window(5,15+n,5+9*n,17+n); Readln;
      Gotoxy(9*j-8,1); Write(j);
      TextAttr := Red+ LightGray*16;;
      FOR i := 1 TO j-1 DO BEGIN
      WINDOW(5,22,75,23);
      GotoXY(1,1);Write('Budu menit U cast (bez diagonaly)        ');
      WINDOW(5,5,5+9*n,5+n);
         sum := a[i,j];
         FOR k := 1 TO i-1 DO
            sum := sum-a[i,k]*a[k,j];
            a[i,j] := sum; Readln;
            Gotoxy(9*j-8,i);Write(a[i,j]:8:3,' ')
      END;
      big := 0.0;
      WINDOW(5,22,75,23);
      GotoXY(1,1);Write('Budu docasne menit L cast vcetne DIAG       ');
      FOR i := j TO n DO BEGIN
         sum := a[i,j];
         FOR k := 1 TO j-1 DO
            sum := sum-a[i,k]*a[k,j];
         a[i,j] := sum; Readln;
         WINDOW(5,5,5+9*n,5+n);
         Gotoxy(9*j-8,i);Write(a[i,j]:8:3,' ');
         Window(57,5,66,5+n); GotoXY(1,i);
         dum := vv^[i]*abs(sum);
         Write(dum:8:3);
         IF dum >= big THEN BEGIN
            big := dum;
            imax := i
         END
      END;
      Window(47,15+n,56,17+n);
      Gotoxy(1,1); Write(imax);Gotoxy(2,1);Write(big:8:3);
      IF j <> imax THEN BEGIN
         WINDOW(5,22,75,23);
         GotoXY(1,1);Write('PREHAZUJI radek na diagonale s radkem IMAX  ');
         WINDOW(5,5,5+9*n,5+n);
         FOR k := 1 TO n DO BEGIN
            dum := a[imax,k];
            a[imax,k] := a[j,k];
            a[j,k] := dum;Readln;
            Gotoxy(9*k-8,imax);Write(a[imax,k]:8:3,' ');
            Gotoxy(9*k-8,j);Write(a[j,k]:8:3,' ')
         END;
         d := -d;
         Readln;Window(47,5,56,5+n);
         vv^[imax] := vv^[j];
         Gotoxy(1,imax); Write(vv^[imax]:8:3)
      END;
      indx[j] := imax;
      Readln;Window(71,5,80,5+n);
      GotoXY(1,j);Write(imax);
      IF a[j,j] = 0.0 THEN a[j,j] := tiny;
      IF j <> n THEN BEGIN
         WINDOW(5,22,75,23);
         GotoXY(1,1);Write('Delim poddiagonalni cast  /u[j,j]             ');
         dum := 1.0/a[j,j];
         FOR i := j+1 TO n DO BEGIN
            a[i,j] := a[i,j]*dum; Readln;
            WINDOW(5,5,5+9*n,5+n);
            Gotoxy(9*j-8,i);Write(a[i,j]:8:3,' ')
            END;
      END
   END;
   dispose(vv);
{   WINDOW(5,15,5+9*n,15+n);
   Gotoxy(1,1);
   FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO Write(a[i,j]:8:3,' '); Writeln
        END;  }
   Repeat
   Until KeyPressed;
   TextMode(mode)
END;

VAR a,apom,left,upper,asouc: RealArrayNPbyNP;
    n,i,j,k: integer;
    indx: IntegerArrayNP;
    b,bpom,bvys: RealArrayNP;
    detzn,dethod,pom: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];

    ludcmp(a,n,indx,detzn);
    Writeln(' Matice dekomponovana');
    FOR i:=1 TO n DO begin
        Write('  |');
        FOR j:=1 TO i-1 DO Begin Write(a[i,j]:8:3); left[i,j] := a[i,j]
        END;
        Write('       1');left[i,i] :=1.;
        FOR j:=i+1 TO n DO BEGIN Write('       0'); left[i,j] := 0.
        END;
        Write('  |      |');
        FOR j:=1 TO i-1 DO BEGIN Write('       0'); upper[i,j] := 0.
        END;
        FOR j:=i TO n DO BEGIN Write(a[i,j]:8:3); upper[i,j] := a[i,j]
        END;
        Write('  |'); Writeln
        END;
    FOR i:=1 TO n DO
        FOR j:=1 TO n DO BEGIN
               asouc[i,j] :=0.;
               for k:=1 TO n DO asouc[i,j]:= asouc[i,j]+left[i,k]*upper[k,j];
        End;
    Readln; Writeln;
    FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO Write(asouc[i,j]:8:3); Writeln
    END;
    for k:=n-1 DOWNTO 1 DO
        if (indx[k] <> k) THEN
                FOR j:= 1 TO n DO BEGIN
                        pom := asouc[k,j];
                        asouc[k,j] := asouc[indx[k],j];
                        asouc[indx[k],j] := pom
                END;
    Writeln; Writeln(' asouc s radky prehozenymi do spravne polohy');
    FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO Write(asouc[i,j]:8:3); Writeln
    END;
    FOR i:=1 TO n DO BEGIN Read(fvst,b[i]); bpom[i]:=b[i]; 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]:8:3,' |');
        END;
    Writeln; Readln;
    lubksb(a,n,indx,b);
    FOR i:=1 TO n DO Write(b[i]:9:3);
    Writeln; Writeln('         KONTROLA Spravnosti');
    FOR i:=1 TO n DO BEGIN bvys[i]:=0.;
        FOR j:=1 TO n DO bvys[i]:=bvys[i]+apom[i,j]*b[j];
        Writeln(' Zadana   ',bpom[i],'  Kontrola  ',bvys[i]);
        End;
    Readln;

END.