我有一个代码,我想循环遍历网格中的所有点,并为每个点检查给定条件是否适用于足够数量的相邻点。此外,我在网格上有周期性边界。
这个问题和生命游戏非常相似。
我当前的代码看起来像这样
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
这既不优雅,也可能不是很有效。有一个更好的方法吗?这段代码会重复很多,所以我想尽快完成。
最佳答案
我会在 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
此处,当 my_array(i,j)
的邻域满足您的条件时,logical_function_of()
返回 true。在列出 N、S、E、W 邻居之后我感到疲倦,对于生产代码,我可能无论如何都会将其编写为索引的函数。以我的经验,forall
很优雅(对某些人来说)但不如嵌套循环那么高效。
关于optimization - 快速计算 3D 数组中相邻点的方法,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/11096887/