Program ResLin3;

USES Crt;

CONST
   np = 50;
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;
BEGIN
   ii := 0;
   FOR i := 1 TO n DO BEGIN
      ip := indx[i];
      sum := b[ip];
      b[ip] := b[i];
      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;
   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
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;
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;
   FOR j := 1 TO n DO BEGIN
      FOR i := 1 TO j-1 DO BEGIN
         sum := a[i,j];
         FOR k := 1 TO i-1 DO
            sum := sum-a[i,k]*a[k,j];
         a[i,j] := sum
      END;
      big := 0.0;
      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;
         dum := vv^[i]*abs(sum);
         IF dum >= big THEN BEGIN
            big := dum;
            imax := i
         END
      END;
      IF j <> imax THEN BEGIN
         FOR k := 1 TO n DO BEGIN
            dum := a[imax,k];
            a[imax,k] := a[j,k];
            a[j,k] := dum
         END;
         d := -d;
         vv^[imax] := vv^[j]
      END;
      indx[j] := imax;
      IF a[j,j] = 0.0 THEN a[j,j] := tiny;
      IF j <> n THEN BEGIN
         dum := 1.0/a[j,j];
         FOR i := j+1 TO n DO
            a[i,j] := a[i,j]*dum
      END
   END;
   dispose(vv)
END;

VAR a,apom,ainv,ajed: RealArrayNPbyNP;
    n,i,j,k: integer;
    indx: IntegerArrayNP;
    b,bpom,bvys: RealArrayNP;
    detzn,dethod:Real;
    fvst,fvys:Text;
    vstnam:String;
Begin
    Writeln(' Zadej dimenzi matice ');
    Readln(n);
    Writeln(' Zadej x(1), ... x(n)');
    FOR i:=1 TO n DO Read(bpom[i]);
    Readln;
    FOR i:=1 TO n DO BEGIN
        a[i,1] := 1.;
        FOR j:=2 TO n DO a[i,j] := a[i,j-1]*bpom[i];
        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');
    Writeln(' Zadej y(1).. y(n)');
    FOR i:=1 TO n DO BEGIN Read(b[i]); bpom[i]:=b[i]; END;
    Writeln; Writeln('      ZADANI Ulohy');
    Writeln;
    FOR i:=1 TO n DO BEGIN
        Write('|   ');
        FOR j:=1 TO n DO Write(apom[i,j]:11:3);
        Writeln('    |   | x(',i:0,') |       | ',b[i]:11:3,' |');
        END;
    Writeln; Readln;
    lubksb(a,n,indx,b);
    FOR i:=1 TO n DO Write(b[i]:11);
    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;
    Writeln;
    dethod := detzn;
    FOR i:=1 TO n DO dethod := dethod*a[i,i];
    Writeln(' DETERMINANT A  = ',dethod); Writeln;
    Readln;

    FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO bpom[j]:=0.;
        bpom[i] := 1.;
        lubksb(a,n,indx,bpom);
        FOR j:=1 TO n DO ainv[j,i]:=bpom[j];
        END;
    FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO Write(ainv[i,j]:11); 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]*ainv[k,j]
        END;
    Writeln;
    FOR i:=1 TO n DO BEGIN
        FOR j:=1 TO n DO Write(ajed[i,j],'   '); Writeln
        END;
    Readln;

END.