arrays - 向量化访问非连续内存位置的循环

Question

我有一个这种结构的循环

引用:Maxwell Code Example

do z=1,zend
    do y=1,yend
        do x=1,xend
            k=arr(x,y,z)
            do while(k.ne.0)
                ix=fooX(k)
                iy=fooY(k)
                iz=fooZ(k)
                x1=x(ix  ,iy  ,iz)
                x2=x(ix+1,iy  ,iz)
                x3=x(ix  ,iy+1,iz)
                x4=x(ix+1,iy+1,iz)
                x5=x(ix  ,iy  ,iz+1)
                x6=x(ix+1,iy  ,iz+1)
                x7=x(ix  ,iy+1,iz+1)
                x8=x(ix+1,iy+1,iz+1)

                y1=y(ix  ,iy  ,iz)
                y2=y(ix+1,iy  ,iz)
                y3=y(ix  ,iy+1,iz)
                y4=y(ix+1,iy+1,iz)
                y5=y(ix  ,iy  ,iz+1)
                y6=y(ix+1,iy  ,iz+1)
                y7=y(ix  ,iy+1,iz+1)
                y8=y(ix+1,iy+1,iz+1)

                z1=z(ix  ,iy  ,iz)
                z2=z(ix+1,iy  ,iz)
                z3=z(ix  ,iy+1,iz)
                z4=z(ix+1,iy+1,iz)
                z5=z(ix  ,iy  ,iz+1)
                z6=z(ix+1,iy  ,iz+1)
                z7=z(ix  ,iy+1,iz+1)
                z8=z(ix+1,iy+1,iz+1)
                sumX+=x1+x2+..x8
                sumY+=y1+y2+..y8
                sumZ+=z1+z2+..z8

                k=linkArr(k)
            enddo
        enddo
    enddo
enddo

x1 到 x8 是长方体的 8 个角。向量化此代码面临三个挑战。一是8个数组元素在内存中不连续。其次是固有的 while 循环结构以及链表访问。第三，从 fooX、fooY、fooZ 返回的 ix、iy、iz 的值不是不连续的。所以循环的每次迭代都有一组完全不同的 ix, iy, iz。因此，即使在迭代中，内存访问也是分散的。我尝试了以下方法:1. 将 3 级 DO 循环展开为:​​

do z=1,zend
    do y=1,yend
        do x=1,xend  
           if(arr(x,y,z).NE.0) then
                kArr(indx)=arr(x,y,z)
                DO WHILE (kArr(indx).NE.0)
                  indx = indx + 1
                  kArr(indx)=linkArr(kArr(indx-1))
                ENDDO
            endif
        enddo
    enddo
enddo

有了这个，我摆脱了 while 循环结构，现在我可以在 kArr 上运行一个大循环，在其中我将 8 个元素分组(假设我的 VPU 可以一次容纳 8 组数据)。它没有提供性能改进。如果有人感兴趣，我可以发布这些细节。我需要有关如何优化此代码的建议。我尝试的另一个选择是将 x、y、z 数据合并到一个数组中，这样当我计算 x1 时，y1 和 z1 也将位于相邻的内存位置。

Answer 1

那个 while 循环会要了你的命。在几年前的类似情况下，我通过执行以下操作获得了适度的性能提升:

! at top of your code, introduce:
integer :: special_index
integer :: ix(1000), iy(1000), iz(1000)  !promoting scalars to arrays.
                                         ! make as big as possibly needed.

! code as usual until you get to your loops, then

! first, make lookup table
special_index=0
do z=1,zend
  do y=1,yend
    do x=1,xend
      k=arr(x,y,z)
      do while(k.ne.0)
        special_index=special_index+1
        ix(special_index)=fooX(k)
        iy(special_index)=fooY(k)
        iz(special_index)=fooZ(k)
        k=linkArr(k)
      enddo
    enddo
  enddo
endoo
! and now we do the calculation, loop over lookup table:
do n=1,special_index
  x1=x(ix(n)  ,iy(n)  ,iz(n))
  x2=x(ix(n)+1,iy(n)  ,iz(n))
  x3=x(ix(n)  ,iy(n)+1,iz(n))
  etc.
enddo

就像我说的，这在几年前帮助了我。你的旅费可能会改变。第一个循环仍然不会矢量化，但第二个循环可能会提供更好的性能。