Possible bug using plan_irfft and plan_rfft with mul! and ldiv!

General Usage
, , ,
1

Hello,
I found a strange behaviour while using precomputed plans in 2D of real fft, plan_rfft and plan_irfft, with the function mul! and ldiv! from the LinearAlbegra package. For example, when computing the inverse transformation using the inverse plan with mul!, the two input matrices are mutated.
Let me show you:

using FFTW, LinearAlgebra
N  = 4
hN = N÷2+1
p_2d  = plan_rfft(zeros(N,N))
pi_2d = plan_irfft(zeros(ComplexF64,hN,N),N)

B = rand(N,N)
b = Matrix{ComplexF64}(undef,hN,N)

mul!(b,p_2d,B) # compute the fourier coefficients

giving the output:

julia> b
3×4 Array{Complex{Float64},2}:
   9.03435+0.0im        0.928178+0.0428595im   -1.26097+0.0im        0.928178-0.0428595im
 -0.169974+0.736463im   -1.19741-1.28333im    -0.491232+0.206306im   -1.59964+0.561067im
  -1.08496+0.0im       -0.511485-1.51638im     0.944023+0.0im       -0.511485+1.51638im

Now, when we compute the inverse using mul!(B,pi_2d,b) both matrices, B and b are mutated, but this is not expected, only B should be mutated. As you can see:

julia> mul!(B,pi_2d,b);
julia> b
3×4 Array{Complex{Float64},2}:
  9.62973+0.0im      10.2096+0.0im        5.91702+0.0im        10.381+0.0im
 -3.45825+0.22051im  2.16565+0.932386im   2.13584+1.66503im  -1.52313+0.127927im
 -1.16391+0.0im      1.00377+0.0im       0.882031+0.0im      -5.06174+0.0im

Whereas B is correctly computed, b is also changed. This also happens when using ldiv!(B,p_2d,b) to compute the inverse.

Since this behaviour doesn’t happen when using real plans in 1D or complex plans in 1D and 2D, I guess we could say there is a bug using real plans in 2D?

I didn’t know where to post this, If there is a better place to write this post, somebody can point it out?

Thanks!

1 Like
2

This is not a bug, it’s a limitation of the underlying FFTW library: inverse multidimensional rfft plans mutate the input.

You can prevent this by creating the inverse (c2r) plan with the FFTW.PRESERVE_INPUT flag, at the expense of some performance. Whoops, no, this flag is not supported in the multidimensional case.

3 Likes
3

Oh ok, so reading the documentation the flag is also not supported when using the direct plan to compute the inverse, i.e. ldiv!(B,p_2d,b). Actually throws an error:

image
image955×217 13.9 KB

Thanks so much for the time :+1:, I’ll compute the inverse transforms using B=pi_2d*b, this way the input is not destroyed, with little sacrifice of performance. Thanks once again!