位操作-有效的z顺序转换在Fortran

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 - 0: 6.2340秒
gfortran - 3: 5.1564秒
[au:]

我也尝试了不同的优化选项的gfortran没有成功。虽然经过优化的代码已经快了很多，但它仍然是我的程序的瓶颈。如果有人能告诉我如何在Fortran中更有效地进行位交错，那将是非常有帮助的。

Answer 1

可以使用类似于这里描述的查找表来优化从3个协词到z顺序的转换。由于您只使用输入值的20位，因此使用具有1024项而不是256项的查找表将更有效，足以索引10位，因此您只需要为3个输入值中的每个值进行2次查找，并针对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。