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

可以使用类似于here 描述的查找表来优化从 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 个值用作数组的索引，并使用移位和相加组合查找的值以交错值一起。在这种情况下，加法等同于按位或运算，因为在每个位位置最多设置一个位的情况下，不会有任何进位。因为表中的值只能设置每 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。