C++ Optimize vectorizing nested loops

Question

I have a program dealing with multiple nested loops, operating on a 3D domain:

unsigned int sX(m_sizeZ*m_sizeY);
unsigned int b(sX+m_sizeZ);
for(unsigned int i(1);i<m_sizeX-1;i++){
    for(unsigned int j(1);j<m_sizeY-1;j++){     
        for(unsigned int k(1);k<m_sizeZ-1;k++){
            m_r[b+k]=m_func[b+k]-m_cX*(m_data[b+k-sX]+m_data[b+k+sX]-2.0*m_data[b+k])
                        -m_cY*(m_data[b+k-m_sizeZ]+m_data[b+k+m_sizeZ]-2.0*m_data[b+k])
                        -m_cZ*(m_data[b+k-1]+m_data[b+k+1]-2.0*m_data[b+k]);
        }
        b+=m_sizeZ;
    }
    b+=2*m_sizeZ;
}

Where my arrays are double of size m_sizeX*m_sizeY*m_sizeZ.

I iterate this way because I don't want to touch the boundaries of the domain.

When compiling with (g++) -msse2 -ftree-vectorizer-verbose=2, I of course get the multiple nested loops note.

Is there any way I could use instead a single loop without (more or less) complicated check operations?

Thanks!

Well you could probably code all 3 coords into one number and then decode it with modulos. Just like you can map 2D array into a single block of memory and vice-versa. However I would discourage such approach. The loop is fairly simple so the compiler should not have any troubles transforming it on its own. — luk32, May 29 '14 at 08:57
"the multiple nested loops note" Are you always this vague? Is it fun? — Lightness Races in Orbit, May 29 '14 at 09:28
Sorry if I looked vague. The complete message I'm getting is : PoissonSolverCG.cpp:93: note: not vectorized: multiple nested loops. PoissonSolverCG.cpp:93: note: bad loop form. Analyzing loop at PoissonSolverCG.cpp:94 — repptilia, May 29 '14 at 09:29

score 2 · Accepted Answer · answered May 29 '14 at 11:24

If your goal is good vectorization, it is probably best to just apply the same calculation to your edge points as to the inner points, only to reset them after you are done calculating all the points. Something like this:

unsigned int sX(m_sizeZ*m_sizeY);
unsigned int start = (1*m_sizeY + 1)*m_sizeZ + 1;
unsigned int end = ((m_sizeX - 1)*m_sizeY - 1)*m_sizeZ - 1;
//Do calculation for everything, including the edges.
for(unsigned int i = start; i < end; i++) {
    m_r[i]=m_func[i]-m_cX*(m_data[i-sX]+m_data[i+sX]-2.0*m_data[i])
                -m_cY*(m_data[i-m_sizeZ]+m_data[i+m_sizeZ]-2.0*m_data[i])
                -m_cZ*(m_data[i-1]+m_data[i+1]-2.0*m_data[i]);
}
//Reset the edges.
for(unsigned x = 0; x < m_sizeX; x++) {
    for(unsigned y = 0; y < m_sizeY; y++) {
        m_r[x*sX + y*m_sizeZ] = m_data[x*sX + y*m_sizeZ];
        m_r[x*sX + y*m_sizeZ + m_sizeZ-1] = m_data[x*sX + y*m_sizeZ + m_sizeZ-1];
    }
}
for(unsigned x = 0; x < m_sizeX; x++) {
    for(unsigned z = 0; z < m_sizeZ; z++) {
        m_r[x*sX + z] = m_data[x*sX + z];
        m_r[x*sX + (m_sizeY-1)*m_sizeZ + z] = m_data[x*sX + (m_sizeY-1)*m_sizeZ + z];
    }
}

This is additional calculations that would be done, but it has two positive effects:

It is now very easy for your compiler to vectorize the first loop (which takes most of the time).
This approach virtually eliminates the edge problem that results from the fixed vector size: Since your vector unit handles several alligned(!) loop iterations in one, every edge in your calculations leads to two special iterations that need to be done. One at the start of a run to get the vector loop aligned, and one at the end to handle the leftovers of the vector loop.

Thank you! this looks like a nice solution indeed. – repptilia May 29 '14 at 13:48 — repptilia, May 29 '14 at 13:48

Pablo Francisco Pérez Hidalgo · Answer 2 · 2016-08-15T12:12:30.933

1

You can iterate, in a single loop, from 1 to m_sizeX*m_sizeY*m_sizeZ (with counter C) and calculate i,j,k as:

i = C / (m_sizeY*m_sizeZ)
j = (C % (m_sizeY*m_sizeZ)) / m_sizeZ
k = (C % (m_sizeY*m_sizeZ)) % m_sizeZ

This method has the constraint that you have to take care of m_sizeX*m_sizeY*m_sizeZ being within C range without overflowing.

EDIT

To control your boundaries without using if-else clauses you can create a function

size_t nextToCalculate(size_t previous)
{
    return previous+1+!condition;
}

And use it in your loop:

