FFTW.jl about twice as slow as C++ for sizes 17, 19, 23, 29, 31, 34, 37, etc

[Marc-3d]

Hi, I am measuring the speed of cross-correlation in the frequency domain and I am comparing different packages, ( my code in Julia 1.3.1, some C++ code and some Python code ).

This is the code for generating the in-place plans, in case its useful.

    corr_size   = (IA, IA);                                 # input/output size
    corr_matrix = zeros( Float64, corr_size );              # output array
    padi        =  rand( Complex{Float64}, corr_size );     # input array 1
    pads        =  rand( Complex{Float64}, corr_size );     # input array 2
    plan        =  FFTW.plan_fft!( padi, flags=FFTW.MEASURE ); 
    iplan       = FFTW.plan_ifft!( padi, flags=FFTW.MEASURE );

After performing the evaluations for different values of IA I realized that Julia matches the speed of calling libfftw3 from C++, except for non-prime sizes, like 17, 19, 23. Curiously, Julia takes twice the time as C++, which suggests to me that maybe the code behind FFTW.jl has some internal overhead in these cases.

I decided to do show( plan ); in Julia and the respective fftw_print_plan in C++ and confirmed that the plans generated from Julia are identical as the ones from C++. This is what a plan looks like for IA = 10:

FFTW in-place forward plan for 10×10 array of Complex{Float64}
(dft-rank>=2/1
  (dft-direct-10-x10 "n1fv_10_avx2")
  (dft-direct-10-x10 "n1fv_10_avx2"))

For IA = 17:

FFTW in-place forward plan for 17×17 array of Complex{Float64}
(dft-rank>=2/1
  (dft-vrank>=1-x17/1
    (dft-generic-17))
  (dft-vrank>=1-x17/1
    (dft-generic-17)))

The plan looks different and seems to be dividing the Fourier Transform into two elements. ...two elements in the plan description, double the execution time... maybe this is hinting to something, but I don't have the knowledge to follow it.

This is basically all info I have gathered.

I can also say that I am running the code on linux, and FFTW_jl has chosen "x86_64-linux-gnu.jl" as its preferred platform. I tried linking my C++ script to the libfftw3.so in .julia/Artifacts and C++ didn't show this delay. Therefore, either I linked it wrong, or it's not a problem of the library.

[stevengj]

How are you timing the plan execution in Julia? Should be something like:

using FFTW, BenchmarkTools
x = zeros(ComplexF64, 17,17);
p = plan_fft!(x, flags=FFTW.MEASURE)
@btime $p * $x;

The plan looks different and seems to be dividing the Fourier Transform into two elements. ...two elements in the plan description, double the execution time

It's just a loop (dft-vrank>=1-x17) over 17 transforms of size 17 (dft-generic-17) for the columns and another for the rows. For the 10x10 case the loop is combined with the pregenerated kernel of size 10, but for size-17 there is no optimized kernel ("codelet") in the standard FFTW build and so it uses a generic routine.

[Marc-3d]

When I run your code I get around 6 microseconds.

I am timing with @timedwithin a loop, and then I take the average of the execution time. When I time your code with this method I get the same result, around 6 microseconds.

From C++ I get 3 microseconds.

Here I copy some of the results I got (in microseconds), so that you can appreciate the pattern:

size            Julia               C++
10             0.6910              0.509     
11             0.9465              1.007
12             0.8708              0.931
13             1.3580              1.032
14             1.2426              1.092
15             1.6779              1.128
16             1.5830              1.002
17            21.2691             10.747
18             6.3690              7.145
19            25.2877             14.888
20             2.5464              2.041
21             9.6394             11.383
22             9.4058             10.693
23            43.3771             25.654
24             8.9971              9.425
25             6.9862              6.292
26            13.4980             12.527
27            15.6364             17.585
28            11.9178             12.464
29            94.4101             48.766

Notice that in this list, for size 17 Julia takes 21.2691 microseconds, not the 6 microseconds I mentioned above when running @btime $p * $x. This is because I am measuring cross-correlation, which performs 3 Fourier Transforms. So the numbers seem to make sense, even with my way of measuring execution time.

[stevengj]

