Здравствуйте! У меня ошибка в алгоритме решения системы линейных алгебраичских уравнений (СЛАУ). Система решается методом LU-факторизации с выбором главного элемента и перестоновкой строк. Ошибка есть потому, что величина вектора невязки (т. е. вектора, в котором записана разница найденных и подставленных в исходную систему А ответов Х и вектора ответов системы В, предназначен для проверки точности) слишком большая (около 155). Быть может кто-то сталкивался, помогите, пожалуйста, найти "жучка".

uses crt;
type
vector=array [1..4] of real;
matryca=array [1..4,1..4] of real;
var
B, Z, Y, X, Nev, Temp, BNB: vector;
A: matryca;
i,j,k: integer;
MP: array[1..4] of integer; {massyv perestanovok}
procedure initmatr(var A:matryca);  {zadannya po4atkovoji matryci}
begin
     A[1,1]:=6;
     A[2,1]:=3;
     A[3,1]:=5;
     A[4,1]:=3;
     A[1,2]:=8;
     A[2,2]:=0;
     A[3,2]:=9;
     A[4,2]:=7;
     A[1,3]:=7;
     A[2,3]:=900;
     A[3,3]:=8;
     A[4,3]:=6;
     A[1,4]:=54;
     A[2,4]:=9;
     A[3,4]:=5;
     A[4,4]:=-1;
end;

procedure initvect(var B:vector); {zadnnya vector-stovpcja vilnyh 4leniv}
begin
     B[1]:=3;
     B[2]:=0;
     B[3]:=7;
     B[4]:=5;
end;
procedure initvecty(var Y:vector);
begin
Y[1]:=0;
Y[2]:=0;
Y[3]:=0;
Y[4]:=0;
end;

procedure initmass(var MR:array of integer);
begin
     MR[1]:=1;
     MR[2]:=2;
     MR[3]:=3;
     MR[4]:=4;
end;

procedure printmatr(var A:matryca); {vyvid matryci}
begin
for i:=1 to 4 do
    begin
    for j:=1 to 4 do
    write(A[i,j]:7:1);
    writeln;
    end;
    readkey;
end;

procedure peretvorennya(var A:matryca); {peremishennya golovnuh elementiv}
var
maxel:real;{maksymalnyj element stroki}
promel:real;{promighyi element - slugbova zminna}
promelint:integer;
l:integer;
begin
     for k:=1 to 4 do
     begin
          for j:=k to 4 do
          begin
               maxel:=0;
               for i:=k to 4 do
               if (abs(A[i,j])>maxel) then
               begin
                    maxel:=A[i,j];
                    if i<>k then
                    begin
                         promelint:=MP[i];
                         MP[i]:=MP[k];
                         MP[k]:=promelint;
                         promel:=B[i];
                         B[i]:=B[k];
                         B[k]:=promel;
                         for l:=1 to 4 do
                         begin
                             promel:=A[l,j];
                             A[l,j]:=A[k,j];
                             A[k,j]:=promel;
                         end;
                    end;

               end;
           end {loop i}
     end; {loop k}
end; {peretvorennya}
procedure factorizacia(var A:matryca);
var R:real;
begin
     for k:=1 to 4 do
     begin
          for i:=k+1 to 4 do
          begin
               R:=A[i,k]/A[k,k];
               A[i,k]:=R;
               for j:=k+1 to 4 do
               begin
               A[i,j]:=A[i,j]-A[k,j]*R;
               end; {loop j}
          end;{loop i}
     end;{loop k}
end;{factorizacia}
procedure rozvL(var A:matryca;var  Y,B:vector);
var S:real;
begin
     for i:=1 to 4 do
     begin
     Temp[i]:=A[i,i];
     A[i,i]:=1;
     end;
     Y[1]:=B[1]/A[1,1];
     for i:=2 to 4 do
     begin

          S:=0;
          for j:=1 to i-1 do
          S:=S+A[i,j]*Y[j];
          Y[i]:=(B[i]-S)/A[i,i];
     end; {loop i}
end; {rozvL}
procedure rozvU(var A:matryca;var  Z,Y:vector);
var S:real;
begin
     for i:=1 to 4 do
     A[i,i]:=Temp[i];
     Z[4]:=Y[4]/A[4,4];
     for i:=4-1 downto 1 do
     begin
          S:=0;
          for j:=i+1 to 4 do
          S:=S+A[i,j]*Z[j];
          Z[i]:=(Y[i]-S)/A[i,i];
     end;{loop i}
end;{rozvU}
procedure perestanovka(var X,Z:vector);
begin
for i:=1 to 4 do
X[MP[i]]:=Z[i];
end;
procedure printvect(var X:vector);
begin
for i:=1 to 4 do
begin
writeln;
write('  ','X[',i,']=',X[i]:3:20);
end;
readkey
end;
procedure nevyazka(var BNB:vector; var A:matryca);
var S:real;
begin
     initmatr(A);
     for i:=1 to 4 do
     begin
     S:=0;
     for j:=1 to 4 do
         S:=S+A[i,j]*X[j];
     BNB[i]:=S-B[i];
     end;
end;
{==================================MAIN======================================}
begin
clrscr;
initmatr(A);
initvect(B);
initvecty(Y);
initmass(MP);
textcolor(white);
textbackground(lightgray);
writeln('                                                                        Matryca:');
textcolor(lightgray);
textbackground(black);
printmatr(A);
peretvorennya(A);
textcolor(white);
textbackground(lightgray);
writeln('                                                            Peretvorena matryca:');
textcolor(lightgray);
textbackground(black);
printmatr(A);
factorizacia(A);
writeln;
rozvL(A,Y,B);
rozvU(A,Z,Y);
perestanovka(X,Z);
textcolor(white);
textbackground(lightgray);
writeln('                                                              Vector rozvjazkiv:');
textcolor(lightgray);
textbackground(black);
printvect(X);
writeln;
nevyazka(BNB,A);
textcolor(white);
textbackground(lightgray);
writeln('                                                                Vector nevjazok:');
textcolor(lightgray);
textbackground(black);
printvect(BNB);
end.

А зачем приводишь код? Лучше бы описал подробнее сам алгоритм.

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