{$N+}{$E+}
Program Demzpr;

USES Crt;

CONST
   np = 70;
TYPE
   Real = Single;
   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;

PROCEDURE mprove(VAR a,alud: RealArrayNPbyNP;
                          n: integer;
                   VAR indx: IntegerArrayNP;
                      VAR b: RealArrayNP;
                      VAR x: RealArrayNP);
{ zlepsuje reseni X systemu N rovnic, na ktere byla pouzita LU dekompozice
  a je matice soustavy, alud je LU matice vznikla pomoci LUDCMP
                        indx jsou indexy radek pri LU dekompozici
  b je prava strana rovnic
  procedura muze byt opakovane pouzita pro iteraci
}
VAR
   j,i: integer;
   sdp: double;
   r: ^RealArrayNP;
BEGIN
   new(r);
   FOR i := 1 TO n DO BEGIN
      sdp := -b[i];
      FOR j := 1 TO n DO
         sdp := sdp+a[i,j]*x[j];
      r^[i] := sdp
   END;
   lubksb(alud,n,indx,r^);
   FOR i := 1 TO n DO
      x[i] := x[i]-r^[i];
   dispose(r)
END;

VAR a,apom: RealArrayNPbyNP;
    n,i,j,k: integer;
    indx: IntegerArrayNP;
    b,bpom,bvys,bold: RealArrayNP;
    detzn,dethod,eps,delt,delmax,bnorm,bpmin,bpmax,dbp:Real;
    fvst,fvys:Text;
    vstnam:String;
Begin
    Writeln(' Zadej dimenzi matice, presnost reseni eps ');
    Readln(n,eps);
    Writeln(' Zadej x(1),x(n)');
    Readln(bpmin,bpmax);
    dbp := (bpmax-bpmin)/(n-1);
    FOR i:=1 TO n DO bpom[i]:= bpmin + (i-1)*dbp;

    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)');
    Readln(bpmin,bpmax);
    dbp := (bpmax-bpmin)/(n-1);
    FOR i:=1 TO n DO BEGIN b[i]:=bpmin+(i-1)*dbp; 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]);
    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;
    REPEAT
       FOR i:=1 TO n DO bold[i]:=b[i];
       mprove(apom,a,n,indx,bpom,b);
       delmax:=0.;
       bnorm:=0;
       FOR i:=1 TO n DO BEGIN
          IF (ABS(b[i]-bold[i]) > delmax ) THEN delmax := ABS(b[i]-bold[i]);
          IF (ABS(b[i]) > bnorm ) THEN bnorm := ABS(b[i]);
          END;
       delt := delmax/bnorm;
       FOR i:=1 TO n DO Write(b[i]);
            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;
       Writeln;
       Writeln(' DELTA = ', delt);
       Readln;
    UNTIL delt <= eps;
    Writeln('Koncim');Readln
End.