// Metoda Jacobiego
// www.algorytm.org
// (c) 2006 Tomasz Lubinski

program jacobi;
{$APPTYPE CONSOLE}

uses
  SysUtils;

var
A: Array [1..100] of Array [1..100] of Double;
M: Array [1..100] of Array [1..100] of Double;
N: Array [1..100] of Double;
b: Array [1..100] of Double;
x1: Array [1..100] of Double;
x2: Array [1..100] of Double;

num, iter: Integer;
i, j, k : Integer;

begin
   // Get n
   Writeln('Metoda Jacobiego');
   Writeln('Rozwiazywanie ukladu n-rownan z n-niewiadomymi Ax=b');
   Writeln('Podaj n');
   Readln(num);
   if (num < 1) or (num > 100) then
      begin
         Writeln('Nieprawidlowa warosc parametru n');
         exit;
      end;

   // Get values of A
   for i:=1 to num do
      for j:=1 to num do
         begin
            writeln('A[' + IntToStr(i) + '][' + IntToStr(j) +'] = ');
            read(A[i][j]);
            if (i = j) and (A[i][j] = 0)  then
               begin
                 writeln('Wartosci na przekatnej musza byc rozne od 0');
                 exit;
               end;
         end;

   // Get values of b
   for i:=1 to num do
      begin
         writeln('b[' + IntToStr(i) + '] = ');
         read(b[i]);
      end;

   // Calculate N = D^-1
   for i:=1 to num do
      N[i] := 1/A[i][i];

   // Calculate M = -D^-1 (L + U)
   for i:=1 to num do
      for j:=1 to num do
         if (i = j) then
            M[i][j] := 0
         else
            M[i][j] := - (A[i][j] * N[i]);

   //Initialize x
   for i:=1 to num do
      x1[i] := 0;

   writeln('Ile iteracji algorytmu wykonac? ');
   readln(iter);

   for k:=1 to iter do
      begin
         for i:=1 to num do
            begin
               x2[i] := N[i]*b[i];
               for j:=1 to num do
                  x2[i] := x2[i] + M[i][j]*x1[j];
            end;
         for i:=1 to num do
            x1[i] := x2[i];
      end;

   writeln('Wynik');
   for i:=1 to num do
      writeln('x[' + IntToStr(i) + '] = ' + FloatToStr(x1[i]));

   readln;

end.