Быстрая индексация симметричной матрицы

Question

В большом коде, написанном на Fortran08 для расчета термодинамических равновесий и фазовых диаграмм, я использую много симметричных матриц, которые я сохраняю как одномерные массивы, и индексирую с помощью небольшой функции

integer function ixsym(i,j)
   if(i.gt.j) then
      ixsym=j+i*(i-1)/2
   else
      ixsym=i+j*(j-1)/2
   endif
   return
end

Это работает отлично, но после улучшенияПо скорости выполнения других частей кода эта процедура теперь занимает 15-20% времени вычислений (она используется очень часто).Я предполагаю, что есть различные способы ускорить это, но я не знаю, C или другой способ заменить эту функцию, поэтому я ищу помощь.Я использую gfortran, но замена должна быть портативной.

Бо Сандман

Answer 1

Единственное, что вы можете рассмотреть, это избавиться от ветвления в этой функции:

Минимум и максимум двух чисел можно вычислить как:

max = (a+b + abs(a-b))/2
min = (a+b - abs(a-b))/2 = a+b - max

Так что вы можете использовать это как

integer function ixsym(i,j)
   integer :: p, q
   q = i+j; p = (q + abs(i-j))/2; q = q - p
   ixsym = q + (p*(p-1))/2        
   return
end

, который вы можете уменьшить как

integer function ixsym(i,j)
   integer :: p
   ixsym = i+j; p = (ixsym + abs(i-j))/2;
   ixsym = ixsym + (p*(p-3))/2        
   return
end

Answer 2

Компиляторы Фортрана раньше имели оптимизацию на уровне или лучше, чем компиляторы Си. Поэтому я не ожидал бы выигрыша, просто переключая язык, и скорее сосредоточился бы на улучшении алгоритмов.

• Как насчет замены расчета преобразования индекса на справочную таблицу? Есть ли у вас память для хранения значений ixsym для заданных индексов i и j?
  Да, он компенсирует вашу память для компромисса процессора, но если у вас много матриц, эта дополнительная может помочь.

• Действительно ли необходимо постоянно вычислять преобразование? Например. если вы перебираете элементы: ixsym(i, j+1) = ixsym(i, j) + 1, если i < j.

• Другая идея, хотя и зависит от оборудования, может заключаться в том, чтобы упорядочить данные по-другому, чтобы они оставались в пределах областей кэша ЦП. ( Link )

О вашем преобразовании индекса:

Сначала я думал, что вы использовали какой-то вариант функции связывания Кантора для перечисления вашего симметричного двумерного массива. Я попросил мою подругу Руби построить несколько пар, и она сказала мне:

(0, 0) ->  0  (0, 1) ->  0  (0, 2) ->  1  (0, 3) ->  3  (0, 4) ->  6
(1, 0) ->  0  (1, 1) ->  1  (1, 2) ->  2  (1, 3) ->  4  (1, 4) ->  7
(2, 0) ->  1  (2, 1) ->  2  (2, 2) ->  3  (2, 3) ->  5  (2, 4) ->  8
(3, 0) ->  3  (3, 1) ->  4  (3, 2) ->  5  (3, 3) ->  6  (3, 4) ->  9
(4, 0) ->  6  (4, 1) ->  7  (4, 2) ->  8  (4, 3) ->  9  (4, 4) -> 10

Я бы ожидал только два вхождения рассчитанного индекса, но я вижу три для некоторых пар. Это предназначено?

Обновление:

Это было начало индекса, как отметил один из пользователей Жан-Клод Арбо в своем комментарии. Вот ответ Руби с индексами, начинающимися с 1:

(1, 1) ->  1  (1, 2) ->  2  (1, 3) ->  4  (1, 4) ->  7  (1, 5) -> 11
(2, 1) ->  2  (2, 2) ->  3  (2, 3) ->  5  (2, 4) ->  8  (2, 5) -> 12
(3, 1) ->  4  (3, 2) ->  5  (3, 3) ->  6  (3, 4) ->  9  (3, 5) -> 13
(4, 1) ->  7  (4, 2) ->  8  (4, 3) ->  9  (4, 4) -> 10  (4, 5) -> 14
(5, 1) -> 11  (5, 2) -> 12  (5, 3) -> 13  (5, 4) -> 14  (5, 5) -> 15