Should the derivative of an array be calculated array by array or element by element in CFD codes?

Question

I am making my own finite difference computational magnetohydrodynamic code in Fortran 90. Looking at other codes they appear to calculate for example their $x$-derivatives, bb of their variables, e.g. aa like the following:

CALL system_clock(tstart, count_rate)
DO ix = 1, nx - 1
  ixp = ix + 1
  DO iy = 1, ny
    DO iz = 1, nz
      bb(ix,iy,iz) = (aa(ixp,iy,iz) - aa(ix,iy,iz)) / dx
    END DO
  END DO
END DO
CALL system_clock(tstop, count_rate)
PRINT*, "Computation time =", REAL(tstop - tstart) / REAL(count_rate)

In other words they calculate the derivative element by element. Another way of calculating the derivative is array by array, i.e.:

CALL system_clock(tstart, count_rate)
bb(1:nx-1,:,:) = (aa(2:nx,:,:) - aa(1:nx-1,:,:)) / dx
CALL system_clock(tstop, count_rate)
PRINT*, "Computation time =", REAL(tstop - tstart) / REAL(count_rate)

The second method seems to be about 37 times quicker. So why do I not see more codes calculating their $x$-derivatives array by array instead of element by element? What is the downside of array by array, is there something I am missing?

Answer 1

From an algorithmic point of view, the two versions are practically identical. What makes the difference is that Fortran90 is very good at large array ops and there are very mature compilers for this type of operation. In your second example, the compiler has more freedom to optimize for efficiency. I'm no expert on the inner workings of compilers but think of it this way: When you write:

b(1:nx-1,:,:) = (aa(2:nx,:,:) - aa(1:nx-1,:,:)) / dx

Then you tell your compiler that you want to do this operation for all the indices targeted with ":". You also kind of tell it that you do not care about the sequence in which the operations are to be done. The compiler can then optimize for your particular architecture, rearrange the operations, do smart prefetching and so on.

In your first example you are a bit more specific in the way you present the operation. You tell it that it should iterate in row/column-major, that the sequence in which to do this has to be the way you stated.

short: In Fortran, most of the times, you profit from giving the compiler reasonable freedom!

Why do so few people use it in this way? Well, it strongly depends on what you are optimizing for. If people have two days to implement a prototype, then these details may not matter. Some people may find the first version more expressive and human readable etc. Some people may only have a rudimentary understanding of how a compiler may optimize your instructions. You made exactly the right call by actually testing and comparing, because that what it often comes down to.

(Beware of index overflow! There might not always be a left/right neighbor for differentiation.)