FFT(W) does not infer

Hello,

Is this behavior normal/expected or is it a bug (I just noticed it doesn’t happen for 2 dimensional arrays)?

fff()=(u=rand(3,3,4);FFTW.fft(u))
@code_warntype optimize=:true fff()
Variables
  #self#::Core.Compiler.Const(fff, false)
  u::Array{Float64,3}

Body::Any
1 ─ %1 = invoke Random.rand(Random.Float64::Type{Float64}, (3, 3, 4)::Tuple{Int64,Int64,Int64})::Array{Float64,3}
│   %2 = FFTW.fft::Core.Compiler.Const(AbstractFFTs.fft, false)
│   %3 = invoke %2(%1::Array{Float64,3}, $(QuoteNode(1:3))::UnitRange{Int64})::Any
└──      return %3

The code instability is affecting some code I am working on, a type assertion helps the rest infer.

I appreciate any help :slight_smile:

1 Like

Make the array complex maybe?

julia> fff()=(u=rand(3,3,4);FFTW.fft(complex.(u)))
fff (generic function with 1 method)

julia> @code_warntype fff()
Variables
  #self#::Core.Compiler.Const(fff, false)
  u::Array{Float64,3}

Body::Array{Complex{Float64},3}
1 ─      (u = Main.rand(3, 3, 4))
│   %2 = FFTW.fft::Core.Compiler.Const(AbstractFFTs.fft, false)
│   %3 = Base.broadcasted(Main.complex, u)::Base.Broadcast.Broadcasted{Base.Broadcast.DefaultArrayStyle{3},Nothing,typeof(complex),Tuple{Array{Float64,3}}}
│   %4 = Base.materialize(%3)::Array{Complex{Float64},3}
│   %5 = (%2)(%4)::Array{Complex{Float64},3}
└──      return %5

The complex conversion allocates, needlessly. Something is off, as this is not needed in 2D:

using FFTW,BenchmarkTools
fff()=(u=rand(3,4);FFTW.fft(u))
@code_warntype optimize=:true fff()

Variables
  #self#::Core.Compiler.Const(fff, false)
  u::Array{Float64,2}

Body::Array{Complex{Float64},2}
1 ─ %1  = invoke Random.rand(Random.Float64::Type{Float64}, (3, 4)::Tuple{Int64,Int64})::Array{Float64,2}
│   %2  = Base.arraysize(%1, 1)::Int64
│   %3  = Base.arraysize(%1, 2)::Int64
│   %4  = Base.slt_int(%2, 0)::Bool
│   %5  = Base.ifelse(%4, 0, %2)::Int64
│   %6  = Base.slt_int(%3, 0)::Bool
│   %7  = Base.ifelse(%6, 0, %3)::Int64
│   %8  = $(Expr(:foreigncall, :(:jl_alloc_array_2d), Array{Complex{Float64},2}, svec(Any, Int64, Int64), 0, :(:ccall), Array{Complex{Float64},2}, :(%5), :(%7), :(%7), :(%5)))::Array{Complex{Float64},2}
│   %9  = invoke circcopy!(%8::Array{Complex{Float64},2}, %1::Array{Float64,2})::Array{Complex{Float64},2}
│         Base.arraysize(%9, 1)
│         Base.arraysize(%9, 2)
│   %12 = AbstractFFTs.plan_fft::typeof(plan_fft)
│   %13 = invoke FFTW.:(var"#plan_fft#7")(FFTW.ESTIMATE::UInt32, FFTW.NO_TIMELIMIT::Float64, %12::typeof(plan_fft), %9::Array{Complex{Float64},2}, $(QuoteNode(1:2))::UnitRange{Int64})::FFTW.cFFTWPlan{Complex{Float64},-1,false,2}
│   %14 = invoke AbstractFFTs.:*(%13::FFTW.cFFTWPlan{Complex{Float64},-1,false,2}, %9::Array{Complex{Float64},2})::Array{Complex{Float64},2}
└──       return %14

This is what fft does under the hood anyway, you just don’t see it.

