2014 dxdy logo

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

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




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


09/02/06
50
Киев
Господа!

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

 Профиль  
                  
 
 Re: Маленький вопрос по С++. Детерминант.
Сообщение22.02.2006, 17:20 
Аватара пользователя


20/01/06
64
оттуда
yog писал(а):
...
Может кто-то знает алгоритм по-проще?
...

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

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

 Профиль  
                  
 
 
Сообщение22.02.2006, 17:57 
Аватара пользователя


20/01/06
64
оттуда
Прошу прощения, 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 


09/02/06
50
Киев
Большое СПАСИБО!!!
Замечательная программа!

 Профиль  
                  
 
 
Сообщение22.02.2006, 18:50 
Заслуженный участник
Аватара пользователя


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

 Профиль  
                  
 
 
Сообщение22.02.2006, 19:35 
Аватара пользователя


20/01/06
64
оттуда
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 
Заслуженный участник
Аватара пользователя


23/07/05
17973
Москва
yog писал(а):
Я ещё прикидывал. Для того, чтобы сосчитать дет. матрицы 10*10 надо создать где-то 2,6 млн. массивов и выделить около 77Мб оперативы :shock:. Может я иду не тем путём? :?:


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

 Профиль  
                  
 
 
Сообщение23.02.2006, 10:28 
Экс-модератор
Аватара пользователя


23/12/05
12046
Надо копать книжки - наверняка есть алгоритмы и пошустрее - у меня в MatLAB только что (на компе c процессором Athlon 1700+ ) определитель для матрицы 1000x1000 считал 860 миллисекунд. Неплохо да? Цикл, в котором я задавал матрицу значительно дольше крутился... Странно.

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

 Профиль  
                  
 
 
Сообщение23.02.2006, 16:10 
Заслуженный участник
Аватара пользователя


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

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

 Профиль  
                  
 
 
Сообщение23.02.2006, 16:22 
Экс-модератор
Аватара пользователя


23/12/05
12046
незванный гость писал(а):
Ничуть не странно. Систему такого типа обычно интерпретируют Вашу программу, и уж, во всяком случае, не оптимизируют ее. В тоже время 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 
Заслуженный участник
Аватара пользователя


17/10/05
3709
:evil:
Дело не в оптимизации -- в интерпретации. Я думаю, что MatLAB все-таки не компилирует, а интерпретирует программу.

 Профиль  
                  
 
 
Сообщение26.02.2006, 13:33 
Заслуженный участник
Аватара пользователя


12/10/05
478
Казань
По-моему, 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 
Экс-модератор
Аватара пользователя


23/12/05
12046
Я давно для себя открыл такую разницу в коде - в 6-ом и 7-ом MatLABe - тоже самое - вы динамически наращиваете размер массива - это очень медленная операция, поэтому если приходится работать с большими массивами, вы заранее можете оценить их размер снизу, лучше сформировать его, используя zeros() или ones(), а потом забивать их нужными значениями.

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

 Профиль  
                  
 
 
Сообщение27.02.2006, 10:40 
Заслуженный участник
Аватара пользователя


12/10/05
478
Казань
:) Спасибо за совет. Бум пользоватся в случае чего.

 Профиль  
                  
 
 
Сообщение27.02.2006, 11:11 
Экс-модератор
Аватара пользователя


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

 Профиль  
                  
Показать сообщения за:  Поле сортировки  
Начать новую тему Ответить на тему  [ Сообщений: 17 ]  На страницу 1, 2  След.

Модераторы: Karan, Toucan, PAV, maxal, Супермодераторы



Кто сейчас на конференции

Сейчас этот форум просматривают: нет зарегистрированных пользователей


Вы не можете начинать темы
Вы не можете отвечать на сообщения
Вы не можете редактировать свои сообщения
Вы не можете удалять свои сообщения
Вы не можете добавлять вложения

Найти:
Powered by phpBB © 2000, 2002, 2005, 2007 phpBB Group