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;

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

[IUnknown]

У тебя как минимум в двух местах неправильно сделан перенос с С++... Оба - циклы по K. И там и там надо ходить до (i - 1), а не до i... Тестировал программу на той же матрице, что и Сишный вариант?

P.S.

Матрицы вот такие должны получаться (это результат работы твоего чуть-чуть подпиленного дельфийского кода):

L-matrix:
    1.0000    0.0000    0.0000    0.0000
   -8.0000    1.0000    0.0000    0.0000
    0.0000    3.0000    1.0000    0.0000
    4.0000    5.1250  -10.4688    1.0000
U-matrix:
    2.0000   -5.0000    6.0000    1.0000
    0.0000    8.0000   -7.0000    3.0000
    0.0000    0.0000   -4.0000    1.0000
    0.0000    0.0000    0.0000   -4.9063

?

Да, кстати, у тебя там балаган с возвращаемыми значениями, ты под Res память-то не выделил. Да и как собрался возвращать и L- и U-матрицы? Проще передать через параметры...

Сообщение отредактировано: IUnknown - 1.11.2011 19:20