Plan FFTs both in julia 0.6.3 and 0.7

fftw

#1

I’m trying to make my package (FourierFlows.jl) compatible with v0.7.0. I have some trouble planning the fft’s. I can plan the ffts and do in-place ffts in both 0.7.0 and 0.6.3. However, I don’t know how to code it up in a manner that would work in both versions. I looked up Compat.jl but I couldn’t find FFTW modules included in it.

I attach here two minimal examples of what I’m trying to do. Both examples work in the corresponding julia versions.

In julia 0.6.3

export AbstractGrid

abstract type AbstractGrid end

export OneDGrid

struct OneDGrid{T} <: AbstractGrid
  nx::Int
  nk::Int
  nkr::Int
  Lx::T
  dx::T
  x::Array{T,1}
  kr::Array{T,1}
  rfftplan::FFTW.rFFTWPlan{T,-1,false,1}
  irfftplan::Base.DFT.ScaledPlan{Complex{T},Base.DFT.FFTW.rFFTWPlan{Complex{T},1,false, 1},T}
end

function OneDGrid(nx, Lx; x0=-0.5*Lx)
  T = typeof(Lx)
  dx = Lx/nx
  x = Array{T}(linspace(x0, x0+Lx-dx, nx))

  nk = nx
  nkr = Int(nx/2+1)

  kr = Array{T}(2π/Lx*cat(1, 0:Int(nx/2)))

  rfftplan  = plan_rfft(Array{T,1}(nx))
  irfftplan = plan_irfft(Array{Complex{T},1}(nkr), nx)

  OneDGrid(nx, nk, nkr, Lx, dx, x, kr, rfftplan, irfftplan)
end

n, L = 2^8, 2π
g = OneDGrid(n, L)

freal = rand(n)
fhr1 = rfft(freal)
fhr2 = deepcopy(0*fhr1)
A_mul_B!(fhr2, g.rfftplan, freal)
freal1 = irfft(fhr1, n)
freal2 = deepcopy(0*freal1)
A_mul_B!(freal2, g.irfftplan, fhr2)

println("tests for rfft")
println(vecnorm(fhr1-fhr2))
println(" ")
println("tests for irfft")
println(vecnorm(freal-freal1))
println(vecnorm(freal-freal2))

while in julia 0.7.0

using FFTW

if VERSION > v"0.6.3-"
    using LinearAlgebra
end

export AbstractGrid

abstract type AbstractGrid end

export OneDGrid

struct OneDGrid{T} <: AbstractGrid
  nx::Int
  nk::Int
  nkr::Int
  Lx::T
  dx::T
  x::Array{T,1}
  kr::Array{T,1}
  rfftplan::FFTW.rFFTWPlan{T,-1,false,1}
end

function OneDGrid(nx, Lx; x0=-0.5*Lx)
  T = typeof(Lx)
  dx = Lx/nx
  x = Array{T}(range(x0, stop=x0+Lx-dx, length=nx))

  nk = nx
  nkr = Int(nx/2+1)

  kr = Array{T}(2π/Lx*cat(0:Int(nx/2); dims = 1))

  rfftplan  = plan_rfft(Array{Float64, 1}(undef, nx))

  OneDGrid(nx, nk, nkr, Lx, dx, x, kr, rfftplan)
end

n, L = 2^8, 2π
g = OneDGrid(n, L)

freal = rand(n)
fhr1 = rfft(freal)
fhr2 = deepcopy(0*fhr1)
mul!(fhr2, g.rfftplan, freal)
freal1 = irfft(fhr1, n)
freal2 = deepcopy(0*freal1)
ldiv!(freal2, g.rfftplan, fhr2)

println("tests for rfft")
println(vecnorm(fhr1-fhr2))
println(" ")
println("tests for irfft")
println(vecnorm(freal-freal1))
println(vecnorm(freal-freal2))

#2

You are having trouble with the inverse plan, right? This is how I did it, at line 17


#3

@favba

Actually no. I did have some trouble with irfft plan in the beginning but then I figure out that the inverse plans are not needed! You can do forward fft using mul! (which is analogous to A_mul_B!) and inverse fft with ldiv!. (Have a look to the last part of my v0.7.0 code above.

My question is: How should I implement all these so that both julia v0.6.3 and v.0.7 work?