View source for C vs Fortran memory order

Everyone knows that in C, multidimensional arrays are indexed with the most rapidly changing index last: [k][j][i]. In Fortran arrays, the most rapidly changing index comes first: (i,j,k). Let's compare C memory ordering and Fortran memory ordering to see if there is a difference in how the data are actually stored in memory or whether the order is the same and it's just the language syntax that is different.

=Printing array values=
First, we're going to need a function to print a linear array with some number of values. Let's write that in C since it is convenient. We can call it later from Fortran too so let's make all of the arguments pass by reference.

<source lang="c">
#include <stdio.h>

void
print_array_values(const int *values, const int *nvalues)
{
    int i, c = 0;
    for(i = 0; i < *nvalues ; i++)
    {
        printf("%d, ", values[i]);
        if(++c >= 10)
        {
            c = 0;
            printf("\n");
        }
    }
    printf("\n");
}

/* Called from Fortran (the trailing underscore makes it easier to link) */
void
print_array_values_(const int *values, const int *nvalues)
{
    print_array_values(values, nvalues);
}
</source>

Our print function will take a linear array (remember that even multidimensional arrays are stored linearly in memory) and print its contents. Now we need to compile our function:

<source lang="bash">
gcc -c -o print_array_values.o print_array_values.c
</source>

=C memory order=
Here is our C example program. We create a 3D array values[NZ][NY][NX]. It's declared that way since that the most rapidly changing index should be last for optimum cache performance. The example program iterates through our indexing space in k,j,i order, storing an increasing index into each array element. When printed, we should see values starting with 0 and going up through 59 in increments of 1.

<source lang="c">
#define NX 3
#define NY 4
#define NZ 5

void print_array_values(const int *values, const int *nvalues);

int
main(int argc, char *argv[])
{
    int index = 0;
    int i,j,k,nvals;

    int values[NZ][NY][NX];
    for(k = 0; k < NZ; ++k)
        for(j = 0; j < NY; ++j)
            for(i = 0; i < NX; ++i)
                values[k][j][i] = index++;

    nvals = NX*NY*NZ;
    print_array_values((const int *)values, &nvals);

    return 0;
}
</source>

Compile the example:
<source lang="bash">
gcc -o indexing indexing.c print_array_values.o
</source>

Output:
<code><pre>
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 
10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 
30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 
</pre></code>

=Fortran memory order=
Now, let's see what a Fortran program will do. Remember that Fortran arrays are indexed using the most rapidly changing variable first. We want to see how that affects the memory layout of the array.

<source lang="fortran">
      program main
      implicit none
      integer NX,NY,NZ
      parameter (NX = 3)
      parameter (NY = 4)
      parameter (NZ = 5)
      integer values(NX,NY,NZ),i,j,k,nvals,index

      index = 0
      do 30 k=1,NZ
          do 20 j=1,NY
              do 10 i=1,NX
                  values(i,j,k) = index
                  index = index + 1
10            continue
20        continue
30    continue

      nvals = NX*NY*NZ
      call print_array_values(values, nvals)

      stop
      end
</source>

Compiling the example:
<source lang="bash">
gfortran -o indexing_f indexing.f print_array_values.o
</source>

Output:
<code><pre>
0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 
10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 
20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 
30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 
40, 41, 42, 43, 44, 45, 46, 47, 48, 49, 
50, 51, 52, 53, 54, 55, 56, 57, 58, 59, 
</pre></code>

=Summary=
Note that the output is exactly the same as in the C program!

This means that when the loops are iterated in the same order as in C (k,j,i) we get the same memory layout as in a C program, even though the language syntax requires (i,j,k) ordering. So, it's really just the indexing syntax that is different if the loops are treated as in C. In fact, if we were storing anything but an increasing index into the array, we'd get the same values in C and Fortran regardless of the order of the loops, but cache performance would vary.

'''Q:''' But why would I want to treat loops as in C? Why wouldn't I want to iterate in i,j,k order?<br>

'''A:''' Well, you could iterate in i,j,k order but you would no longer access memory in a contiguous manner. Our program demonstrates that indexing in k,j,i order, as in C, lets you access adjacent memory locations with each store so you'll have a better memory access pattern. If you were iterating in i,j,k order, the i index would no longer change most rapidly and your memory accesses would skip around inside the array, leading to poorer cache performance.

In summary:
*The "shape" of the memory is exactly the same in C and Fortran even though language differences make it look different due to reversed array indexing.
*If you don't iterate through a Fortran array in k,j,i order, you'll access memory out of order and negatively impact your cache performance.

[[Category: Developer documentation]]