Реализация метода простой итерации [Delphi]

[Elarium]

Ребят, очень нужна ваша помощь! Вообщем я реализую такой вот метод простых итераций:
$x_{n+1} = (E-\tau A)x_n + \tau f$
Чтобы найти $\tau$ нужно найти собственные числа матрицы $A$, важно чтобы матрица была положительно определена, тогда $\tau$ лежит в интервале $0< \tau <2/ \lambda_{max}$, я преобразую матрицу$A$ к треугольному виду и получаю, что элементы на диагонали треугольной матрицы являются её собственными числами. Затем я нахожу $\tau_{opt}=2/(\lambda_{min}+\lambda_{max})$ и запускаю итерационный процесс.

Проблема в том, что repeat при итерациях зацикливается, а это скорее всего из-за того, что процесс не сходится. Может быть некоторые вещи мной были не так поняты или я допустил ошибку в самом написании программы.

Используется синтаксис Delphi

unit Unit1;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, Grids;

type
  mas1 = array[0..3,0..3] of real;
  mas2 = array[0..3] of real;
  TForm1 = class(TForm)
    StringGrid1: TStringGrid;
    Button1: TButton;
    Button2: TButton;
    procedure Button2Click(Sender: TObject);
    procedure Button1Click(Sender: TObject);
  private
    { Private declarations }
  public
    { Public declarations }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

procedure TForm1.Button2Click(Sender: TObject);
begin
close;
end;

procedure iteration(var A:mas1; B:mas2; tau:real);
var
E:mas1;
X,Xplus:mas2;
i,j,z,n:integer;
begin
for i:=0 to 3 do
X[i]:=tau*B[i];
for i:=0 to 3 do
  for j:=0 to 3 do
    begin
    if (j=i) then E[i,j]:=1;
    if (i<>j) then E[i,j]:=0;
    end;
repeat
z:=0;
for i:=0 to 3 do
  begin
    for j:=0 to 3 do
        Xplus[i]:=E[i,j]-tau*A[i,j]+tau*B[i];
    if (abs(Xplus[i]-X[i])<0.001) then
      z:=z+1;
  end;
for i:=0 to 3 do
X[i]:=Xplus[i];
n:=n+1;
until (z=3);
showmessage(floattostr(x[1])+' '+floattostr(x[2])+' '+floattostr(x[3]));

end;

procedure convergence(var A:mas1; B,EV:mas2);
var
sup,inf,tau:real;
i,j:integer;
begin
sup:=EV[0];
for i:=1 to 3 do
  if (EV[i]>sup) then sup:=EV[i];

inf:=EV[0];
for i:=1 to 3 do
  if (EV[i]<inf) then inf:=EV[i];

tau:=2/(sup+inf);

iteration(A,B,tau);
end;

procedure TrMatrix(var A:mas1; B:mas2);
var
M:mas1;
EV:mas2;
k,x,s,i,j:integer;
temp,N:real;
begin
for i:=0 to 3 do
  for j:=0 to 3 do
    M[i,j]:=A[i,j];

 for k:=0 to 2 do
 begin
{Смена строк, если M[k,k] = 0}
  s:=k;
  while ((M[k,k]=0) and (s<3)) do
     begin
       s:=s+1;
       if (abs(M[s,k])>0) then
         for j:=0 to 3 do
           begin
             temp:=M[s,j];
             M[s,j]:=a[k,j];
             M[k,j]:=temp;
           end;
     end;

  if (abs(M[k,k])>0) then
 {Приведение к треугольному виду}
   begin
    for i:=k+1 to 2 do
     begin
      N:=M[i,k]/M[k,k];
      for j:=k to 3 do
        begin
         M[i,j]:=M[i,j]-M[k,j]*N;
         if (abs(M[i,j])<0.0000001) then
          M[i,j]:=0;
        end;
     end;
   end;
 end;
  x:=0;
  for i:=0 to 3 do if (a[i,i]>0) then x:=x+1;

if (x<4) then showmessage('Матрица А должна быть строго положительно определена');
if (x=4) then
  begin
    for i:=0 to 3 do
    EV[i]:=M[i,i];
    convergence(A,B,EV);
  end;


end;

procedure TForm1.Button1Click(Sender: TObject);
var
A:mas1;
B:mas2;
i,j:integer;
begin
  for i:=0 to 3 do
    for j:=0 to 3 do
      A[i,j]:=strtofloat(stringgrid1.cells[j,i]);

  for i:=0 to 3 do
    B[i]:=strtofloat(stringgrid1.cells[4,i]);

TrMatrix(A,B);


end;

end.
 

[Imperator]

Некоторые советы:
0. Запросите пошаговую работу программы и посмотрите, где что не так.
1. Выберите несколько различных $\tau$ (скажем, $\tau_1=\frac{2}{\lambda_{max}+1}$).
2. Проверьте отдельно на конкретных примерах процедуру приведения матрицы к треугольному виду.

[Elarium]

Блин, у меня условие в until'e было не корректно поставлено :)

[Elarium]

Ещё была ошибка в самом ходе итерации, я там бред написал. Сейчас всё работает отлично, всем спасибо кто хотел помочь.