JuliaMath / NFFT.jl

Julia implementation of the Non-equidistant Fast Fourier Transform (NFFT)
Other
153 stars 28 forks source link

How to create `Float32` plan? #99

Closed JeffFessler closed 2 years ago

JeffFessler commented 2 years ago

For large-scale problems I'd like to be able to save how much data is moved around by using Float32 types. The following code produces 64-bit output even though I tried hard to make it 32-bit:

using NFFT
M, N = 30, 20
T = Float32
w = LinRange(0,7π/8,M)
w = T.(w)
p = plan_nfft(w/(2π), N)
# p = plan_nfft(Array{T}, w/(2π), N) # does not help
M = length(w)
Tx = ComplexF32
x = zeros(Tx, N); x[N÷2+1] = 1
y = p * x
@show eltype(y), T # y is ComplexF64 !?

I could not determine from the docs or docstrings how to make the plan 32-bit. typeof(p) returns NFFTPlan{Float64, 1, 1}

I have a hunch that the culprit may be that the following lines of code seem to ignore the type T, but it is just a guess: https://github.com/JuliaMath/NFFT.jl/blob/02910a5be7bbfceb5e77631953d3f15678534f92/src/NFFT.jl#L46

This is slightly related to #33.

Is there a kwarg I'm overlooking for specifying the plan element type?

tknopp commented 2 years ago

just quick answer:

p = plan_nfft(Float32.(w/(2π)), N)

i.e. it is derived from the type of the sampling nodes. Documentation certainly needs to improved.

JeffFessler commented 2 years ago

Oh gosh do I feel dumb. I thought I took care of that with the line w = T.(w) in the code above, but I forgot that w/(2π) would convert it back to 64-bit. Doh! Thanks for the quick reply and it works fine now so I'll close the issue.