Alternative to messy grid node indexing within multiple layers of loops

Question

Recently,I dive into a set of somehow ancient Fortran codes and try to fully understand them. A large fraction of these codes are multiple layers of loops over many state variable dimensions, which are aiming to find policy function(fixed point). That's OK, but what's really annoyed me is the messy expression of the grid node index. For example,

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 

It seems the author want to stack all index into one array and keep track each of them. However, it's naturally a multiple array. Such expression is both confusing and error-prone. Moreover, in later case, the author also write something like asset(index,age) or even define three dimension array to represent 1 state/node. So why he bother to write in the above way

What I want to know is whether is just a matter of programming style or this kind of style can gain some benefit (at least w.r.t Fortran)? If the second reason holds, then can I have some alternative way to "trace node" when I switch to some more modern scientific computation language, like Julia?

Answer 1

This is not a good way to do modern programming for many reasons. First of all, as you pointed out, this kind of code is hard to read and maintain. Secondly, this tends to be done in old versions of Fortran for reasons which are no longer a problem. Namely, old versions of Fortran didn't have ways to "bundle together" objects. If you made different arrays for everything, then you'd have to pass them all together into functions, leading to large maintainable (and prone to error) function signatures. But lastly and most importantly, this can interfere with compiler optimizations. There are a lot of compiler optimizations that occur if you loop over contiguous arrays (loop unrollling, SIMD, etc) in an easy fashion. By using such odd indexing, there's a good chance that this will cause the compiler to not add these (modern) optimizations, leading to worse performance.

So you probably shouldn't design large programs in this style if you are using a modern langaguage like Julia. Instead, you should exploit Julia's type system and multiple dispatch. If you naturally have multiple arrays associated with one item, put them together into a type. Write the assert function to dispatch on this type. Etc. This will be much easier to read and will could allow you to loop in an easier fashion over your contiguous arrays, which will allow more compiler optimizations and better performance (indeed, I have beat many classic Fortran codes in benchmarks with DifferentialEquations.jl, so properly written Julia code will do this!)

Note that in many cases where a Fortran loop is used, Julia's broadcasting . syntax can be used instead (in any case where you "loop over all indices of an array"). Unlike MATLAB, Julia's broadcast will fuse loops and essentially perform the optimizations I mentioned to make the as optimal as the loops (and allow for in-place updates with no temporary arrays unlike is normally seen in MATLAB/NumPy-style vectorized code). This could greatly enhance readability without sacrificing performance.