Динамические матрицы в процедуре, Есть проблема с динамическими данными... Опции

[P-Tigr]

Столкнулся с одной проблемой, 2 дня уже бьюсь... что делать - не знаю...
Если коротко, то дело вот в чем.
Я написал процедуру для подсчета определителя матрицы любого порядка (методом Гаусса). И там над матрицей выполняются различные преобразования (складывания строк и т.д.).
Причем - внимание - по условию задания матрица должна быть обязательно динамической!
Процедура работает нормально, определитель находит, НО в то же время изменяет исходную матрицу, а этого нельзя допустить!! Происходит это, как мне кажется, из-за того, что такие матрицы - ссылки на память, поэтому передать ее в процедуру строго по значению нельзя...

Замучился я с этими динамическими структурами...
Please, умные люди, help me!

[volvo]

P-Tigr, давайте так: вы можете показать код процедуры? тогда покажите, как именно вы работаете с матрицей, и можно будет сказать, как избежать ее изменения. Теоретически (не зная деталей реализации) очень трудно давать советы...

[P-Tigr]

Вот выдержки из кода:

type
  TDMas = array of real;
  TDMatr = array of TDMas;
............

{<--------------------------------------------->}
{ф-ция ищет определитель}

function Opred_Gauss(A:TDMAtr):real;
var i,j,jk,ik,new:integer;
    tmp,mnoj,mn1,mn2:real;
    found:boolean;
    AC:TDMAtr;
begin
  AC:=A;
  A:=nil;
  // проверили равность определителя нулю
  for j:=0 to high(AC) do
    if AC[j,j]=0 then
      begin
        found:=false;
        For new:=0 to high(AC) do
          if AC[new,j]<>0 then
            begin
              for jk:=0 to high(AC[new]) do
                begin
                  tmp:=AC[j,jk];
                  AC[j,jk]:=AC[new,jk];
                  AC[new,jk]:=tmp;
                end;
              found:=true;
              break;
            end;
        if not found then
          begin
            Result:=0;
            Exit;
          End;
      end;
  // приволим матрицу у верхне-треугольному виду
  {(все эл-ты под главной диагональю сделаем = 0)}
  for j:=0 to high(AC)-1 do
    begin
      for i:=j+1 to high(AC) do
        begin
          mn1:=AC[i,j]; mn2:=AC[j,j];
          mnoj:=mn1/mn2;
          for jk:=0 to high(AC) do
            begin
              AC[i,jk]:=AC[i,jk]-AC[j,jk]*mnoj;
            end;
        end;
    end;
  Result:=1;
  For i:=0 to high(AC) do
    Result:=Result*aC[i,i];
end;

{<--------------------------------------------->}
{процедура обработки: поместить определитель в глав. диагональ и 
потом отразить все эл-ты симметрич. относ. глав. диагонали}

procedure Work_Matr(a:TDMatr; var b:TDMatr);
var opred,tmp:real;
    i,j,tmpi:integer;
begin
  b:=nil;
  SetLength(b,length(a),length(a));
  for i:=0 to high(a) do
      for j:=0 to high(A) do
        begin
          tmp:=a[i,j];
          b[i,j]:=tmp;
        end;

  opred:=Opred_Gauss(a);
  b[0,0]:=opred;
  {for i:=1 to high( B ) do
    begin
      for j:=0 to i-1 do
        begin
          tmp:=b[i,j];
          tmpi:=i-j;
          b[i,j]:=b[i-tmpi,j+tmpi];
          b[i-tmpi,j+tmpi]:=tmp;
          b[i,i]:=opred;
        end;
    end;}
end;

{<--------------------------------------------->} 

вроде все... здесь А - исходная матрица
В - полученная

[P-Tigr]

Причем кроме вышеописанной проблемы, при данных допустим
1 2 3 4
1 2 3 4
1 2 3 4
1 2 3 4

генерируется странная ошибка - "Floating operation error"

[volvo]

Компилятор какой? Delphi или FPC? Турбо-Паскаль открытые массивы не поддерживает.

[P-Tigr]

Естессно, Delphi ver 6

[volvo]

Естественно, что я переношу тему в раздел Дельфи

[P-Tigr]

:D Just OK. I'm sorry, первый раз на форуме...

[P-Tigr]

Thanks, Volvo!

[volvo]

Значит, предложение такое: при входе в функцию вычисления определителя ты не присваиваешь адрес массива A в AC (как это происходит сейчас, а копируешь массив A в AC полностью и работаешь с его копией). Копирование производишь так:

    setlength(ac, length(a));
    for i := 0 to pred(length(a)) do
      begin
        setlength(ac[i], length(a[i]));
        for j := 0 to pred(length(a[i])) do
          ac[i, j] := a[i, j]
      end;

И не забудь убрать строки

 AC:=A;
A:=nil;

