2014 dxdy logo

Научный форум dxdy

Математика, Физика, Computer Science, Machine Learning, LaTeX, Механика и Техника, Химия,
Биология и Медицина, Экономика и Финансовая Математика, Гуманитарные науки




На страницу 1, 2  След.
 
 Маленький вопрос по С++. Детерминант.
Сообщение22.02.2006, 16:55 
Господа!

Вот скажите как лучше написать прогу по вычислению детерминанта матрицы.
Задачка школьная, но всё таки. :oops:
Пробовал рекурсией, но замучился с динамическими массивами.
Может кто-то знает алгоритм по-проще?
Я ещё прикидывал. Для того, чтобы сосчитать дет. матрицы 10*10 надо создать где-то 2,6 млн. массивов и выделить около 77Мб оперативы :shock:. Может я иду не тем путём? :?:

 
 
 
 Re: Маленький вопрос по С++. Детерминант.
Сообщение22.02.2006, 17:20 
Аватара пользователя
yog писал(а):
...
Может кто-то знает алгоритм по-проще?
...

Так с рекурсией и есть самый простой - в плане реализации. Только неэффективный.
А чтобы быстрее посчитать, надо сделать треугольную матрицу и перемножить элементы на главной диагонали.
Цитата:
...
Я ещё прикидывал. Для того, чтобы сосчитать дет. матрицы 10*10 надо создать где-то 2,6 млн. массивов и выделить около 77Мб оперативы :shock:. Может я иду не тем путём? :?:

Это, наверное, если все промежуточные расчёты в памяти хранить, и то, думаю, слишком завышенная оценка. Я делал на BP процедуру вычисления определителя прямым расчётом - она у меня под DOS и для 100x100 матрицы определитель считала, правда, долго.

 
 
 
 
Сообщение22.02.2006, 17:57 
Аватара пользователя
Прошу прощения, 100x100 я, наверное, не [до]считал :)
Однако, всё равно, даже если реализовывать прямой расчёт, то не нужно передавать рекурсивно в подпрограмму матрицы меньшего размера. Достаточно запоминать номера уже использованных строк/столбцов.
А вот моя программа (прямой расчёт):
Код:
const n = 10;
type matrix = array [1..n, 1..n] of real;
type tSP = ^tST;
     tST = record
       str, col, cnt, sign: integer;
       arg: real;
       prev: tSP;
     end;
var SP: tSP;
procedure push (elem: tST);
  var t: tSP;
begin
  t := SP;
  new (SP);
  SP^ := elem;
  SP^.prev := t;
end;
procedure pop (var elem: tST);
  var t: tSP;
begin
  elem := SP^;
  t := SP;
  SP := SP^.prev;
  dispose (t);
end;
function stEmpty: boolean;
begin
  stEmpty := (SP = nil);
end;
function InStackStr (str: integer): boolean;
  var t: tSP;
begin
  t := SP;
  while t <> nil do begin
    if t^.str = str then begin
      InStackStr := true;
      exit;
    end;
    t := t^.prev;
  end;
  InStackStr := false;
end;
function InStackCol (col: integer): boolean;
  var t: tSP;
begin
  t := SP;
  while t <> nil do begin
    if t^.col = col then begin
      InStackCol := true;
      exit;
    end;
    t := t^.prev;
  end;
  InStackCol := false;
end;
var matr: matrix;
    str0: integer;
{ const matr0: matrix = ((1,2,4,6),(8,5,4,9),(3,1,7,2),(4,7,9,2)); }
function definitor: real;
  var level: integer;
      s: tST;
      res: real;
  procedure findCol;
    var i, j: integer;
  begin
    i := s.cnt;
    j := 1;
    repeat
      if not (InStackCol (j)) then
        Dec (i);
      if i = 0 then
        break;
      Inc (j);
    until false;
    s.col := j;
  end;
  procedure findStr;
    var i: integer;
  begin
    i := 1;
    while InStackStr (i) do
      Inc (i);
    s.str := i;
    findCol;
  end;
