From 87530a6277ee6443bacdd908b70a04be5148f53f Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Wed, 4 Sep 2024 12:14:41 -0500 Subject: [PATCH 1/2] start with fmm_leg2cheb --- Project.toml | 2 +- src/FastTransforms.jl | 13 +++++------ src/libfasttransforms.jl | 50 ++++++++++++++++++++++++++++++++++++++-- 3 files changed, 55 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index a9269fe0..001ec41f 100644 --- a/Project.toml +++ b/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" diff --git a/src/FastTransforms.jl b/src/FastTransforms.jl index 82c34990..c8365ac3 100644 --- a/src/FastTransforms.jl +++ b/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 diff --git a/src/libfasttransforms.jl b/src/libfasttransforms.jl index ea724ce9..62bcf274 100644 --- a/src/libfasttransforms.jl +++ b/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 From a08c90a9467ac003eca7d100f0d90bc01dc3359b Mon Sep 17 00:00:00 2001 From: MikaelSlevinsky Date: Thu, 5 Sep 2024 09:44:19 -0500 Subject: [PATCH 2/2] mul! -> lmul! --- src/libfasttransforms.jl | 30 +++++++++--------------------- 1 file changed, 9 insertions(+), 21 deletions(-) diff --git a/src/libfasttransforms.jl b/src/libfasttransforms.jl index 62bcf274..8a5eb705 100644 --- a/src/libfasttransforms.jl +++ b/src/libfasttransforms.jl @@ -324,42 +324,30 @@ function plan_fmm_leg2cheb(::Type{Float64}, n::Integer) M = 18 BOTH = 2 lagrange = 0 - verbose = 2 + verbose = 0 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}) +function lmul!(p::FTPlan{Float64, 1, FMMLEG2CHEB}, x::StridedVector{Float64}) checksize(p, x) - checksize(p, y) checkstride(p, x) - checkstride(p, y) + y = zero(x) 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 + x .= y + return x end -function div!(y::StridedVector{Float64}, p::FTPlan{Float64, 1, FMMLEG2CHEB}, x::StridedVector{Float64}) +function ldiv!(p::FTPlan{Float64, 1, FMMLEG2CHEB}, x::StridedVector{Float64}) checksize(p, x) - checksize(p, y) checkstride(p, x) - checkstride(p, y) + y = zero(x) 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) + x .= y + return x end - struct AdjointFTPlan{T, S, R} parent::S adjoint::R