иначе все копирование будет лишено смысла
я только что прогнал этот вариант - исходный массив не меняется (правда я не использовал процедуру Work_Matr, а проверил только вычисление определителя...)

[P-Tigr]

Извините пожалуйста, еще 1 вопрос, помогите разобраться!
Почему при вводе некоторых данных, например
1 1 1
1 1 1
1 1 1
в блоке функции нахождения определителя:

for j:=0 to high(AC)-1 do
    begin
      for i:=j+1 to high(AC) do
        begin
          mn1:=AC[i,j]; mn2:=AC[j,j];
// ОШИБКА - в следующей строке   генерируется       
          mnoj:=mn1/mn2;
          for jk:=0 to high(AC) do
            begin
              AC[i,jk]:=AC[i,jk]-AC[j,jk]*mnoj;
            end;
        end;
    end;

в строке
mnoj:=mn1/mn2;

появляется ошибка, связанная с плавающей запятой:
"Invalid floating point operation" или "Floating divizion by zero".
ПОЧЕМУ?

P.S.1. Ошибка появляется и при вводе разных чисел, но правда не всегда - очень странно...
P.S.2. Я присоединил к письму exe-файл в rar-архиве - чтоб удобней было тестировать.

[volvo]

появляется ошибка, связанная с плавающей запятой:
"Invalid floating point operation" или "Floating divizion by zero".
ПОЧЕМУ?

Это не связано с вводимыми данными... Только с тем, что у тебя в отмеченных строках просто напросто происходит выход за границы массива (не забывай, что динамические массивы индексируются с 0 до N-1), и кто знает, какие значения там хранятся...

for j:=0 to high(AC)-1 do
   begin
     for i:=j+1 to high(AC) do  // **** Здесь - i принимает значение на 1 больше чем нужно
       begin
         mn1:=AC[i,j]; mn2:=AC[j,j];

         mnoj:=mn1/mn2; // ОШИБКА

         for jk:=0 to high(AC) do
         // **** Здесь - jk тоже принимает значение на 1 больше, чем допустимо,
         // это тоже может привести к проблемам
           begin
             AC[i,jk]:=AC[i,jk]-AC[j,jk]*mnoj;
           end;
       end;
   end;

Я совсем забыл написать, что я менял границы индексов когда проверял... Привычка...

P.S. Смотри и выше в тексте функции. Там тоже есть выход за пределы массива...

[volvo]

P-Tigr, небольшая поправочка:
Выход за пределы массива тут ни при чем...

Сейчас еще раз прогнал твою программу... Вот посмотри на скриншот, который я сделал в определенный момент. Как ты думаешь, что произойдет на следующем шаге?

Я думаю, что деление на 0...

Это в принципе ясно, ведь если у тебя матрица:
1, 1, 1
1, 1, 1
1, 1, 1
то после вычитания первой строки из второй и из третьей, она становится такой:
1, 1, 1
0, 0, 0
0, 0, 0
и далее как только ты пытаешься найти множитель для вычитания 3-ей строки из второй - Oops "Деление на ноль". Но ведь на самом-то деле его и искать не надо... Матрица уже приведена к треугольному виду... Я думаю, надо добавить проверку, не являются ли первые N символов строки N нулевыми, т.е. не является ли данная строка уже строкой подходящей для треугольной матрицы...

Добавлено:
Ну или вот так (взято из реализации trminator-а отсюда: Определитель матрицы ):

function Opred_Gauss(A:TDMAtr):real;
// Описание переменных...
const
  eps = 1e-10;
begin
  { Это то самое копирование матрицы }
  setlength(ac, length(a));
  for i := 0 to pred(length(a)) do
    begin
      setlength(ac[i], length(a[i]));
      for j := 0 to pred(length(a[i])) do
        ac[i, j] := a[i, j]
    end;

  { А теперь - сам алгоритм }
  for i := 0 to high(AC) do
    begin
      if abs(aс[i,i])<eps then
        begin result :=0.0; exit end;
      for j := succ(i) to high(AC) do
        begin
          mnoj:=a[j,i]/a[i,i];
          for jk:=i to high(AC) do
            a[j,jk]:=a[j,jk]-d*a[i,jk];
        end;
    end;
  result := 1.0;
  for i:=1 to high(AC) do result:=result*a[i,i];
end;

[P-Tigr]

ОК, ошибка "Floating divizion by zero" устранена... Мне надо ставить 2 за внимательность...

Но осталась 2-я ошибка - "Invalid floating point operation".
Итак, код такой:

function Opred_Gauss(A:TDMAtr):real;
// Описание переменных...
var i,j,jk,ik,new:integer;
    tmp,mnoj,mn1,mn2:real;
    found:boolean;
    AC:TDMAtr;
