Macros for 3D loops in C

Question

I'm developing a C (C99) program that loops heavily over 3-D arrays in many places. So naturally, the following access pattern is ubiquitous in the code:

for (int i=0; i<i_size, i++) {
    for (int j=0; j<j_size, j++) {
        for (int k=0; k<k_size, k++) {
            ...
        }
    }
}

Naturally, this fills many lines of code with clutter and requires extensive copypasting. So I was wondering whether it would make sense to use macros to make it more compact, like this:

#define BEGIN_LOOP_3D(i,j,k,i_size,j_size,k_size) \
for (int i=0; i<(i_size), i++) { \
    for (int j=0; j<(j_size), j++) { \
        for (int k=0; k<(k_size), k++) { 

and

#define END_LOOP_3D }}}

On one hand, from a DRY principle standpoint, this seems great: it makes the code a lot more compact, and allows you to indent the contents of the loop by just one block instead of three. On the other hand, the practice of introducing new language constructs seems hideously ugly and, even though I can't think of any obvious problems with it right now, seems alarmingly prone to creating bugs that are a nightmare to debug.

So what do you think: do the compactness and reduced repetition justify this despite the ugliness and the potential drawbacks?

Answer 1

Never put open or close {} inside macros. C programmers are not used to this so code gets difficult to read.

In your case this is even completely superfluous, you just don't need them. If you do such a thing do

FOR3D(I, J, K, ISIZE, JSIZE, KSIZE)      \
for (size_t I=0; I<ISIIZE, I++)          \
    for (size_t J=0; J<JSIZE, J++)       \
        for (size_t K=0; K<KSIZE, K++)

no need for a terminating macro. The programmer can place the {} directly.

Also, above I have used size_t as the correct type in C for loop indices. 3D matrices easily get large, int arithmetic overflows when you don't think of it.

Answer 2

The best way is to use a function. Let the compiler worry about performance and optimization, though if you are concerned you can always declare functions as inline.

Here's a simple example:

#include <stdio.h>
#include <stdint.h>

typedef void(*func_t)(int* item_ptr);


void traverse_3D (size_t  x, 
                  size_t  y, 
                  size_t  z, 
                  int     array[x][y][z],
                  func_t  function)
{
  for(size_t ix=0; ix<x; ix++)
  {
    for(size_t iy=0; iy<y; iy++)
    {
      for(size_t iz=0; iz<z; iz++)
      {
        function(&array[ix][iy][iz]);
      }
    }
  }
}


void fill_up (int* item_ptr)  // fill array with some random numbers
{
  static uint8_t counter = 0;
  *item_ptr = counter;
  counter++;
}

void print (int* item_ptr)
{
  printf("%d ", *item_ptr);
}


int main()
{
  int arr [2][3][4];

  traverse_3D(2, 3, 4, arr, fill_up);
  traverse_3D(2, 3, 4, arr, print);
}

EDIT

To shut up all speculations, here are some benchmarking results from Windows. Tests were done with a matrix of size [20][30][40]. The fill_up function was called either from traverse_3D or from a 3-level nested loop directly in main(). Benchmarking was done with QueryPerformanceCounter().

Case 1: gcc -std=c99 -pedantic-errors -Wall

With function, time in us:      255.371402
Without function, time in us:   254.465830

Case 2: gcc -std=c99 -pedantic-errors -Wall -O2

With function, time in us:      115.913261
Without function, time in us:   48.599049

Case 3: gcc -std=c99 -pedantic-errors -Wall -O2, traverse_3D function inlined

With function, time in us:      37.732181
Without function, time in us:   37.430324

Why the "without function" case performs somewhat better with the function inlined, I have no idea. I can comment out the call to it and still get the same benchmarking results for the "without function" case.

The conclusion however, is that with proper optimization, performance is most likely a non-issue.

Answer 3

If these 3D arrays are “small”, you can ignore me. If your 3D arrays are large, but you don't much care about performance, you can ignore me. If you subscribe to the (common but false) doctrine that compilers are quasi-magical tools that can poop out optimal code almost irrespective of the input, you can ignore me.

You are probably aware of the general caveats regarding macros, how they can frustrate debugging, etc., but if your 3D arrays are “large” (whatever that means), and your algorithms are performance-oriented, there may be drawbacks of your strategy that you may not have considered.

First: if you are doing linear algebra, you almost certainly want to use dedicated linear algebra libraries, such as BLAS, LAPACK, etc., rather than “rolling your own”. OpenBLAS (from GotoBLAS) will totally smoke any equivalent you write, probably by at least an order of magnitude. This is doubly true if your matrices are sparse and triply true if your matrices are sparse and structured (such as tridiagonal).

Second: if your 3D arrays represent Cartesian grids for some kind of simulation (like a finite-difference method), and/or are intended to be fed to any numerical library, you absolutely do not want to represent them as C 3D arrays. You will want, instead, to use a 1D C array and use library functions where possible and perform index computations yourself (see this answer for details) where necessary.

Third: if you really do have to write your own triple-nested loops, the nesting order of the loops is a serious performance consideration. It might well be that the data-access pattern for ijk order (rather than ikj or kji) yields poor cache behavior for your algorithm, as is the case for dense matrix-matrix multiplication, for example. Your compiler might be able to do some limited loop exchange (last time I checked, icc would produce reasonably fast code for naive xGEMM, but gcc wouldn't). As you implement more and more triple-nested loops, and your proposed solution becomes more and more attractive, it becomes less and less likely that a “one loop-order fits all” strategy will give reasonable performance in all cases.

Fourth: any “one loop-order fits all” strategy that iterates over the full range of every dimension will not be tiled, and may exhibit poor performance.

Fifth (and with reference to another answer with which I disagree): I believe, in general, that the “best” data type for any object is the set with the smallest size and the least algebraic structure, but if you decide to indulge your inner pedant and use size_t or another unsigned integer type for matrix indices, you will regret it. I wrote my first naive linear algebra library in C++ in 1994. I've written maybe a half dozen in C over the last 8 years and, every time, I've started off trying to use unsigned integers and, every time, I've regretted it. I've finally decided that size_t is for sizes of things and a matrix index is not the size of anything.

Sixth (and with reference to another answer with which I disagree): a cardinal rule of HPC for deeply nested loops is to avoid function calls and branches in the innermost loop. This is particularly important where the op-count in the innermost loop is small. If you're doing a handful of operations, as is the case more often than not, you don't want to add a function call overhead in there. If you're doing hundreds or thousands of operations in there, you probably don't care about a handful of instructions for a function call/return and, therefore, they're OK.

Finally, if none of the above are considerations that jibe with what you're trying to implement, then there's nothing wrong with what you're proposing, but I would carefully consider what Jens said about braces.