Mul! modifying input arrays when using FFTW.jl to compute 2D inverse RFFT

General Usage
,
1

I have just encountered what I believe is a bug when using FFTW.jl to compute inverse 2D FFT. As reproduced in the MWE below, the use of mul!(A,P,B) modifies the input array B.

using FFTW, LinearAlgebra;
N = 1;
B = rand(ComplexF64,(N+1,2*N));
P = plan_irfft(B,2*N);
A = zeros(Float64,(2*N,2*N));
#####
v0 = B[1,1];
mul!(A,P,B);
v1 = B[1,1];
#####
println(v0 == v1)  # --> Returns false

Playing around with these functions, the bug seems to go away if (i) the FFT is 1D rather than 2D; or if (ii) one uses the forward transform plan_rfft rather than the inverse one.

Unfortunately, I am not versed enough in julia to dig further the origin of this behaviour. Maybe someone else would be more lucky?

If useful:

Julia Version 1.11.1
Commit 8f5b7ca12ad (2024-10-16 10:53 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 12 × Apple M2 Pro
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m2)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)
1 Like
2

That is by design:

From the fftw manual:

FFTW includes multi-dimensional transforms for real data of any rank. As with the one-dimensional real transforms, they save roughly a factor of two in time and storage over complex transforms of the same size. Also as in one dimension, these gains come at the expense of some increase in complexity–the output format is different from the one-dimensional RFFTW (and is more similar to that of the complex FFTW) and the inverse (complex to real) transforms have the side-effect of overwriting their input data.

5 Likes
3

Well, a prime example of why I should have Read The Fine Manual more carefully :slight_smile:
Thank you very much, you truly have saved me a lot of trouble!