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

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

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.

2 Likes

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:

Thanks so much for the time , 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!