for(int C = 0; C < m_sizeX*m_sizeY*m_sizeZ; C = nextToCalculate(C))
{
  int z = (C % (m_sizeY*m_sizeZ)) % m_sizeZ;
  int y = (C % (m_sizeY*m_sizeZ)) / m_sizeZ;
  int x = C / (m_sizeY*m_sizeZ);
  ...
  ...
  ...
}

Or even include its implementation in line:

for(int C = 0; C < m_sizeX*m_sizeY*m_sizeZ; C = C+1+!CONDITION(C+1))
{
  int z = (C % (m_sizeY*m_sizeZ)) % m_sizeZ;
  int y = (C % (m_sizeY*m_sizeZ)) / m_sizeZ;
  int x = C / (m_sizeY*m_sizeZ);
  ...
  ...
  ...
}

edited Aug 15 '16 at 12:12

answered May 29 '14 at 08:56

Pablo Francisco Pérez Hidalgo

27,044
8
36
62

I'd add a note to make sure that `C` can hold `m_sizeX*m_sizeY*m_sizeZ`. Especially after below comment. – luk32 May 29 '14 at 08:58
@luk32 I took it as obvious but you are right, I'll add that. – Pablo Francisco Pérez Hidalgo May 29 '14 at 08:59
@repptilia As long as `m_sizeX*m_sizeY*m_sizeZ` fits within `C` range you should not have any problem, neither performance issues. If that value is too big, you can use the same strategy to flat only two levels, mapping either `m_sizeX*m_sizeY` or `m_sizeY*m_sizeZ` – Pablo Francisco Pérez Hidalgo May 29 '14 at 09:04
I'm doing that. How do I take care of the boundaries then? If statements cause a note: not vectorized: control flow in loop. – repptilia May 29 '14 at 09:13
@repptilia Did you try the ternary operator `condition?ifyes:ifno`? On the other hand, you can change how k,j and i are calculated to filter some values. – Pablo Francisco Pérez Hidalgo May 29 '14 at 09:18
@repptilia I got an idea: if your loop is like `for(int C = 0; C < sx*sy*sz; C++)`, remove the `C++` changing it for: `for(int C = 0; C < sx*sy*sz; C = nextToCalculate(C))` where `nextToCalculate` skips thos values out of your boundaries. – Pablo Francisco Pérez Hidalgo May 29 '14 at 09:29
I'm getting note: not vectorized: latch block not empty. – repptilia May 29 '14 at 09:35
@repptilia The only remaining idea I have is: make `nextToCalculate` inline by `inline` reserved word or including its implementation within the for loop increment group. – Pablo Francisco Pérez Hidalgo May 29 '14 at 09:40
okay, still not working. Looks like i'll have to keep this un-vectorized. – repptilia May 29 '14 at 09:43
The `!` operator is just as bad to performance as an `if()` clause (it does the same thing under the hood). Consequently, this "trick" only serves to obfuscate the code. – cmaster - reinstate monica May 29 '14 at 09:44
@cmaster You don't have to use `!` operator, just write the contrary condition. I used `!` to denote that. – Pablo Francisco Pérez Hidalgo May 29 '14 at 09:45
1

@PabloFranciscoPérezHidalgo It doesn't matter which logical operator you use, all logical operators will compile to a comparison plus a conditional branch. This applies to `==`, `<=`, `!`, `&&`, and so on. CPUs don't like any of these. – cmaster - reinstate monica May 29 '14 at 09:50

score 0 · Answer 3 · answered May 29 '14 at 08:56

Put code as entire function which can be compiled. At first look:

use const int m_sizeX_1 = m_sizeX-1 and use it in for loop
use temporary variable for b+k
even use temporary variable for m_data+b+k
generaly advice for vectorization - simplify code with extracing variables and then you will se easier what to do next

score 0 · Answer 4 · answered May 29 '14 at 09:15

0

you may try this (single loop, as you requested):

unsigned int sX(m_sizeZ*m_sizeY);
unsigned int b(sX+m_sizeZ);
unsigned int i, j, k;
for (i = 1, j = 1, k = 1;
     i < m_sizeX-1 && j < m_sizeY - 1 && k < m_sizeZ - 1;
     k++) {
    m_r[b+k]=m_func[b+k]-m_cX*(m_data[b+k-sX]+m_data[b+k+sX]-2.0*m_data[b+k])
             -m_cY*(m_data[b+k-m_sizeZ]+m_data[b+k+m_sizeZ]-2.0*m_data[b+k])
             -m_cZ*(m_data[b+k-1]+m_data[b+k+1]-2.0*m_data[b+k]);
    if (k == (m_sizeZ - 2)) {
      if (j == (m_sizeY - 2)) {
         b+=2*m_sizeZ;
         j = 0;
         i++;
     }
     k = 0;
     b+=m_sizeZ;
     j++;
    }
}

answered May 29 '14 at 09:15

Dr. Debasish Jana

6,980
4
30
69

Looks like it's better as I don't get the bad loop form message, however I get this one : not vectorized: not enough data-refs in basic block. – repptilia May 29 '14 at 09:20
This is likely not a good idea because you introduce additional `if()` statements into the inner loop. CPUs don't like those. – cmaster - reinstate monica May 29 '14 at 09:42

C++ Optimize vectorizing nested loops

4 Answers4