LU разложение. Опции

[Krjuger]

Собственно есть задача, полную ее суть писать не стану,но проблема заключается в следующем нужно с помощью LU разложения решить систему.Рабочий код на С++ у меня был и я его попытался перенести на Delphi.Но к несчастью моя затея обернулась неудачей.

type
  TVec = array of Extended; // векторный тип
  TMatr = array of TVec; // матричный тип

function  LUrazl(matr: Tmatr): Tvec;
var
  i, j,k, row, col: Word;
  res: TVec;
  _matr_: TMatr;
  L,U: TMatr;
  Sum: double;
begin
  _matr_ := CopyMatr(matr); // Copy не работает как надо
	row := Length(_matr_);
  col := Length(_matr_[0]);
  SetLength(L,row,col);
  SetLength(U,row,col);
	if row <> col-1 then
    ShowMessage('Ошибка: некорректные размерности матрицы! (МГ)'); //error
  for i:=0 to  row-1 do
    for j:= 0  to col-1 do
   begin
      L [i][j] := 0;
      U [i][j] := 0;

      if i = j then
        L [i][j] := 1;
   end;



//==============================================

//находим первый столбец L[][] и первую строку U[][]

  for i:= 0 to row-1 do
  begin
    L [i][0]:= _matr_ [i][0]/ _matr_ [0][0];
    U [0][i]:= _matr_ [0][i] / L [0][0];
  end;

    for i:= 0 to row-1 do
     begin
          for j:= 0 to col-1 do
            begin
                   U [0][ i] := _matr_[0][ i];
                   L [i][ 0] := _matr_[i][ 0] / U[0][ 0];
                   sum := 0;
                   for k:= 0 to i  do
                   begin
                       sum :=sum+ L[i][ k] * U[k][ j];
                   end;
                   U[i][ j] := _matr_[i][ j] - sum;
                   if (i > j) then
                     L[j][ i] := 0
                   else
                    begin
                       sum := 0;
                       for  k:= 0 to i do
                       begin
                           sum:=sum+L[j][ k] * U[k][ i];
                       end;
                       L[j][ i] := (_matr_[j][ i] - sum) / U[i][ i];//В этой строке выдает деление на ноль так же в дебаге видно,что значение i=65000+
                    end;

            end;
          res[i]:= U[i][ i];
     end;
end;

Непонятна причина такого поведения,ведь на С++ все это работало.Если надо могу добавить файл целиком.

[Krjuger]

Вот С++,так как оно мне требовалось, чтобы посчитать стопочку типовых расчетов,особого энтузиазма в оформлении нет.

#include <iostream>
using namespace std;

int main ()
{
  const int n=4;
  double sum = 0;
  double A [ n ][ n ]= { 
{ 2, -5,6,1},
{ -16, 48,-55,-5},
{ 0, 24,-25,10},
{8,21,30,4 }};
  double L [ n ][ n ];
  double U [ n ][ n ];

//задаем матрицу A[][] ...

  for (int i  = 0; i < n; i++)
  {
   for (int j = 0; j < n; j++)
   {
 //     cout << "\na[" << i+1 << "][" << j+1 << "] = ";
   //   cin >> A [i][j];
     L [i][j] = 0;
     U [i][j] = 0;

      if (i == j)
        L [i][j] = 1;
   }
  }

//==============================================

//находим первый столбец L[][] и первую строку U[][]

  for (int i = 0; i < n; i++)
  {
    L [i][0] = A [i][0]/ A [0][0];
    U [0][i] = A [0][i] / L [0][0];
  }

for (int i = 0; i < n; i++)
           {
               for (int j = 0; j < n; j++)
               {
                   U [0][ i] = A[0][ i];
                   L [i][ 0] = A[i][ 0] / U[0][ 0];
                   double sum = 0;
                   for (int k = 0; k < i; k++)
                   {
                       sum += L[i][ k] * U[k][ j];
                   }
                   U[i][ j] = A[i][ j] - sum;
                   if (i > j)
                   {
                       L[j][ i] = 0;
                   }
                   else
                   {
                       sum = 0;
                       for (int k = 0; k < i; k++)
                       {
                           sum += L[j][ k] * U[k][ i];
                       }
                       L[j][ i] = (A[j][ i] - sum) / U[i][ i];
                   }
               }
           }


//====================================================
   cout << "\n\n";
   for (int i = 0; i < n; i++)
   {
     for (int j = 0; j < n; j++)
	 {	
		 if (i == j)
          L [i][j] = 1;
       cout << "  " << L [i][j] << "   ";
	 }
     cout << "\n\n";
   }
   

   cout << "\n\n";

   for (int i  = 0; i < n; i++)
   {
     for (int j = 0; j < n; j++)
       cout << "  " << U [i][j] << "   ";
     cout << "\n\n";
   }

   return 0;
}

Так как типовые были сдаты, и пара вариантов были посчитаны руками и результаты были одинаковыми,поэтому я счел,что он работает верно.

Да и забыл упомянуть для текущей задачи мне требуеться,чтобы оно считалось для максимальных размерностей 101х100.А то у меня уже был случай,когда Гаусс верно считал для любой матрицы меньше 20,а для больших переставал работать,пришлось полностью пересмотреть концепцию)))

Сообщение отредактировано: Krjuger - 1.11.2011 18:18