begin
  level := n;
  res := 0.0;
  with s do begin
    str := str0;
    col := 1;
    cnt := 1;
    sign := 1;
    arg := 1.0;
    repeat
      if level = n then
        arg := matr [str, col] * sign
      else
        arg := arg * matr [str, col] * sign;
      if level = 1 then begin
        res := res + arg;
        while not (stEmpty) and (cnt = level) do begin
          pop (s);
          Inc (level);
        end;
        if cnt <> level then begin
          Inc (cnt);
          sign := -sign;
          findCol;
          if level <> n then
            arg := prev^.arg;
        end else
        if cnt = n then
          break;
      end else begin
        push (s);
        cnt := 1;
        sign := 1;
        Dec (level);
        findStr;
      end;
    until false;
  end;
  definitor := res;
end;
var i, j: integer;
begin
  randomize;
  str0 := 1;
  for i := 1 to n do
    for j := 1 to n do
      matr [i, j] := round (random (10));
  { matr := matr0; }
  writeln ('Start !!!');
  writeln (definitor: 12: 0);
  for i := 1 to n do begin
    for j := 1 to n do
      write (matr [i, j]: 2: 0);
    writeln;
  end;
end.

 
 
 
 
Сообщение22.02.2006, 18:48 
Большое СПАСИБО!!!
Замечательная программа!

 
 
 
 
Сообщение22.02.2006, 18:50 
Аватара пользователя
:evil:
Пркиньте сложность рекурсивного определения -- не O(n!) ли? Тогда 100 x 100 Вы будете долгохонько считать. Проще, конечно, начать приводить матрицу к треугольной (a la метод Гаусса), и, вместо обратного прохода перемножить диагональ. Деструктивно, конечно, но одной копии хватит.

 
 
 
 
Сообщение22.02.2006, 19:35 
Аватара пользователя
yog писал(а):
...
Замечательная программа!

Ничего замечательного.

Вот это: http://alglib.sources.ru/matrixops/general/det.php для решения задачи подойдёт ?
Вот ещё: http://www.srcc.msu.su/num_anal/lib_na/cat/cat532.htm

 
 
 
 Re: Маленький вопрос по С++. Детерминант.
Сообщение22.02.2006, 21:56 
Аватара пользователя
yog писал(а):
Я ещё прикидывал. Для того, чтобы сосчитать дет. матрицы 10*10 надо создать где-то 2,6 млн. массивов и выделить около 77Мб оперативы :shock:. Может я иду не тем путём? :?:


Господь с Вами, страсти какие! Его "руками" можно за приемлемое время вычислить.

 
 
 
 
Сообщение23.02.2006, 10:28 
Аватара пользователя
Надо копать книжки - наверняка есть алгоритмы и пошустрее - у меня в MatLAB только что (на компе c процессором Athlon 1700+ ) определитель для матрицы 1000x1000 считал 860 миллисекунд. Неплохо да? Цикл, в котором я задавал матрицу значительно дольше крутился... Странно.

Правда функция det является built-in - внутрь никак не влезть :cry:

 
 
 
 
Сообщение23.02.2006, 16:10 
Аватара пользователя
:evil:
photon писал(а):
... у меня в MatLAB только что (на компе c процессором Athlon 1700+ ) определитель для матрицы 1000x1000 считал 860 миллисекунд. Неплохо да? Цикл, в котором я задавал матрицу значительно дольше крутился...

Ничуть не странно. Систему такого типа обычно интерпретируют Вашу программу, и уж, во всяком случае, не оптимизируют ее. В тоже время built-in функции обычно выписаны вручную на языке низкого уровня (С++), и оптимизированны дважды -- компилятором и человеком. Отсюда и разница в скорости.

 
 
 
 
