Skip to content

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

@VanillaBrooks

Description

@VanillaBrooks

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

Metadata

Metadata

Assignees

No one assigned

    Labels

    No labels
    No labels

    Type

    No type
    No fields configured for issues without a type.

    Projects

    No projects

    Milestone

    No milestone

    Relationships

    None yet

    Development

    No branches or pull requests

    Issue actions