That I know, and it shows in the typed code in 2D. But the conversion to complex would allocate another array that is not used. Doing what FFT should be doing by hand, infers:

fff2()=(u=rand(3,4,5);fu=complex(u);fft!(fu))
@code_warntype optimize=:true fff2()
Variables
  #self#::Core.Compiler.Const(fff2, false)
  u::Array{Float64,3}
  fu::Array{Complex{Float64},3}

Body::Array{Complex{Float64},3}
1 ─ %1  = invoke Random.rand(Random.Float64::Type{Float64}, (3, 4, 5)::Tuple{Int64,Int64,Int64})::Array{Float64,3}
│   %2  = Base.arraysize(%1, 1)::Int64
│   %3  = Base.arraysize(%1, 2)::Int64
│   %4  = Base.arraysize(%1, 3)::Int64
│   %5  = Base.slt_int(%2, 0)::Bool
│   %6  = Base.ifelse(%5, 0, %2)::Int64
│   %7  = Base.slt_int(%3, 0)::Bool
│   %8  = Base.ifelse(%7, 0, %3)::Int64
│   %9  = Base.slt_int(%4, 0)::Bool
│   %10 = Base.ifelse(%9, 0, %4)::Int64
│   %11 = $(Expr(:foreigncall, :(:jl_alloc_array_3d), Array{Complex{Float64},3}, svec(Any, Int64, Int64, Int64), 0, :(:ccall), Array{Complex{Float64},3}, :(%6), :(%8), :(%10), :(%10), :(%8), :(%6)))::Array{Complex{Float64},3}
│   %12 = invoke Base.copyto!($(QuoteNode(IndexLinear()))::IndexLinear, %11::Array{Complex{Float64},3}, $(QuoteNode(IndexLinear()))::IndexLinear, %1::Array{Float64,3})::Array{Complex{Float64},3}
│         Base.arraysize(%12, 1)
│         Base.arraysize(%12, 2)
│         Base.arraysize(%12, 3)
│   %16 = FFTW.ESTIMATE::UInt32
│   %17 = FFTW.NO_TIMELIMIT::Float64
│   %18 = invoke FFTW.cFFTWPlan{Complex{Float64},-1,true,3}(%12::Array{Complex{Float64},3}, %12::Array{Complex{Float64},3}, $(QuoteNode(1:3))::UnitRange{Int64}, %16::UInt32, %17::Float64)::FFTW.cFFTWPlan{Complex{Float64},-1,true,3}
│   %19 = invoke AbstractFFTs.:*(%18::FFTW.cFFTWPlan{Complex{Float64},-1,true,3}, %12::Array{Complex{Float64},3})::Array{Complex{Float64},3}
└──       return %19

I think the problem is that AbstractFFTs.complexfloat is not type stable on a 3-d array, this is what is used to convert the real to the complex array:

julia> @code_warntype AbstractFFTs.complexfloat(u)
Variables
  #self#::Core.Compiler.Const(AbstractFFTs.complexfloat, false)
  x::Array{Float64,3}

Body::Any
1 ─ %1 = AbstractFFTs.zero($(Expr(:static_parameter, 1)))::Core.Compiler.Const(0.0, false)
│   %2 = AbstractFFTs.fftfloat(%1)::Core.Compiler.Const(0.0, false)
│   %3 = AbstractFFTs.complex(%2)::Core.Compiler.Const(0.0 + 0.0im, false)
│   %4 = AbstractFFTs.typeof(%3)::Core.Compiler.Const(Complex{Float64}, false)
│   %5 = AbstractFFTs.copy1(%4, x)::Any
└──      return %5
1 Like

Which leads all the way to Base.circcopy! for the real source of the instability, still chasing it down…

My guess is more on the

    y = Array{T}(undef, map(length, axes(x)))

line, the size of the y array, which is part of its type, can’t probably be inferred.

length isn’t part of array types. Only number of dimensions.

In fact I said size :slight_smile:

Oops. Reading is hard

I suspect the real culprit here to be circcopy!, see Circcopy! does not infer for 3D arrays.

1 Like