// Metoda Gaussa-Seidela
// www.algorytm.org
// (c) 2006 Tomasz Lubinski

program gs;
{$APPTYPE CONSOLE}

uses
  SysUtils;

var
A: Array [1..100] of Array [1..100] of Double;
L: Array [1..100] of Array [1..100] of Double;
D: Array [1..100] of Array [1..100] of Double;
U: Array [1..100] of Array [1..100] of Double;
b: Array [1..100] of Double;
x: Array [1..100] of Double;

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

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

   // Get values of A
   for i:=1 to n do
      for j:=1 to n 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 n do
      begin
         writeln('b[' + IntToStr(i) + '] = ');
         read(b[i]);
      end;

   // Divide A into L + D + U
   for i:=1 to n do
      for j:=1 to n do
         begin
            if (i < j) then
               U[i][j] := A[i][j]
            else if (i > j) then
               L[i][j] := A[i][j]
            else
               D[i][j] := A[i][j];
         end;

   // Calculate D^-1
   for i:=1 to n do
      D[i][i] := 1/D[i][i];

   // Calculate D^-1 * b
   for i:=1 to n do
      b[i] := D[i][i] * b[i];

   //Calculate D^-1 * L
   for i:=1 to n do
      for j:=1 to i-1 do
         L[i][j] := D[i][i] * L[i][j];

   //Calculate D^-1 * U
   for i:=1 to n do
      for j:=i+1 to n do
         U[i][j] := D[i][i] * U[i][j];

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

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

   for k:=1 to iter do
      for i:=1 to n do
         begin
            x[i] := b[i];                       // x = D^-1*b -
            for j:=1 to i-1 do
               x[i] := x[i] - L[i][j]*x[j];    // D^-1*L * x -
            for j:=i+1 to n do
               x[i] := x[i] - U[i][j]*x[j];    // D^-1*U * x
         end;

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

   readln;

end.