Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Inplace irfft calculations with mul! mutate the fourier-space array #101

Open
VanillaBrooks opened this issue May 12, 2023 · 0 comments
Open

Comments

@VanillaBrooks
Copy link

VanillaBrooks commented May 12, 2023

I ran into some unexpected behavior with both FFTW.jl and FastTransforms.jl but not CUFFT.jl. I am thinking this repo is the appropriate location to raise this issue, but its possible it is an issue with both FFTW.jl and FastTransforms.jl

Summary

When applying a irfft plan to a Fourier space matrix uhat with mul! and a preallocated output array u, uhat is unexpectedly mutated.

Given a uhat matrix of elements ComplexF64 and shape (3,4), these work as expected:

u = irfft(uhat, 4, 1:2)
plan = plan_irfft(uhat, 4, 1:2)
u = plan * uhat

this generates the correct u, but mutates uhat:

plan = plan_irfft(uhat, 4, 1:2)
u = zeros(4,4)
mul!(u, plan, uhat)

Reproduction

using FFTW
using LinearAlgebra: mul!

# (3,4) array
uhat = ComplexF64.([
    1 2 3 4;
    1 2 3 4;
    1 2 3 4;
])

#
# first, perform an irfft without a plan
#
u = irfft(uhat, 4, 1:2)

println("\n\n########\n########")
println("unplanned irfft:")
display(u)
println("û: unchanged:")
display(uhat)

# generate an irfft plan
plan = plan_irfft(uhat, 4, 1:2)

#
# calculate the irfft with the plan (allocates)
#
u_plan = plan * uhat

println("\n\n########\n########")
println("planned irfft with allocation:")
display(u_plan)
println("uhat: unchanged:")
display(uhat)

#
# calculate the irfft with the plan (reuse allocation)
#

u_plan_preallocation = zeros(4,4)
mul!(u_plan_preallocation, plan, uhat)

println("\n\n########\n########")
println("planned irfft with preallocation:")
display(u_plan)
println("uhat: HAS changed:")
display(uhat)

This outputs:

########
########
unplanned irfft:
4×4 Matrix{Float64}:
 2.5  -0.5   -0.5  -0.5
 0.0   0.25   0.0  -0.25
 0.0   0.0    0.0   0.0
 0.0  -0.25   0.0   0.25
û: unchanged:
3×4 Matrix{ComplexF64}:
 1.0+0.0im  2.0+0.0im  3.0+0.0im  4.0+0.0im
 1.0+0.0im  2.0+0.0im  3.0+0.0im  4.0+0.0im
 1.0+0.0im  2.0+0.0im  3.0+0.0im  4.0+0.0im


########
########
planned irfft with allocation:
4×4 Matrix{Float64}:
 2.5  -0.5   -0.5  -0.5
 0.0   0.25   0.0  -0.25
 0.0   0.0    0.0   0.0
 0.0  -0.25   0.0   0.25
uhat: unchanged:
3×4 Matrix{ComplexF64}:
 1.0+0.0im  2.0+0.0im  3.0+0.0im  4.0+0.0im
 1.0+0.0im  2.0+0.0im  3.0+0.0im  4.0+0.0im
 1.0+0.0im  2.0+0.0im  3.0+0.0im  4.0+0.0im


########
########
planned irfft with preallocation:
4×4 Matrix{Float64}:
 2.5  -0.5   -0.5  -0.5
 0.0   0.25   0.0  -0.25
 0.0   0.0    0.0   0.0
 0.0  -0.25   0.0   0.25
uhat: HAS changed:
3×4 Matrix{ComplexF64}:
 10.0+0.0im  -2.0-2.0im  -2.0+0.0im  -2.0+2.0im
 10.0+0.0im  -2.0-2.0im  -2.0+0.0im  -2.0+2.0im
 10.0+0.0im  -2.0-2.0im  -2.0+0.0im  -2.0+2.0im
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Labels
None yet
Projects
None yet
Development

No branches or pull requests

1 participant