Why can GCC not vectorize this function and loop?

Question

I'm attempting to make a function SIMD-enabled and vectorize the loop with a function call.

#include <cmath>

#pragma omp declare simd
double BlackBoxFunction(const double x) {
    return 1.0/sqrt(x);
}

double ComputeIntegral(const int n, const double a, const double b) {
    const double dx = (b - a)/n;
    double I = 0.0;
    #pragma omp simd reduction(+: I)

    for (int i = 0; i < n; i++) {
      const double xip12 = a + dx*(double(i) + 0.5);
      const double yip12 = BlackBoxFunction(xip12);
      const double dI = yip12*dx;
      I += dI; 
  }
  return I;
}

For the code above, if I compile it with icpc:

icpc worker.cc -qopenmp -qopt-report=5 -c

The opt-report shows that the function and loop are both vectorized. However, if I try to compile it with g++ 6.5:

g++ worker.cc -O3 -fopenmp -fopt-info-vec-missed -funsafe-math-optimizations -c

The output shows note:not vectorized: control flow in loop. and note: bad loop form, and the loop cannot be vectorized.

How can I vectorize the loop with GCC?

EDIT :

If I write the function into a separate file,

worker.cc:

#include "library.h"

double ComputeIntegral(const int n, const double a, const double b) {
    const double dx = (b - a)/n;
    double I = 0.0;
    #pragma omp simd reduction(+: I)

    for (int i = 0; i < n; i++) {
      const double xip12 = a + dx*(double(i) + 0.5);
      const double yip12 = BlackBoxFunction(xip12);
      const double dI = yip12*dx;
      I += dI; 
  }
  return I;
}

library.h:

#ifndef __INCLUDED_LIBRARY_H__
#define __INCLUDED_LIBRARY_H__

#pragma omp declare simd
double BlackBoxFunction(const double x); 

#endif

and library.cc:

#include <cmath>

#pragma omp declare simd
double BlackBoxFunction(const double x) {
  return 1.0/sqrt(x);
}

Then I compile it with GCC:

g++ worker.cc library.cc -O3 -fopenmp -fopt-info-vec-missed -funsafe-math-optimizations -c

It shows:

worker.cc:9:31: note: loop vectorized

but

library.cc:5:18: note:not vectorized: control flow in loop.
library.cc:5:18: note:bad loop form.

It makes me confused. I wonder whether it is already vectorized.

Relevant GCC source. You can see that it raises the control flow warning if there are more than two basic blocks inside the loop, counting the loop machinery as one of them. I'd guess that means the function call is splitting the loop contents into two or three blocks; if you manually inline BlackBoxFunction does it work? — Rup
– Rup, Commented Jan 11, 2019 at 12:51
Possible duplicate of GCC: vectorization difference between two similar loops — Damian
– Damian, Commented Jan 11, 2019 at 12:53
@Rup I tried but GCC still gives the hint that control flow in loop and bad loop form — pangbryant
– pangbryant, Commented Jan 11, 2019 at 13:11

wim · Accepted Answer · 2019-01-11 14:10:58Z

7

Vectorization is possible with gcc, after some slight modifications of the code:

#include <cmath>

double BlackBoxFunction(const double x) {
    return 1.0/sqrt(x);
}

double ComputeIntegral(const int n, const double a, const double b) {
    const double dx = (b - a)/n;
    double I = 0.0;
    double d_i = 0.0;
    for (int i = 0; i < n; i++) {
      const double xip12 = a + dx*(d_i + 0.5);
      d_i = d_i + 1.0;
      const double yip12 = BlackBoxFunction(xip12);
      const double dI = yip12*dx;
      I += dI; 
  }
  return I;
}

This was compiled with the compiler options: -Ofast -march=haswell -fopt-info-vec-missed -funsafe-math-optimizations. The main loop compiles to

.L7:
    vaddpd  ymm2, ymm4, ymm7
    inc     eax
    vaddpd  ymm4, ymm4, ymm8
    vfmadd132pd     ymm2, ymm9, ymm5
    vsqrtpd ymm2, ymm2
    vdivpd  ymm2, ymm6, ymm2
    vfmadd231pd     ymm3, ymm5, ymm2
    cmp     eax, edx
    jne     .L7

See the following Godbolt link

I removed the #pragma omp ..., because they didn't improve the vectorization, but they did not made the vectorization worse either.

Note that only changing the compiler option from -O3 to -Ofast is sufficient to enable vectorization. Nevertheless, it is more efficient to use a double counter than an int counter which is converted to double each iteration.

Note also that the vectorization reports are quite misleading. Inspect the generated assembly code to see whether or not the vectorization was successful.

edited Jan 11, 2019 at 14:10

answered Jan 11, 2019 at 13:41

wim

4,10823 silver badges28 bronze badges

Sign up to request clarification or add additional context in comments.

11 Comments

wim Over a year ago

You can also consider starting with double xip12 = a + 0.5*dx; outside the loop, and increment with xip12 = xip12 + dx;, which would be slightly more efficient.

wim Over a year ago

Option -Ofast enables both -funsafe-math-optimizations and -fno-math-errno. In order to vectorize the sqrt, you need -fno-math-errno. In the vectorized case, the additions are ordered in a different way than in the scalar case. Therefore you also need the -funsafe-math-optimizations, because floating point additions are not associative. With a more simple BlackBoxFunction, the code already vectorizes with only -O2 -funsafe-math-optimizations -ftree-vectorize, see this Godbolt link .

Peter Cordes Over a year ago

@pangbryant: worth mentioning that ICPC enables fast-math by default, somewhat similar to gcc -O3 -ffast-math. OpenMP pragmas can allow auto-vectorization of FP reductions without -ffast-math (just -O3), but not in this case. Perhaps if you only used -fno-math-errno to allow sqrt to fully inline, gcc would be able to auto-vectorize the reduction thanks to the pragma enabling optimizations that change the order of summing the result. (FP math is not strictly associative because of different rounding from different summation orders.)

Peter Cordes Over a year ago

@Zboson: KNL is the exception to "most", but yes look at Agner Fog's insn table: SKX's divps throughput is scalar/xmm: 3c, ymm: 5c, zmm: 10c. But xmm and ymm both have 11c latency, with ZMM having 18c latency. The double numbers for [v]divpd throughput are 4, 8, and 16c throughput. So it's normal for the dividers to be about half width, not fully pipelined. On KNL, SSE xmm is 32c latency / 20c throughput, while AVX-anything is apparently 32/32 regardless of float/double. Xeon Phi has AVX512ER because divide throughput is garbage.

wim Over a year ago

@Zboson: Peter's comment is illustrated by this Godbolt link. The loop vectorizes with #pragma omp simd reduction(+: I), but not with #pragma omp simd (without the reduction(+: I), unless the -fassociative-math, -fno-signed-zeros, and -fno-trapping-math options are added to gcc.

|

Collectives™ on Stack Overflow

Why can GCC not vectorize this function and loop?

1 Answer 1

11 Comments

Your Answer

Linked

Hot Network Questions

Collectives™ on Stack Overflow

1 Answer 1

11 Comments

Your Answer

Sign up or log in

Post as a guest

Linked

Related