Just in case it is a difference in how you compiled FFTW, can you try linking your C++ code to the version of FFTW that Julia is using? It can be found at FFTW.FFTW_jll.libfftw3_path

[Marc-3d]

I have tried linking the C++ code to FFTW's artifact library, but the C++ code doesn't show the delay in non-prime numbers. To be explicity, I replaced LIBS = -lfftw3 from the C++ Makefile to LIBS = /home/BrucePBon/.julia/artifacts/1f7802d864a42a8958dafb40427722f8a7abbe07/lib/libfftw3.so, which is the directory returned by FFTW.FFTW_jll.libfftw3_path

I also tried to change FFTW_jll's library path, by settings global artifact_dir = "/usr/lib/x86_64-linux-gnu/" in ~/.julia/packages/FFTW_jll/version/src/wrappers/x86_64-linux-gnu.jl (line 41). After opening a new julia terminal, running using FFTW and compiling the package without errors, FFTW.jl still shows the delay in non-prime numbers.

[stevengj]

Can you post your C++ code in a gist?

[Marc-3d]

Sure, I will actually write a minimal version because I am evaluating code from a package with a GUI, Makefile and stuff, and there are many non-essential files and functionalities. I'll get back to you

[Marc-3d]

okay, here is the core of what is going on in C++. I basically get the same results as my last C++ evaluation, so it seems I didn't break anything.

Compilation is quite simple: g++ main.cpp -lfftw3

If you see any obvious mistakes please let me know. I am mostly self-taught, so I might get the stuff people learn at computer science 101 wrong.

[Marc-3d]

Also, I share the julia code. I just wrote it manually from my original file, simplifying some names. There shouldn't be any typos though.

using FFTW

function crossCorrelationFFT!( in1::Array{ComplexF64,2}, in2::Array{ComplexF64,2}, corr::Array{ComplexF64,2}, plan, iplan )
    
    plan * in1
    plan * in2
    
    for e in 1:length(in1)
        @inbounds in1[e] = conj( in1[e] ) * in2[e]; 
    end
    
    iplan * in1; 

    for e in 1:length(in1) 
        @inbounds corr[e] = real( in1[e] ) 
    end
end

IA0, IAmax, repeats = 10, 60, 1000; 

for IA in IA0:IAmax
   t     = 0.0;
   corr  = zeros( Float64, IA, IA )
   in1   = zeros( ComplexF64, IA, IA )
   in2   = zeros( ComplexF64, IA, IA )
   plan  = FFTW.plan_fft!( in1, flags=FFTW_MEASURE ); 
   iplan = FFTW.plan_ifft!( in2, flags=FFTW_MEASURE ); 

   for r in repeats
       temp = @timed crossCorrelationFFT!( in1, in2, corr, plan, iplan )
       t += temp[2]
   end

    println( IA, ": ", 1000000*t/repeats ) # 1000000* to convert from seconds to microseconds
end

[stevengj]

(Usually benchmarking in Julia is done with the BenchmarkTools package and either @btime or @belapsed. If you are benchmarking in global scope you should also generally use $ to interpolate global variables, as in my example above, so that you don't time dynamic dispatch. i.e. temp = @belapsed crossCorrelationFFT!( $in1, $in2, $corr, $plan, $iplan ). Note that BenchmarkTools runs things several times and collects statistics for you, so you don't need your repeats loop if you use @belapsed.)

FFTW.plan_ifft! creates a plan that includes a 1/N scale factor, which isn't present in your C++ code. The analogue of a FFTW_BACKWARD plan (unscaled inverse) is created by FFTW.plan_bfft!.

[Marc-3d]

I totally oversaw the normalization of ifft. I will read into BenchmarkTools cause I really want to understand whats going on behind it, and how it differs from taking an average.

Anyway, the small execution time spikes at non-prime sizes remain, even after changing FFTW.plan_ifft! by FFTW.plan_bfft! and using @belapsed to measure execution time. Have you tried executing the C++ or Julia code?

[giordano]

and how it differs from taking an average.

@belapsed and @btime in BenchmarkTools.jl give you the minimum time, not the average, which would also be misleading because of a few large outliers. You may be interested in reading their e-print, too: https://arxiv.org/abs/1608.04295