start with fmm_leg2cheb

Author: MikaelSlevinsky

Project.toml

@@ -21,7 +21,7 @@ AbstractFFTs = "1.0"
 BandedMatrices = "1.5"
 FFTW = "1.7"
 FastGaussQuadrature = "0.4, 0.5, 1"
-FastTransforms_jll = "0.6.2"
+FastTransforms_jll = "0.6.3"
 FillArrays = "0.9, 0.10, 0.11, 0.12, 0.13, 1"
 GenericFFT = "0.1"
 Reexport = "0.2, 1.0"

src/FastTransforms.jl

@@ -117,13 +117,12 @@ include("toeplitzhankel.jl")
 include("SymmetricToeplitzPlusHankel.jl")
 
 # following use libfasttransforms by default
-for f in (:jac2jac,
-    :lag2lag, :jac2ultra, :ultra2jac, :jac2cheb,
-    :cheb2jac, :ultra2cheb, :cheb2ultra, :associatedjac2jac,
-    :modifiedjac2jac, :modifiedlag2lag, :modifiedherm2herm,
-    :sph2fourier, :sphv2fourier, :disk2cxf, :ann2cxf,
-    :rectdisk2cheb, :tri2cheb, :tet2cheb,
-    :leg2cheb, :cheb2leg, :ultra2ultra)
+for f in (:leg2cheb, :cheb2leg, :ultra2ultra, :jac2jac,
+          :lag2lag, :jac2ultra, :ultra2jac, :jac2cheb,
+          :cheb2jac, :ultra2cheb, :cheb2ultra, :associatedjac2jac,
+          :modifiedjac2jac, :modifiedlag2lag, :modifiedherm2herm,
+          :sph2fourier, :sphv2fourier, :disk2cxf, :ann2cxf,
+          :rectdisk2cheb, :tri2cheb, :tet2cheb, :fmm_leg2cheb)
     lib_f = Symbol("lib_", f)
     @eval $f(x::AbstractArray, y...; z...) = $lib_f(x, y...; z...)
 end

src/libfasttransforms.jl

@@ -144,6 +144,7 @@ end
     SPINSPHERESYNTHESIS
     SPINSPHEREANALYSIS
     SPHERICALISOMETRY
+    FMMLEG2CHEB
 end
 
 Transforms(t::Transforms) = t
@@ -187,7 +188,8 @@ let k2s = Dict(LEG2CHEB             => "Legendre--Chebyshev",
                TETRAHEDRONANALYSIS  => "FFTW Chebyshev analysis on the tetrahedron",
                SPINSPHERESYNTHESIS  => "FFTW Fourier synthesis on the sphere (spin-weighted)",
                SPINSPHEREANALYSIS   => "FFTW Fourier analysis on the sphere (spin-weighted)",
-               SPHERICALISOMETRY    => "Spherical isometry")
+               SPHERICALISOMETRY    => "Spherical isometry",
+               FMMLEG2CHEB          => "FMM Legendre--Chebyshev")
     global kind2string
     kind2string(k::Union{Integer, Transforms}) = k2s[Transforms(k)]
 end
@@ -314,6 +316,50 @@ destroy_plan(p::FTPlan{Complex{Float64}, 2, SPINSPHERESYNTHESIS}) = ccall((:ft_d
 destroy_plan(p::FTPlan{Complex{Float64}, 2, SPINSPHEREANALYSIS}) = ccall((:ft_destroy_spinsphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 destroy_plan(p::FTPlan{Float64, 2, SPHERICALISOMETRY}) = ccall((:ft_destroy_sph_isometry_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
 
+destroy_plan(p::FTPlan{Float32, 1, FMMLEG2CHEB}) = ccall((:ft_free_fmmf, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
+destroy_plan(p::FTPlan{Float64, 1, FMMLEG2CHEB}) = ccall((:ft_free_fmm, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
+
+function plan_fmm_leg2cheb(::Type{Float64}, n::Integer)
+    maxs = 64
+    M = 18
+    BOTH = 2
+    lagrange = 0
+    verbose = 2
+    plan = ccall((:ft_create_fmm, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint, Cint, Cint, Cint), n, maxs, M, BOTH, lagrange, verbose)
+    return FTPlan{Float64, 1, FMMLEG2CHEB}(plan, n)
+end
+
+function mul!(y::StridedVector{Float64}, p::FTPlan{Float64, 1, FMMLEG2CHEB}, x::StridedVector{Float64})
+    checksize(p, x)
+    checksize(p, y)
+    checkstride(p, x)
+    checkstride(p, y)
+    L2C = 0
+    flops = ccall((:ft_execute, libfasttransforms), Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{ft_plan_struct}, Cint, Cint), x, y, p, L2C, 1)
+    return y
+end
+function div!(y::StridedVector{Float64}, p::FTPlan{Float64, 1, FMMLEG2CHEB}, x::StridedVector{Float64})
+    checksize(p, x)
+    checksize(p, y)
+    checkstride(p, x)
+    checkstride(p, y)
+    C2L = 1
+    flops = ccall((:ft_execute, libfasttransforms), Cint, (Ptr{Float64}, Ptr{Float64}, Ptr{ft_plan_struct}, Cint, Cint), x, y, p, C2L, 1)
+    return y
+end
+
+function *(p::FTPlan{T, 1, FMMLEG2CHEB}, x::AbstractArray{T}) where T
+    Ax = Array(x)
+    y = zero(Ax)
+    mul!(y, p, Ax)
+end
+function \(p::FTPlan{T, 1, FMMLEG2CHEB}, x::AbstractArray{T}) where T
+    Ax = Array(x)
+    y = zero(Ax)
+    div!(y, p, Ax)
+end
+
+
 struct AdjointFTPlan{T, S, R}
     parent::S
     adjoint::R
@@ -433,7 +479,7 @@ for f in (:leg2cheb, :cheb2leg, :ultra2ultra, :jac2jac,
           :cheb2jac, :ultra2cheb, :cheb2ultra, :associatedjac2jac,
           :modifiedjac2jac, :modifiedlag2lag, :modifiedherm2herm,
           :sph2fourier, :sphv2fourier, :disk2cxf, :ann2cxf,
-          :rectdisk2cheb, :tri2cheb, :tet2cheb)
+          :rectdisk2cheb, :tri2cheb, :tet2cheb, :fmm_leg2cheb)
     plan_f = Symbol("plan_", f)
     lib_f = Symbol("lib_", f)
     @eval begin