# Possible bug using plan\_irfft and plan\_rfft with mul! and ldiv!

**URL:** <https://discourse.julialang.org/t/possible-bug-using-plan-irfft-and-plan-rfft-with-mul-and-ldiv/57374>\
**Category:** General Usage\
**Tags:** question, fftw, bug, linearalgebra\
**Created:** [March 17, 2021, 12:02pm UTC](https://discourse.julialang.org/t/possible-bug-using-plan-irfft-and-plan-rfft-with-mul-and-ldiv/57374 "2021-03-17T12:02:58Z")\
**Posts on this page:** 3\
**Page:** 1

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [March 17, 2021, 12:02pm UTC](https://discourse.julialang.org/t/possible-bug-using-plan-irfft-and-plan-rfft-with-mul-and-ldiv/57374/1 "2021-03-17T12:02:58Z")

</div>

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:

```julia
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
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
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!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [March 17, 2021, 12:11pm UTC](https://discourse.julialang.org/t/possible-bug-using-plan-irfft-and-plan-rfft-with-mul-and-ldiv/57374/2 "2021-03-17T12:11:18Z")

</div>

> [@tiago](#):
>
> For example, when computing the inverse transformation using the inverse plan with `mul!` , the two input matrices are mutated.

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](http://fftw.org/fftw3_doc/Planner-Flags.html#index-FFTW_005fPRESERVE_005fINPUT-1) in the multidimensional case.

---

<div class="post-metadata">

**Author:** ![tiago](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tiago/32/22567_2.png) [@tiago](https://discourse.julialang.org/u/tiago)\
**Post date:** [March 17, 2021, 1:54pm UTC](https://discourse.julialang.org/t/possible-bug-using-plan-irfft-and-plan-rfft-with-mul-and-ldiv/57374/3 "2021-03-17T13:54:32Z")

</div>

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](https://global.discourse-cdn.com/julialang/original/3X/8/2/829479642890d60f06c63b6736c0c4300ceb61af.png)

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!