2

我有一个代码，我想遍历网格中的所有点，并为每个点检查给定条件是否适用于足够数量的相邻点。此外，我在网格上有周期性边界。

这个问题与生命游戏非常相似。

我当前的代码看起来像这样

do k=1,ksize; do j=1,jsize; do i=1,isize ! Loop over all points
  ncount = 0
  kkloop: do kk=k-1,k+1 ! Loop over neighbours
    ktmp = kk
    if(kk>ksize) ktmp = 1 ! Handle periodic boundary
    if(kk<1) ktmp = ksize
    do jj=j-1,j+1
      jtmp = jj
      if(jj>jsize) jtmp = 1
      if(jj<1) jtmp = jsize
      do ii=i-1,i+1
        if(ii == 0 .and. jj == 0 .and. kk == 0) cycle ! Skip self
        itmp = ii
        if(ii>isize) itmp = 1
        if(ii<1) itmp = isize

        if(grid(itmp,jtmp,ktmp)) ncount = ncount + 1 ! Check condition for neighbour
        if(ncount > threshold) then ! Enough neigbours with condition?
          do_stuff(i,j,k)
          break kkloop
        end if
      end do
    end do
  end do
end do; end do; end do

这既不优雅，也不可能非常有效。有一个更好的方法吗？这段代码会重复很多，所以我想让它尽可能快。

4

3 回答 3

5

我会在 2D 中解决这个问题，让你自己膨胀到 3D。

我要做的第一件事是用深度光晕填充数组，该光晕的深度等于您感兴趣的邻域的深度。所以，如果您的数组被声明为，比如说

real, dimension(100,100) :: my_array

你对每个单元格的 8 个直接邻居感兴趣，

real, dimension(0:101,0:101) :: halo_array
.
.
.
halo_array(1:100,1:100) = my_array
halo_array(0,:) = my_array(100,:)
! repeat for each border, mutatis mutandis

这将节省大量检查边界的时间，无论您是否遵循下一个建议，都值得这样做。如果您愿意，您可以“就地”执行此操作，我的意思是扩展my_array而不是复制它。

对于一个优雅的解决方案，你可以写这样的东西

forall (i=1:100,j=1:100)
   if (logical_function_of(my_array(i-1,j),my_array(i+1,j),my_array(i,j-1),my_array(i,j+1),...) then
      do_stuff(my_array(i,j))
   end if
end forall

在这里，logical_function_of()当邻域my_array(i,j)满足您的条件时返回 true。列出 N、S、E、W 邻居后我感到很累，对于生产代码，无论如何我可能会将其编写为索引的函数。以我的经验forall（对某些人来说）是优雅的，但不如嵌套循环那么高性能。

于 2012-06-19T10:51:22.800 回答
1

您可以使用 kd-tree 或 octree 来细分 3d 数组。像 az morton 阶数曲线这样的空间填充曲线在 3d 中很有用，可以为八个立方体创建一个键。但这最适用于 2 个 3d 阵列的功率。

于 2012-06-19T10:04:46.237 回答
1

您为grid网格上的每个点调用该函数 27 次。如果这是一个昂贵的电话，你想打电话给它更少的次数。

例如，如果grid返回真值的概率很低，您可以为网格上的每个点调用它，并将条件成立的点存储在 kd-tree 中。然后更容易迭代 kd 树内的点，计算每个点的邻居。

否则，您可以使用维度的位矩阵isize*jsize*3来缓存所有点的网格值kk=k-1,k+1。如果这太大了，您可以使用一个大小矩阵isize*3*3进行缓存的中间解决方案。

于 2012-06-19T11:16:04.573 回答