Type stability

[milankl]

Counting the float operations, we currently aren't as type-stable as I was expecting

julia> using GenericFFT, GFlops
julia> rfft_plan = plan_rfft(zeros(Float16,1024))
GenericFFT.DummyrFFTPlan{ComplexF16, false, UnitRange{Int64}}(1024, 1:1, #undef)

julia> @count_ops rfft_plan*randn(Float16,1024)
Flop Counter: 61382 flop
┌────────┬─────────┬─────────┬─────────┐
│        │ Float16 │ Float32 │ Float64 │
├────────┼─────────┼─────────┼─────────┤
│    fma │       0 │       0 │       4 │
│ muladd │       0 │       0 │     153 │
│    add │   18449 │       0 │      83 │
│    sub │   16447 │       0 │      39 │
│    mul │   24622 │       0 │    1250 │
│    div │      40 │      23 │       2 │
│    abs │      40 │      19 │      17 │
│    neg │      15 │       0 │       3 │
│   sqrt │       0 │      19 │       0 │
└────────┴─────────┴─────────┴─────────┘

The Float32 operations might be miscounted similar to triscale-innov/GFlops.jl#40 but I doubt the Float64 operations are. Just raising this as we may want to ensure type stability and explicit conversions rather than relying on promotions (which can easily cascade into Float64s where we don't actually want to use them). The output may still be of eltype T but users of this package probably want a Fourier transform fully in T when they provide an input vector of eltype T?

[daanhb]

Heh, I did not know about @count_ops. I'd say you're right that it is fair to expect a transform in T if T is provided (up to the Float16 vs Float32 issue).

[milankl]

The culprit might be the sinpi function. Mimicking generic_fft_pow2!

julia> function sinpi_test(::Type{T},n::Integer) where T
           big2 = 2one(T)
           logn = 2
           while logn < n
               θ = -big2/logn
               wtemp = sinpi(θ/2)
               wpi = sinpi(θ)
               logn <<= 1
           end
       end
sinpi_test (generic function with 1 method)

julia> @count_ops sinpi_test(Float16,1024)
Flop Counter: 431 flop
┌────────┬─────────┬─────────┬─────────┐
│        │ Float16 │ Float32 │ Float64 │
├────────┼─────────┼─────────┼─────────┤
│ muladd │       0 │       0 │      28 │
│    add │      18 │       0 │      32 │
│    sub │      58 │       0 │       0 │
│    mul │      36 │       0 │      90 │
│    div │      36 │      21 │       0 │
│    abs │      36 │      17 │       0 │
│    neg │      14 │       0 │       0 │
│   sqrt │       0 │      17 │       0 │
└────────┴─────────┴─────────┴─────────┘

So maybe it's okay because the call of sinpi does not depend on the input vector x only on its length? Are these sinpi factors usually something that would be contained in a plan?

[milankl]

As far as I understand this, for an $2^n$-length vector only $n$ calls to sinpi are actually needed. generic_fft_pow2! currently performs $2n$.

[daanhb]

Are the sinpi calls expensive? It seems like reducing the number of calls requires some analytical thinking here about the algorithm. Overall, anything could be put into a plan, as long as the plan itself does not become too large. FFTW precomputes lots of things. To get an idea, see their Documentation, for example The design and implementation of FFTW3 (link to pdf).

On topic, if the issue arises from the sinpi calls, then the current code of GenericFFT actually is type-stable?