fortran - Fortran 中高效的 z 顺序转换

Question

对于我目前在网格生成算法上的工作，我需要一种有效的方法来将三维坐标转换为 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 中更有效地进行位交错，那将非常有帮助。

Answer 1

可以使用类似于描述的查找表来优化从 3 个坐标到 z 顺序的转换 here .由于您只使用输入值的 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 个值用作数组中的索引，并使用移位和加法组合查找的值以交错这些值一起。在这种情况下，加法相当于按位或操作，因为在每个位位置最多设置一个位的情况下，不会有任何进位。因为在表中的值中只能设置每第 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 位)的查找表，用于存储重构输入值的 5 位。您将必须进行 12 次查找，一次为您的三个 20 位值中的每一个获取 5 位。屏蔽掉低 15 位，然后右移 0、1 和 2 位以获得 k、j 和 i 的查找索引。然后移位和掩码得到15-29、30-44和45-59位，每次都这样做，移位和相加重构k、j和i。