Эффективное преобразование z-порядка в Фортране

Question

Для моей текущей работы над алгоритмом генерации сетки мне нужен эффективный способ преобразования трехмерных координат в z-порядок (точнее: три 4-байтовых целых числа в одно 8-байтовое целое число) и наоборот. Эта статья в Википедии описывает это довольно хорошо: кривая Z-порядка. Поскольку я не программист, то решение, которое я придумал, делает то, что должно, но может быть довольно наивным, используя встроенную функцию mvbits для явного чередования битов:

SUBROUTINE pos_to_z(i, j, k, zval)

use types

INTEGER(I4B), INTENT(IN)  :: i, j, k
INTEGER(I8B), INTENT(OUT) :: zval
INTEGER(I8B) :: i8, j8, k8
INTEGER(I4B) :: b

zval = 0
i8 = i-1
j8 = j-1
k8 = k-1

do b=0, 19
    call mvbits(i8,b,1,zval,3*b+2)
    call mvbits(j8,b,1,zval,3*b+1)
    call mvbits(k8,b,1,zval,3*b  )
end do

zval = zval+1

END SUBROUTINE pos_to_z


SUBROUTINE z_to_pos(zval, i, j, k)

use types

INTEGER(I8B), INTENT(IN)  :: zval
INTEGER(I4B), INTENT(OUT) :: i, j, k
INTEGER(I8B) :: i8, j8, k8, z_order
INTEGER(I4B) :: b

z_order = zval-1
i8 = 0
j8 = 0
k8 = 0

do b=0, 19
    call mvbits(z_order,3*b+2,1,i8,b)
    call mvbits(z_order,3*b+1,1,j8,b)
    call mvbits(z_order,3*b  ,1,k8,b)
end do

i = int(i8,kind=I4B) + 1
j = int(j8,kind=I4B) + 1
k = int(k8,kind=I4B) + 1

END SUBROUTINE z_to_pos

Обратите внимание, что я предпочитаю, чтобы диапазоны ввода и вывода начинались с 1 вместо 0, что приводит к некоторым дополнительным вычислениям. Оказывается, эта реализация довольно медленная. Я измерил время, необходимое для преобразования и повторного преобразования 10^7 позиций:
gfortran -O0: 6,2340 секунды
gfortran -O3: 5,1564 секунды
ifort -O0: 4.2058 секунд
ifort -O3: 0,9793 секунды

Я также безуспешно пробовал разные варианты оптимизации для gfortran. Хотя оптимизированный код с помощью ifort уже намного быстрее, он все еще является узким местом моей программы. Было бы очень полезно, если бы кто-то мог указать мне правильное направление, как сделать чередование чередования более эффективным в Фортране.

Answer 1

Преобразование из 3-х координат в z-порядок может быть оптимизировано с использованием справочной таблицы, аналогичной описанной здесь. Поскольку вы используете только 20 бит входных значений, было бы более эффективно использовать справочную таблицу с 1024 записями, а не с 256, что достаточно для индексации 10 битов, так что вам нужно всего лишь сделать 2 поиска для каждого из ваши 3 входных значения и изменены для случая чередования 3 значений вместо 2.

Запись n массива хранит целое число n с разложенными битами, так что бит 0 находится в бите 0, бит 1 перемещается в бит 3, бит 2 перемещается в бит 6 и так далее, со всеми оставшимися битами, установленными в нуль. Массив таблицы поиска может быть инициализирован так:

subroutine init_morton_table(morton_table)
    integer(kind=8), dimension (0:1023), intent (out) :: morton_table
    integer :: b, v, z
    do v=0, 1023
        z = 0
        do b=0, 9
            call mvbits(v,b,1,z,3*b)
        end do
        morton_table(v) = z
    end do
end subroutine init_morton_table

Чтобы фактически чередовать значения, разделите ваши 3 входных значения на их младшие 10 битов и их старшие 10 битов, затем используйте эти 6 значений в качестве индексов в массиве и объедините искомые значения, используя сдвиги и сложения, чтобы чередовать значения вместе. Добавления эквивалентны побитовым операциям ИЛИ в этом случае, потому что не будет никаких переносов, учитывая, что в каждой битовой позиции будет установлен максимум один бит. Поскольку в значениях таблиц может быть установлен только каждый третий бит, смещение одного из значений на 1 бит, а другого на 2 означает, что коллизий не будет.

subroutine pos_to_z(i, j, k, zval, morton_table)
    integer, intent(in) :: i, j, k
    integer(kind=8), dimension (0:1023), intent (in) :: morton_table
    integer(kind=8), intent (out) :: zval
    integer(kind=8) :: z, i8, j8, k8

    i8 = i-1
    j8 = j-1
    k8 = k-1

    z = morton_table(iand(k8, 1023))
    z = z + ishft(morton_table(iand(j8, 1023)),1)
    z = z + ishft(morton_table(iand(i8, 1023)),2)
    z = z + ishft(morton_table(iand(ishft(k8,-10), 1023)),30)
    z = z + ishft(morton_table(iand(ishft(j8,-10), 1023)),31)
    zval = z + ishft(morton_table(iand(ishft(i8,-10), 1023)),32) + 1

end subroutine pos_to_z

Вы можете использовать аналогичную технику, чтобы пойти другим путем, но я не думаю, что она будет столь же эффективной. Создайте таблицу поиска из 32768 значений (15 бит), в которых хранятся 5 бит восстановленного входного значения. Вам нужно будет выполнить 12 поисков, получая по 5 бит за раз для каждого из ваших трех 20-битных значений. Замаскируйте 15 младших битов, затем сдвиньте вправо на 0, 1 и 2 бита, чтобы получить ваши индексы поиска для k, j и i. Затем сдвиг и маску, чтобы получить биты 15-29, 30-44 и 45-59 и делать то же самое каждый раз, сдвигая и добавляя, чтобы восстановить k, j и i.