Fortran Array to C array. Stupid macro tricks wanted

Question

I have this 'simplified' fortran code

real B(100, 200) 
real A(100,200)

... initialize B array code.

do I = 1, 100
  do J = 1, 200
    A(J,I) = B(J,I)
  end do
end do

One of the programming gurus warned me, that fortran accesses data efficiently in column order, while c accesses data efficiently in row order. He suggested that I take a good hard look at the code, and be prepared to switch loops around to maintain the speed of the old program.

Being the lazy programmer that I am, and recognizing the days of effort involved, and the mistakes I am likely to make, I started wondering if there might a #define technique that would let me convert this code safely, and easily.

Do you have any suggestions?

Answer 1

In C, multi-dimensional arrays work like this:

#define array_length(a) (sizeof(a)/sizeof((a)[0]))
float a[100][200];
a[x][y] == ((float *)a)[array_length(a[0])*x + y];

In other words, they're really flat arrays and [][] is just syntactic sugar.

Suppose you do this:

#define at(a, i, j) ((typeof(**(a)) *)a)[(i) + array_length((a)[0])*(j)]
float a[100][200];
float b[100][200];
for (i = 0; i < 100; i++)
    for (j = 0; j < 200; j++)
        at(a, j, i) = at(b, j, i);

You're walking sequentially through memory, and pretending that a and b are actually laid out in column-major order. It's kind of horrible in that a[x][y] != at(a, x, y) != a[y][x], but as long as you remember that it's tricked out like this, you'll be fine.

Edit

Man, I feel dumb. The intention of this definition is to make at(a, x, y) == at[y][x], and it does. So the much simpler and easier to understand

#define at(a, i, j) (a)[j][i]

would be better that what I suggested above.

Answer 2

Are you sure your FORTRAN guys did things right?

The code snippet you originally posted is already accessing the arrays in row-major order (which is 'inefficient' for FORTRAN, 'efficient' for C).

As illustrated by the snippet of code and as mentioned in your question, getting this 'correct' can be error prone. Worry about getting the FORTRAN code ported to C first without worrying about details like this. When the port is working - then you can worry about changing column-order accesses to row-order accesses (if it even really matters after the port is working).