const  eps = 1e-10;
begin
  { Это то самое копирование матрицы }
  setlength(ac, length(a));
  for i := 0 to pred(length(a)) do
    begin
      setlength(ac[i], length(a[i]));
      for j := 0 to pred(length(a[i])) do
        ac[i, j] := a[i, j]
    end;
  { А теперь - сам алгоритм }
  for i := 0 to high(AC) do
    begin
      if abs(ac[i,i])<eps then
        begin
          result :=0.0;
          exit
        end;
      for j := succ(i) to high(AC) do
        begin
          mnoj:=a[j,i]/a[i,i];
          for jk:=i to high(AC) do
            a[j,jk]:=a[j,jk]-mnoj*a[i,jk];
        end;
    end;
  result := 1.0;
  for i:=1 to high(AC) do
  result:=result*a[i,i];
end;

При данных

1 1 1
1 1 1
1 1 1

возникает вышеописанная ошибка в той же строке
mnoj:=a[j,i]/a[i,i];

Немного помучившись и покапавшись в справочниках, нашел, что это бывает связано с проблемами сопроцессора, и устраняется путем добавления строчки

asm FINIT end;

тогда вместо ошибки возвращается константа NAN.

Окончательный вариант кода выглядит так:

function Opred_Gauss(A:TDMAtr):real;
// Описание переменных...
var i,j,jk,ik,new:integer;
    tmp,mnoj,mn1,mn2:real;
    found:boolean;
    AC:TDMAtr;
const  eps = 1e-10;
begin
  asm FINIT end;

  { Это то самое копирование матрицы }
  setlength(ac, length(a));
  for i := 0 to pred(length(a)) do
    begin
      setlength(ac[i], length(a[i]));
      for j := 0 to pred(length(a[i])) do
        ac[i, j] := a[i, j]
    end;
  { А теперь - сам алгоритм }
  for i := 0 to high(AC) do
    begin
      if abs(ac[i,i])<eps then
        begin
          result :=0.0;
          exit
        end;
      for j := succ(i) to high(AC) do
        begin
          mnoj:=a[j,i]/a[i,i];
          for jk:=i to high(AC) do
            a[j,jk]:=a[j,jk]-mnoj*a[i,jk];
        end;
    end;
  result := 1.0;
  for i:=1 to high(AC) do
  result:=result*a[i,i];
  If Result = NAN then
    Result:=0;
end;

Вот теперь (!) все работает без ошибок... УРА!

Но интересно, есть ли другие пути решения проблемы, без асмы? И отчего она все-таки, эта ошибка?

[volvo]

ОК, ошибка "Floating divizion by zero" устранена... Мне надо ставить 2 за внимательность...

Но осталась 2-я ошибка - "Invalid floating point operation".
Итак, код такой:

function Opred_Gauss(A:TDMAtr):real;
// Описание переменных...
begin
  { Это то самое копирование матрицы }
  ...
  { А теперь - сам алгоритм }
  for i := 0 to high(AC) do
    begin
      if abs(ac[i,i])<eps then
        begin
          result :=0.0;
          exit
        end;
      for j := succ(i) to high(AC) do
        begin
          mnoj:=a[j,i]/a[i,i];  // **** !!! ****
          for jk:=i to high(AC) do
            a[j,jk]:=a[j,jk]-mnoj*a[i,jk]; // **** !!! ****
        end;
    end;
  result := 1.0;
  for i:=1 to high(AC) do
  result:=result*a[i,i];
end;

Вам надо еще раз поставить 2 за "внимательность"... Неужели нельзя было проверить все ли A заменены на AC? Ведь NaN появляется при той же самой ошибке - "деление на 0"... Почему? Проверяется на 0 значение AC, а операция производится с A !!!

В общем, окончательный вариант функции без всякого Асма:

Function Opred_Gauss(A:TDMatr):double;
  const eps=1E-10;
  var
    i,j,jk:integer;
    mnoj:real;
    AC:TDMAtr;
  begin
    setlength(ac, length(a));
    for i := 0 to pred(length(a)) do
      begin
        setlength(ac[i], length(a[i]));
        for j := 0 to length(a[i]) do
          ac[i, j] := a[i, j]
      end;

  for i := 0 to high(AC) do
   begin
     if abs(ac[i,i])<eps then
       begin
         result :=0.0;
         exit
       end;
     for j := succ(i) to high(AC) do
       begin
         mnoj:=ac[j,i]/ac[i,i];
         for jk:=i to high(AC) do
           ac[j,jk]:=ac[j,jk]-mnoj*ac[i,jk];
       end;
   end;
    Result:=1;
    For i:=0 to high(AC) do
      Result:=Result*aC[i,i];
  end;

[P-Tigr]

Блин...
Спасибо за помощь! Буду учиться

[$ad!st]

вопрос не в тему...
а что такое матрица???
вобщем объясните...
спасибо

[klem4]

http://forum.pascalnet.ru/index.php?showtopic=2694