Fortran 中高效的 z 顺序转换

Efficient z-order transformation in Fortran

对于我目前在网格生成算法上的工作,我需要一种有效的方法来将三维坐标转换为 z 顺序(更准确地说:三个 4 字节整数转换为一个 8 字节整数),反之亦然.这篇维基百科文章描述得相当好: Z-order curve。 由于我不是程序员,我想出的解决方案做了它应该做的事情,但使用 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 优化后的代码已经快了很多,但它仍然是我程序的瓶颈。 如果有人能指出正确的方向如何在 Fortran 中更有效地进行位交错,那将非常有帮助。

从 3 坐标到 z 顺序的转换可以使用与描述的 here 类似的查找 table 进行优化。由于您只使用输入值的 20 位,因此使用具有 1024 个条目而不是 256 个条目的查找 table 会更有效,足以索引 10 位,因此您只需执行 2 次查找 - ups 为你的 3 个输入值中的每一个输入值,并针对交错 3 个值而不是 2 个值的情况进行了修改。

数组的条目 n 存储整数 n,其位展开以便位 0 在位 0 中,位1 被移动到位 3,位 2 被移动到位 6 等等,所有剩余的位都设置为零。查找 table 数组可以这样初始化:

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 个值用作数组中的索引,并使用移位组合查找的值并添加到交错价值观在一起。在这种情况下,加法相当于按位或操作,因为在每个位位置最多设置一个位的情况下,不会有任何进位。因为只能在 table 中的值中设置每第 3 位,所以将一个值偏移 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 位)的查找 table,用于存储重构输入值的 5 位。您将必须进行 12 次查找,一次为您的三个 20 位值中的每一个获取 5 位。屏蔽掉低 15 位,然后右移 0、1 和 2 位以获得 k、j 和 i 的查找索引。然后移位和掩码得到位15-29、30-44和45-59并且每次都做同样的事情,移位和添加以重建k,j和i。