Решение СЛАУ матодом LU-разложения.

[gotex]

Вобщем то есть прога, но она вылетает при решении матрицы 7*7, 8*8..(может еще в каких то случаях) :(

может кто подскажет в чем косяк:

вот вводимые данные:
7(размерность)
1 3 2 0 1 1 4 (матрица)
1 3 1 2 2 4 6
1 3 5 3 3 1 9
3 2 8 7 6 5 2
2 5 7 6 6 5 3
1 2 2 4 6 3 2
2 4 5 9 7 6 3

3 5 1 9 8 1 3(столбец свободных членов)

Код:

program LU_Main;

{$APPTYPE CONSOLE}

uses
SysUtils;

type
matrica=array [1..15,1..15] of real; //тип для матрицы коэффициэнтов системы
mas=array [1..15] of real; //тип для матриц свободных членов и решений

var
matr:text; //файлы из которых считываются матрицы
i,j:byte; //счетчики
flag:boolean; //флаг успешного выполнения разбиения
n:byte; //размерность матриц
A:matrica; //матрица коэффициентов
U:matrica; //матрица коэффициентов и матрица U LU-разложения
D:matrica; //матрица с замененым столбцом
F:mas; //матрица свободных членов
X:mas; //матрица решений
det_main,det_dop:real; //определители для оригинальной матрицы A, и для матрицы с заменой столбца на свободный член соответственно


//функция для расчета дополнительной суммы произведений элементов матриц
//v - перменная показатель варианта нахождения суммы
function sum(const i,j:byte;const L,U:matrica;const v:byte):real;
var
k:byte; //счетчик
dop:real; //дополнительная переменная для ресчета суммы
begin
dop:=0;

for k:=1 to i-1 do
case v of
1:dop:=dop+L[i,k]*U[k,j];
2:dop:=dop+L[j,k]*U[k,i];
end;
sum:=dop;
end;

//функция раскладывающая матрицу А(размерности n) LU-разложением, в матрицу LU
//в случае успешного выполнения возвращает True
function LUrazlozhenie(A:matrica;const n:byte;var LU:matrica):boolean;
var
i_dop:integer;
dop:real; //дополнительная
i,j: byte; //счетчики
L,U:matrica; //L и U матрицы для LU-разложения
begin
for i:=1 to n do
for j:=1 to n do
U[i,j]:=0;


i_dop:=0;
if A[1,1]=0 then
begin
for i:=2 to n do
begin
if A[i,1]<>0 then
i_dop:=i;
break;
end;
if i_dop=0 then
begin
Writeln('Error!');
exit;
end; //при первом члене равном нулю
for j:=1 to n do //осуществлям перестановку строк
begin
dop:=A[i_dop,j];
A[i_dop,j]:=A[1,j];
A[1,j]:=dop;
end;
end;

for j:=1 to n do {**********************}
U[1,j]:=A[1,j]; { }
{ }
for j:=2 to n do { }
L[j,1]:=A[j,1]/U[1,1]; { } { }
{ Выполняем алгоритм }
for i:=2 to n do { LU-разложения }
begin { }
for j:=i to n do { }
U[i,j]:=A[i,j]-sum(i,j,L,U,1); { }
{ }
for j:=i+1 to n do { }
L[j,i]:=(A[j,i]-sum(i,j,L,U,2))/U[i,i]; { }
end; { }
{ }
for i:=1 to n do { }
for j:=1 to n do { }
LU[i,j]:=U[i,j]; { }
{ }
LUrazlozhenie:=True; {**********************}

Writeln;
for i:=1 to n do
begin
for j:=1 to n do
Write(U[i,j]:6:2,' ');
Writeln;
end;

end;

[gotex]

//функция нахождения определителя LU-разложенной матрицы
function det(const LU:matrica;const n:byte):real;
var
i:byte; //счетчик
dop:real; //дополнительная переменная для расчетов определителя
begin
dop:=1; {********************************** ***************}
{ }
for i:=1 to n do {находим определитель, как произведение элементов }
dop:=dop*LU[i,i]; {главной диагонали LU-разложенной матрицы }
{ }
det:=dop; {********************************** ***************}
end;

//процедура производящая замену столбца основной матрицы на столбец свободных коэффициентов
procedure Zamena(const A:matrica;const F:mas;const n:byte;var D:matrica;const k:byte);
var
i,j:byte; //счетчики
begin
for i:=1 to n do
for j:=1 to n do
D[i,j]:=A[i,j];

for i:=1 to n do
D[i,k]:=F[i];
end;

begin
flag:=False;

for i:=1 to n do
for j:=1 to n do
U[i,j]:=0;

Assign(matr,'matr3.txt');
Reset(matr);

Readln(matr,n); //ввод размерности

//ввод значений матрицы коэффициентов
for i:=1 to n do
begin
for j:=1 to n do
Read(matr,A[i,j]);
Readln(matr);
end;

Readln(matr);

//вывод матрицы коэффициентов на экран
Writeln('A:');
for i:=1 to n do
begin
for j:=1 to n do
Write(A[i,j]:6:2,' ');
Writeln;
end;

//ввод матрицы свободных членов
for i:=1 to n do
Read(matr,F[i]);
Readln(matr);

//вывод матрицы свободных членов
Writeln('F:');
for i:=1 to n do
Write(F[i]:6:2,' ');
Writeln;

flag:=LUrazlozhenie(A,n,U); //находим LU-разложенную матрицу, для матрицы A, поднимаем или опускаем флаг в зависимости от результата выполнения

if flag then //если LU-разбиение успешно
begin
det_main:=det(U,n); //находим определитель оригинальной мартицы


for i:=1 to n do //для каждого решения
begin
Zamena(A,F,n,D,i); //создаем матрицу с замененным столбцом
flag:=LUrazlozhenie(D,n,U); //находим U матрицу LU-разложения для этой матрицы
if flag then //если LU-разбиение успешно
begin
det_dop:=det(U,n); //вычисляем ее определитель
X[i]:=det_dop/det_main; //вычисляем решение
end
else //иначе выходим из цикла
break;
end;

if flag then
Writeln('X:');
for i:=1 to n do
Writeln('X[',i,']=',X[i]:6:2,' '); //выводим полученое значение на экран
end;

Close(matr);
Readln;
readln;
end.