Сообщение23.02.2006, 16:22 
Аватара пользователя
незванный гость писал(а):
Ничуть не странно. Систему такого типа обычно интерпретируют Вашу программу, и уж, во всяком случае, не оптимизируют ее. В тоже время built-in функции обычно выписаны вручную на языке низкого уровня (С++), и оптимизированны дважды -- компилятором и человеком. Отсюда и разница в скорости.

Различие в уровне оптимизации - это понятно, я просто не ожидал, что это будет настолько - тупое задание в цикле random, как
Код:
for j=1:1000, for j=1:1000, a(i,j)=rand(1);end,end

потребовало 4.597 с.

Без циклов, как
Код:
a=rand([1000 1000])
- потебовало, конечно, гораздо меньше времени (rand - тоже встроенная функция)- всего 60 мс.

 
 
 
 
Сообщение23.02.2006, 20:04 
Аватара пользователя
:evil:
Дело не в оптимизации -- в интерпретации. Я думаю, что MatLAB все-таки не компилирует, а интерпретирует программу.

 
 
 
 
Сообщение26.02.2006, 13:33 
Аватара пользователя
По-моему, MATLAB работает подобно Яве - сначала транслирует файлы в байт-код, а потом их выполняет.
Вот такое выражение:
Код:
for j=1:1000;
   for i=1:1000;
      a(i,j)=rand(1);
   end;
end;

у меня выполнялось аж 50 секунд с хвостиком (на Pentium-III, 933 MHz, MATLAB v. 5.3);
А вот "почти" то же самое:
Код:
a = zeros(1000, 1000);
for j=1:1000;
   for i=1:1000;
      a(i,j)=rand(1);
   end;
end;

выполнялось всего около 10 сек. (засекал на наручных часах по секундной стрелке. В этом преимущество древнего железа - для профилировки достаточно простейшего ПО и секундомера :) ) Связано это с тем, что при выполнении первого варианта внутри цикла каждый раз идет неявное перераспределение памяти (к массиву добавляются новые эл-ты), а во втором варианте память уже выделена.
Хотя, возможно что при прогоне этого кода под шестым MATLAB-ом такой заметной разницы не будет.

 
 
 
 
Сообщение27.02.2006, 10:07 
Аватара пользователя
Я давно для себя открыл такую разницу в коде - в 6-ом и 7-ом MatLABe - тоже самое - вы динамически наращиваете размер массива - это очень медленная операция, поэтому если приходится работать с большими массивами, вы заранее можете оценить их размер снизу, лучше сформировать его, используя zeros() или ones(), а потом забивать их нужными значениями.

Маленький совет - не засекайте по наручным часам. В MatLAB для этого есть специальные функции: tic - запускает сеундомер, toc - выводит время от момента его запуска.

 
 
 
 
Сообщение27.02.2006, 10:40 
Аватара пользователя
:) Спасибо за совет. Бум пользоватся в случае чего.

 
 
 
 
Сообщение27.02.2006, 11:11 
Аватара пользователя
Кстати в 6-ом Матлабе (не помню, было ли в 5.3) есть такая штука Profiler - запускаете через него свой код и потом можете посмотреть: сколько раз вызывалась какая функция и какие (по времени и в процентах от общего времени) были затраты на каждую строку. Иногда, если вычисления длительные очень помогает:1) Бывает, что увидишь, что простая операция присвоения съедает основную цасть времени - глядь, а ты как раз занят тем самым динамическим наращиванием массива, вставляешь одну строку, инициализирующую массив, и раз в 10 выигрываешь в скорости + выиграть как правило можно, уходя, если это возможно, от циклов к работе с массивами; 2) Иногда видно, что нет смысла биться над оптимизацией, потому что сколько-то твоих строк кода занимают какой-то 1% от общих временных затрат, несмотря на то, что они не оптимизированы, а основное время занимает какой-нибудь eig() или femeig, в которые влезть и улучшить нельзя

 
 
 [ Сообщений: 17 ]  На страницу 1, 2  След.


Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group