在多层循环中替代凌乱的网格节点索引

计算科学 正则 网格 朱莉娅
2021-11-30 05:49:59

最近,我深入研究了一组不知何故古老的 Fortran 代码,并试图完全理解它们。这些代码中有很大一部分是在许多状态变量维度上的多层循环,旨在找到策略函数(不动点)。没关系,但真正让我恼火的是网格节点索引的凌乱表达。例如,

do z=1,notimes
   do i=1,weignosims
      do t=1,nt
        ! main body of operations 
        ! nodes index below
           node = (i-1)*(nt+ntr) + t                               ! To keep track of an individual's life
           index = (begnosims -1)*notimes + weignosims*(z-1) + i   ! To keep track of which individual in which economy 
      enddo
   enddo
endo 

似乎作者希望将所有索引堆叠到一个数组中并跟踪它们中的每一个。但是,它自然是一个多重数组。这样的表达既令人困惑又容易出错。此外,在后面的情况下,作者还编写了类似asset(index,age)甚至定义三维数组来表示 1 个状态/节点。所以他为什么费心用上面的方式写

我想知道的是仅仅是编程风格还是这种风格可以获得一些好处(至少对于Fortran而言)?如果第二个原因成立,那么当我切换到一些更现代的科学计算语言(如 Julia)时,我可以有一些替代方法来“跟踪节点”吗?

1个回答

由于许多原因,这不是进行现代编程的好方法。首先,正如您所指出的,这种代码很难阅读和维护。其次,由于不再是问题的原因,这往往是在旧版本的 Fortran 中完成的。也就是说,旧版本的 Fortran 没有办法将对象“捆绑在一起”。如果您为所有内容制作了不同的数组,那么您必须将它们一起传递到函数中,从而导致大量可维护(并且容易出错)的函数签名。但最后也是最重要的一点,这会干扰编译器优化。如果您以简单的方式循环连续数组(循环展开、SIMD 等),则会发生许多编译器优化。通过使用这种奇怪的索引,有'

因此,如果您使用像 Julia 这样的现代语言,您可能不应该设计这种风格的大型程序。相反,您应该利用 Julia 的类型系统和多重分派。如果您自然有多个数组与一个项目相关联,请将它们组合成一个类型。编写assert要在此类型上调度的函数。等等。这将更容易阅读,并且可以让你以更简单的方式循环你的连续数组,这将允许更多的编译器优化和更好的性能(事实上,我已经在使用 DifferentialEquations.jl 的基准测试中击败了许多经典的 Fortran 代码,所以正确编写的 Julia 代码可以做到这一点!)

请注意,在许多使用 Fortran 循环的情况下,.可以改用 Julia 的广播语法(在任何情况下,您“循环遍历数组的所有索引”)。与 MATLAB 不同,Julia 的广播将融合循环并基本上执行我提到的优化以使循环与循环一样优化(并允许就地更新而没有临时数组,这与 MATLAB/NumPy 风格的矢量化代码中通常看到的不同)。这可以在不牺牲性能的情况下大大提高可读性。