Skip to content

Commit c3da771

Browse files
MikaelSlevinskyMikaelSlevinsky
committed
start 2D adjoint plans
tests will fail unless built from source.
1 parent ef06b36 commit c3da771

File tree

1 file changed

+92
-25
lines changed

1 file changed

+92
-25
lines changed

‎src/libfasttransforms.jl

Lines changed: 92 additions & 25 deletions
Original file line numberDiff line numberDiff line change
@@ -248,8 +248,7 @@ destroy_plan(p::FTPlan{Float64, 1}) = ccall((:ft_destroy_tb_eigen_FMM, libfasttr
248248
destroy_plan(p::FTPlan{BigFloat, 1}) = ccall((:ft_mpfr_destroy_plan, libfasttransforms), Cvoid, (Ptr{mpfr_t}, Cint), p, p.n)
249249
destroy_plan(p::FTPlan{Float32, 1, ASSOCIATEDJAC2JAC}) = ccall((:ft_destroy_btb_eigen_FMMf, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
250250
destroy_plan(p::FTPlan{Float64, 1, ASSOCIATEDJAC2JAC}) = ccall((:ft_destroy_btb_eigen_FMM, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
251-
destroy_plan(p::FTPlan{Float64, 2}) = ccall((:ft_destroy_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
252-
destroy_plan(p::FTPlan{Float64, 3}) = ccall((:ft_destroy_tetrahedral_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
251+
destroy_plan(p::FTPlan{Float64}) = ccall((:ft_destroy_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
253252
destroy_plan(p::FTPlan{Complex{Float64}, 2, SPINSPHERE}) = ccall((:ft_destroy_spin_harmonic_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
254253
destroy_plan(p::FTPlan{Float64, 2, SPHERESYNTHESIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
255254
destroy_plan(p::FTPlan{Float64, 2, SPHEREANALYSIS}) = ccall((:ft_destroy_sphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
@@ -267,47 +266,63 @@ destroy_plan(p::FTPlan{Complex{Float64}, 2, SPINSPHERESYNTHESIS}) = ccall((:ft_d
267266
destroy_plan(p::FTPlan{Complex{Float64}, 2, SPINSPHEREANALYSIS}) = ccall((:ft_destroy_spinsphere_fftw_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
268267
destroy_plan(p::FTPlan{Float64, 2, SPHERICALISOMETRY}) = ccall((:ft_destroy_sph_isometry_plan, libfasttransforms), Cvoid, (Ptr{ft_plan_struct}, ), p)
269268

270-
struct AdjointFTPlan{T, S}
269+
struct AdjointFTPlan{T, S, R}
271270
parent::S
271+
adjoint::R
272+
function AdjointFTPlan{T, S, R}(parent::S) where {T, S, R}
273+
new(parent)
274+
end
275+
function AdjointFTPlan{T, S, R}(parent::S, adjoint::R) where {T, S, R}
276+
new(parent, adjoint)
277+
end
272278
end
273279

274-
AdjointFTPlan(p::FTPlan) = AdjointFTPlan{eltype(p), typeof(p)}(p)
280+
AdjointFTPlan(p::FTPlan) = AdjointFTPlan{eltype(p), typeof(p), typeof(p)}(p)
281+
AdjointFTPlan(p::FTPlan, q::FTPlan) = AdjointFTPlan{eltype(q), typeof(p), typeof(q)}(p, q)
275282

276283
adjoint(p::FTPlan) = AdjointFTPlan(p)
277284
adjoint(p::AdjointFTPlan) = p.parent
278285

279-
eltype(p::AdjointFTPlan{T, S}) where {T, S} = T
280-
ndims(p::AdjointFTPlan{T, S}) where {T, S} = ndims(p.parent)
281-
function show(io::IO, p::AdjointFTPlan{T, S}) where {T, S}
286+
eltype(p::AdjointFTPlan{T}) where T = T
287+
ndims(p::AdjointFTPlan) = ndims(p.parent)
288+
function show(io::IO, p::AdjointFTPlan)
282289
print(io, "Adjoint ")
283290
show(io, p.parent)
284291
end
285292

286293
checksize(p::AdjointFTPlan, x) = checksize(p.parent, x)
287294

288-
unsafe_convert(::Type{Ptr{ft_plan_struct}}, p::AdjointFTPlan{T, FTPlan{T, N, K}}) where {T, N, K} = unsafe_convert(Ptr{ft_plan_struct}, p.parent)
289-
unsafe_convert(::Type{Ptr{mpfr_t}}, p::AdjointFTPlan{T, FTPlan{T, N, K}}) where {T, N, K} = unsafe_convert(Ptr{mpfr_t}, p.parent)
295+
unsafe_convert(::Type{Ptr{ft_plan_struct}}, p::AdjointFTPlan) = unsafe_convert(Ptr{ft_plan_struct}, p.parent)
296+
unsafe_convert(::Type{Ptr{mpfr_t}}, p::AdjointFTPlan) = unsafe_convert(Ptr{mpfr_t}, p.parent)
290297

291-
struct TransposeFTPlan{T, S}
298+
struct TransposeFTPlan{T, S, R}
292299
parent::S
300+
transpose::R
301+
function TransposeFTPlan{T, S, R}(parent::S) where {T, S, R}
302+
new(parent)
303+
end
304+
function TransposeFTPlan{T, S, R}(parent::S, transpose::R) where {T, S, R}
305+
new(parent, transpose)
306+
end
293307
end
294308

295-
TransposeFTPlan(p::FTPlan) = TransposeFTPlan{eltype(p), typeof(p)}(p)
309+
TransposeFTPlan(p::FTPlan) = TransposeFTPlan{eltype(p), typeof(p), typeof(p)}(p)
310+
TransposeFTPlan(p::FTPlan, q::FTPlan) = TransposeFTPlan{eltype(q), typeof(p), typeof(q)}(p, q)
296311

297312
transpose(p::FTPlan) = TransposeFTPlan(p)
298313
transpose(p::TransposeFTPlan) = p.parent
299314

300-
eltype(p::TransposeFTPlan{T, S}) where {T, S} = T
301-
ndims(p::TransposeFTPlan{T, S}) where {T, S} = ndims(p.parent)
302-
function show(io::IO, p::TransposeFTPlan{T, S}) where {T, S}
315+
eltype(p::TransposeFTPlan{T}) where T = T
316+
ndims(p::TransposeFTPlan) = ndims(p.parent)
317+
function show(io::IO, p::TransposeFTPlan)
303318
print(io, "Transpose ")
304319
show(io, p.parent)
305320
end
306321

307322
checksize(p::TransposeFTPlan, x) = checksize(p.parent, x)
308323

309-
unsafe_convert(::Type{Ptr{ft_plan_struct}}, p::TransposeFTPlan{T, FTPlan{T, N, K}}) where {T, N, K} = unsafe_convert(Ptr{ft_plan_struct}, p.parent)
310-
unsafe_convert(::Type{Ptr{mpfr_t}}, p::TransposeFTPlan{T, FTPlan{T, N, K}}) where {T, N, K} = unsafe_convert(Ptr{mpfr_t}, p.parent)
324+
unsafe_convert(::Type{Ptr{ft_plan_struct}}, p::TransposeFTPlan) = unsafe_convert(Ptr{ft_plan_struct}, p.parent)
325+
unsafe_convert(::Type{Ptr{mpfr_t}}, p::TransposeFTPlan) = unsafe_convert(Ptr{mpfr_t}, p.parent)
311326

312327
for f in (:leg2cheb, :cheb2leg, :ultra2ultra, :jac2jac,
313328
:lag2lag, :jac2ultra, :ultra2jac, :jac2cheb,
@@ -547,11 +562,44 @@ function plan_spinsph2fourier(::Type{Complex{Float64}}, n::Integer, s::Integer)
547562
return FTPlan{Complex{Float64}, 2, SPINSPHERE}(plan, n)
548563
end
549564

550-
for (fJ, fC, fE, K) in ((:plan_sph_synthesis, :ft_plan_sph_synthesis, :ft_execute_sph_synthesis, SPHERESYNTHESIS),
551-
(:plan_sph_analysis, :ft_plan_sph_analysis, :ft_execute_sph_analysis, SPHEREANALYSIS),
552-
(:plan_sphv_synthesis, :ft_plan_sphv_synthesis, :ft_execute_sphv_synthesis, SPHEREVSYNTHESIS),
553-
(:plan_sphv_analysis, :ft_plan_sphv_analysis, :ft_execute_sphv_analysis, SPHEREVANALYSIS),
554-
(:plan_disk_synthesis, :ft_plan_disk_synthesis, :ft_execute_disk_synthesis, DISKSYNTHESIS),
565+
for (fJ, fadJ, fC, fE, K) in ((:plan_sph_synthesis, :plan_sph_analysis, :ft_plan_sph_synthesis, :ft_execute_sph_synthesis, SPHERESYNTHESIS),
566+
(:plan_sph_analysis, :plan_sph_synthesis, :ft_plan_sph_analysis, :ft_execute_sph_analysis, SPHEREANALYSIS),
567+
(:plan_sphv_synthesis, :plan_sphv_analysis, :ft_plan_sphv_synthesis, :ft_execute_sphv_synthesis, SPHEREVSYNTHESIS),
568+
(:plan_sphv_analysis, :plan_sphv_synthesis, :ft_plan_sphv_analysis, :ft_execute_sphv_analysis, SPHEREVANALYSIS))
569+
@eval begin
570+
$fJ(x::Matrix{T}) where T = $fJ(T, size(x, 1), size(x, 2))
571+
$fJ(::Type{Complex{T}}, x...) where T <: Real = $fJ(T, x...)
572+
function $fJ(::Type{Float64}, n::Integer, m::Integer)
573+
plan = ccall(($(string(fC)), libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint), n, m)
574+
return FTPlan{Float64, 2, $K}(plan, n, m)
575+
end
576+
adjoint(p::FTPlan{T, 2, $K}) where T = AdjointFTPlan(p, $fadJ(T, p.n, p.m))
577+
transpose(p::FTPlan{T, 2, $K}) where T = TransposeFTPlan(p, $fadJ(T, p.n, p.m))
578+
function lmul!(p::FTPlan{Float64, 2, $K}, x::Matrix{Float64})
579+
if p.n != size(x, 1) || p.m != size(x, 2)
580+
throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.m), x has dimensions $(size(x, 1)) × $(size(x, 2))"))
581+
end
582+
ccall(($(string(fE)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), 'N', p, x, size(x, 1), size(x, 2))
583+
return x
584+
end
585+
function lmul!(p::AdjointFTPlan{Float64, FTPlan{Float64, 2, $K}}, x::Matrix{Float64})
586+
if p.adjoint.n != size(x, 1) || p.adjoint.m != size(x, 2)
587+
throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.m), x has dimensions $(size(x, 1)) × $(size(x, 2))"))
588+
end
589+
ccall(($(string(fE)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), 'T', p.adjoint, x, size(x, 1), size(x, 2))
590+
return x
591+
end
592+
function lmul!(p::TransposeFTPlan{Float64, FTPlan{Float64, 2, $K}}, x::Matrix{Float64})
593+
if p.transpose.n != size(x, 1) || p.transpose.m != size(x, 2)
594+
throw(DimensionMismatch("FTPlan has dimensions $(p.n) × $(p.m), x has dimensions $(size(x, 1)) × $(size(x, 2))"))
595+
end
596+
ccall(($(string(fE)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), 'T', p.transpose, x, size(x, 1), size(x, 2))
597+
return x
598+
end
599+
end
600+
end
601+
602+
for (fJ, fC, fE, K) in ((:plan_disk_synthesis, :ft_plan_disk_synthesis, :ft_execute_disk_synthesis, DISKSYNTHESIS),
555603
(:plan_disk_analysis, :ft_plan_disk_analysis, :ft_execute_disk_analysis, DISKANALYSIS),
556604
(:plan_rectdisk_synthesis, :ft_plan_rectdisk_synthesis, :ft_execute_rectdisk_synthesis, RECTDISKSYNTHESIS),
557605
(:plan_rectdisk_analysis, :ft_plan_rectdisk_analysis, :ft_execute_rectdisk_analysis, RECTDISKANALYSIS),
@@ -788,8 +836,27 @@ end
788836
for (fJ, fC, K) in ((:lmul!, :ft_execute_sph2fourier, SPHERE),
789837
(:ldiv!, :ft_execute_fourier2sph, SPHERE),
790838
(:lmul!, :ft_execute_sphv2fourier, SPHEREV),
791-
(:ldiv!, :ft_execute_fourier2sphv, SPHEREV),
792-
(:lmul!, :ft_execute_disk2cxf, DISK),
839+
(:ldiv!, :ft_execute_fourier2sphv, SPHEREV))
840+
@eval begin
841+
function $fJ(p::FTPlan{Float64, 2, $K}, x::Matrix{Float64})
842+
checksize(p, x)
843+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), 'N', p, x, size(x, 1), size(x, 2))
844+
return x
845+
end
846+
function $fJ(p::AdjointFTPlan{Float64, FTPlan{Float64, 2, $K}}, x::Matrix{Float64})
847+
checksize(p, x)
848+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), 'T', p, x, size(x, 1), size(x, 2))
849+
return x
850+
end
851+
function $fJ(p::TransposeFTPlan{Float64, FTPlan{Float64, 2, $K}}, x::Matrix{Float64})
852+
checksize(p, x)
853+
ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Ptr{ft_plan_struct}, Ptr{Float64}, Cint, Cint), 'T', p, x, size(x, 1), size(x, 2))
854+
return x
855+
end
856+
end
857+
end
858+
859+
for (fJ, fC, K) in ((:lmul!, :ft_execute_disk2cxf, DISK),
793860
(:ldiv!, :ft_execute_cxf2disk, DISK),
794861
(:lmul!, :ft_execute_rectdisk2cheb, RECTDISK),
795862
(:ldiv!, :ft_execute_cheb2rectdisk, RECTDISK),
@@ -896,11 +963,11 @@ for fJ in (:lmul!, :ldiv!)
896963
x .= complex.($fJ(p, real(x)), $fJ(p, imag(x)))
897964
return x
898965
end
899-
function $fJ(p::AdjointFTPlan{T, FTPlan{T, N, K}}, x::AbstractArray{Complex{T}}) where {T, N, K}
966+
function $fJ(p::AdjointFTPlan{T}, x::AbstractArray{Complex{T}}) where T
900967
x .= complex.($fJ(p, real(x)), $fJ(p, imag(x)))
901968
return x
902969
end
903-
function $fJ(p::TransposeFTPlan{T, FTPlan{T, N, K}}, x::AbstractArray{Complex{T}}) where {T, N, K}
970+
function $fJ(p::TransposeFTPlan{T}, x::AbstractArray{Complex{T}}) where T
904971
x .= complex.($fJ(p, real(x)), $fJ(p, imag(x)))
905972
return x
906973
end

0 commit comments

Comments
 (0)