Eliminacja Seidla

program Eliminacja_Seidla;
uses
   crt;
const
   N = 3;
type
   vector = array[1..N] of real;
   matrix = array[1..N,1..N] of real;
const
   a: matrix = ((6.10, 2.20, 1.20),
                (2.20, 5.50,-1.50),
                (1.20,-1.50, 7.20));

   b: vector = (16.55, 10.55, 16.80);
var
   x, blad: vector;

procedure wynik(n, q: integer;
                x, blad: vector);
   var
      i: integer;
   begin
      writeln('  -- ',q,' --               blad');
      for i:= 1 to n do
         writeln('x',i,'= ',x[i]:13:10,'  ',blad[i]:13:10);
      readkey;
   end;

procedure saidl(n: integer; a: matrix; b: vector;
                var x, blad: vector);
   var
       d, xp: vector;
      i, j, q, odp: integer;
      s: real;
   begin
      for i:= 1 to n do
         begin
            d[i]:= -1/a[i,i];{ macierz diagonalna -D^-1 }
            a[i,i]:= 0;      { macierz bez przekatnej = L+U }
            x[i]:= 0;        { vector Xi }
            xp[i]:= 0;       { vector Xi-1 }
         end;
      write('podaj ilosc iteracji: ');
      readln(odp);
      q:= 0;                 { liczba iteracji }
      repeat
         for i:= 1 to n do
            begin
               s:= 0;
               for j:= 1 to n do
                  begin                    { dla zmodyfikowanej metody }
                     if i<j then           { gdy nad przekatna }
                        s:= s+a[i,j]*xp[j];{ mnozenie mac U*Xi-1 }
                     if i>j then           { gdy pod przekatna }
                        s:= s+a[i,j]*x[j]; { mnozenie mac L*Xi }
                  end;
               x[i]:= d[i]*(s-b[i]);       { Xi:= -D^-1*(L*Xi+U*Xi-1-B) }
               blad[i]:= abs(x[i]-xp[i]);
               xp[i]:= x[i];
            end;
         wynik(n,q,x,blad);
         inc(q);
      until q=odp+1;
   end;

begin
   clrscr;
   saidl(n,a,b,x,blad);
end.