Hi, I just started to use Julia recently, and when I try to use mul! and plan_fft to get better performance, I always get a bug that I can’t deal with
My input is

u0=zeros(nz,nr)
for i in 1:nz
for j in 1:nr
R=sqrt((6*((i*h)-1))^2+((j-0.5)*h)^2)
u0[i,j]=tanh((1.2-R)/(sqrt(2)*e))
end
end
FFT=plan_fft(u0)
up=similar(u0);
mul!(up,FFT,u0);

It lacks some parameters to be sure but can you try with setting up as Complex ? The output of the FFT is indeed a complex vector and you initialise a Float64.

To apply P to an array A, use P * A; in general, the syntax for applying plans is much like that of matrices. (A plan can only be applied to arrays of the same size as the A for which the plan was created.) You can also apply a plan with a preallocated output array Â by calling mul!(Â, plan, A). (For mul!, however, the input array A must be a complex floating-point array like the output Â.) You can compute the inverse-transform plan by inv(P) and apply the inverse plan with P \ Â (the inverse plan is cached and reused for subsequent calls to inv or ), and apply the inverse plan to a pre-allocated output array A with ldiv!(A, P, Â).

Note that for a real-valued input array, you may be able to achieve better performance with rfft/plan_rfft, since a FFT on real input produces Hermitian-symmetric output with redundant negative frequency terms.

Thank you for your advice, after turning the input array into a complex array, everything seems to work fine. I also realized that rift is much faster than fft, but since I need to get full wave numbers for my simulation, I have to stick on fft this time.