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!