Optimizing a 3D-Loop (C++)

Question

I'm working on a multigrid solver in C++ and right now I'm trying to improve the serial performance. Most time-consuming part in this is the smoother, which in my case is a successive over-relaxation solver. This looks as follows (I hope it is quite self explanatory):

int idx; 
int strideY = stride_[level][0];
int strideZ = stride_[level][1];
for(int i = 0; i < steps; ++i) {

    for(int z = 1; z <= innerGridpoints_[level][2]; ++z) {
        for(int y = 1; y <= innerGridpoints_[level][1]; ++y) {
            idx = getIndexInner(level, 1,y,z);
            for(int x = 1; x <= innerGridpoints_[level][0]; ++x, ++idx) {
                grid[idx] = (1. - omega)  * grid[idx] + omega * 1./6. * (grid[idx+1] + grid[idx-1] +
                                    grid[idx + strideY]  + grid[idx - strideY] + 
                                    grid[idx + strideZ]  + grid[idx - strideZ] - 
                                    spacing_[level] * spacing_[level] * rhs[idx]);
            }
        }
    }
}

I already made a few optimizations: The loops are positioned such that the inner loop gives the most local entries (i.e. neighboring elements are along the x dimension), and the precalculation of idx (even though this is an inline function, it saved quite some time this way). I also tried blocking, i.e., not iterating over the whole grid, but only over small chunks to increase locality, but that did not have any influence. Last idea I have is to try some loop unrolling, but I actually don't expect big improvements from that. I was thinking that maybe there are some possible improvements towards memory access. Any tipps welcome :)

Just FYI: The gridsize will vary from very small up to 255x255x255. Also, the grid has some boundaries in every dimension, consisting of a small number of rows, i.e., the iteration is not over the whole grid.

Answer 1

A good optimizing compiler will do most of the easy things for you anyway, so always measure if the changes you're making actually improve things. And, inspect (and learn to understand) the generated assembly code to see what the compiler is actually doing.

But there are a few things I'd try, as the expression is complex, and even good optimizers sometimes need a bit of help:-

First, hoisting subexpressions that are invariant within the inner loop out to the surrounding loop. In you example, the obvious ones are spacing_[level] * spacing_[level] and omega * 1./6.

The other thing to try is to make idx be a pointer rather than an array index, and increment the pointer in your loop.

 int *idx = &grid[getIndexInner(level, 1,y,z)];  // assuming grid is array of ints.

Your expression then starts to look like this

*idx = (1. - omega)  * *idx + omega * 1./6. * (idx[1] + idx[-1] +
                                idx[strideY]  + idx[- strideY] + // etc...

Your optimizer (assuming it's turned on???) may well be doing this already. But it's worth a shot. As I said, without measurements this is a pointless exercise.

And, as @AkiSuihkonen mentions in the comments above "First make it work". Debugging highly optimized code is much more difficult, so make sure you have your algorithm performing exactly how it should be before you start worrying about performance.