cleanup C - Julia BigFloat interop (#104)

* cleanup C - Julia BigFloat interop

* rework help text

Author: MikaelSlevinsky

Project.toml

@@ -1,6 +1,6 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.9.1"
+version = "0.9.2"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
@@ -22,7 +22,7 @@ BinaryProvider = "0.5"
 DSP = "0.6"
 FFTW = "1"
 FastGaussQuadrature = "0.4"
-FastTransforms_jll = "0.3.1"
+FastTransforms_jll = "0.3.2"
 Reexport = "0.2"
 SpecialFunctions = "0.8, 0.9, 0.10"
 ToeplitzMatrices = "0.6"

deps/build.jl

@@ -1,7 +1,7 @@
 using BinaryProvider
 import Libdl
 
-version = v"0.3.1"
+version = v"0.3.2"
 
 if arch(platform_key_abi()) != :x86_64
     @warn "FastTransforms has only been tested on x86_64 architectures."

src/FastTransforms.jl

@@ -12,7 +12,6 @@ import Base: unsafe_convert, eltype, ndims, adjoint, transpose, show, *, \,
              inv, size, view
 
 import Base.GMP: Limb
-import Base.MPFR: BigFloat, _BigFloat
 
 import AbstractFFTs: Plan, ScaledPlan,
                      fft, ifft, bfft, fft!, ifft!, bfft!,

src/libfasttransforms.jl

@@ -32,13 +32,21 @@ struct mpfr_t <: AbstractFloat
     d::Ptr{Limb}
 end
 
-mpfr_t(x::BigFloat) = mpfr_t(x.prec, x.sign, x.exp, x.d)
+"""
+`BigFloat` is a mutable struct and there is no guarantee that each entry in
+an `Array{BigFloat}` has unique pointers. For example, looking at the `Limb`s,
+
+    Id = Matrix{BigFloat}(I, 3, 3)
+    map(x->x.d, Id)
 
-function BigFloat(x::mpfr_t)
-    nb = ccall((:mpfr_custom_get_size,:libmpfr), Csize_t, (Clong,), precision(BigFloat))
-    nb = (nb + Core.sizeof(Limb) - 1) ÷ Core.sizeof(Limb) # align to number of Limb allocations required for this
-    str = unsafe_string(Ptr{UInt8}(x.d), nb * Core.sizeof(Limb))
-    _BigFloat(x.prec, x.sign, x.exp, str)
+shows that the ones and the zeros all share the same pointers. If a C function
+assumes unicity of each datum, then the array must be renewed with a `deepcopy`.
+"""
+function renew!(x::Array{BigFloat})
+    for i in eachindex(x)
+        @inbounds x[i] = deepcopy(x[i])
+    end
+    return x
 end
 
 set_num_threads(n::Integer) = ccall((:ft_set_num_threads, libfasttransforms), Cvoid, (Cint, ), n)
@@ -602,31 +610,22 @@ for (fJ, fC, elty) in ((:lmul!, :ft_bfmvf, :Float32),
     end
 end
 
-for (fJ, fC) in ((:lmul!, :ft_mpfr_trmv),
-                 (:ldiv!, :ft_mpfr_trsv))
+for (fJ, fC) in ((:lmul!, :ft_mpfr_trmv_ptr),
+                 (:ldiv!, :ft_mpfr_trsv_ptr))
     @eval begin
         function $fJ(p::FTPlan{BigFloat, 1}, x::Vector{BigFloat})
             checksize(p, x)
-            xt = deepcopy.(x)
-            xc = mpfr_t.(xt)
-            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Int32), 'N', p.n, p, p.n, xc, Base.MPFR.ROUNDING_MODE[])
-            x .= BigFloat.(xc)
+            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{BigFloat}, Int32), 'N', p.n, p, p.n, renew!(x), Base.MPFR.ROUNDING_MODE[])
             return x
         end
         function $fJ(p::AdjointFTPlan{BigFloat, FTPlan{BigFloat, 1, K}}, x::Vector{BigFloat}) where K
             checksize(p, x)
-            xt = deepcopy.(x)
-            xc = mpfr_t.(xt)
-            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Int32), 'T', p.parent.n, p, p.parent.n, xc, Base.MPFR.ROUNDING_MODE[])
-            x .= BigFloat.(xc)
+            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{BigFloat}, Int32), 'T', p.parent.n, p, p.parent.n, renew!(x), Base.MPFR.ROUNDING_MODE[])
             return x
         end
         function $fJ(p::TransposeFTPlan{BigFloat, FTPlan{BigFloat, 1, K}}, x::Vector{BigFloat}) where K
             checksize(p, x)
-            xt = deepcopy.(x)
-            xc = mpfr_t.(xt)
-            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Int32), 'T', p.parent.n, p, p.parent.n, xc, Base.MPFR.ROUNDING_MODE[])
-            x .= BigFloat.(xc)
+            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{BigFloat}, Int32), 'T', p.parent.n, p, p.parent.n, renew!(x), Base.MPFR.ROUNDING_MODE[])
             return x
         end
     end
@@ -655,31 +654,22 @@ for (fJ, fC, elty) in ((:lmul!, :ft_bfmmf, :Float32),
     end
 end
 
-for (fJ, fC) in ((:lmul!, :ft_mpfr_trmm),
-                 (:ldiv!, :ft_mpfr_trsm))
+for (fJ, fC) in ((:lmul!, :ft_mpfr_trmm_ptr),
+                 (:ldiv!, :ft_mpfr_trsm_ptr))
     @eval begin
         function $fJ(p::FTPlan{BigFloat, 1}, x::Matrix{BigFloat})
             checksize(p, x)
-            xt = deepcopy.(x)
-            xc = mpfr_t.(xt)
-            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Cint, Cint, Int32), 'N', p.n, p, p.n, xc, size(x, 1), size(x, 2), Base.MPFR.ROUNDING_MODE[])
-            x .= BigFloat.(xc)
+            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{BigFloat}, Cint, Cint, Int32), 'N', p.n, p, p.n, renew!(x), size(x, 1), size(x, 2), Base.MPFR.ROUNDING_MODE[])
             return x
         end
         function $fJ(p::AdjointFTPlan{BigFloat, FTPlan{BigFloat, 1, K}}, x::Matrix{BigFloat}) where K
             checksize(p, x)
-            xt = deepcopy.(x)
-            xc = mpfr_t.(xt)
-            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Cint, Cint, Int32), 'T', p.parent.n, p, p.parent.n, xc, size(x, 1), size(x, 2), Base.MPFR.ROUNDING_MODE[])
-            x .= BigFloat.(xc)
+            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{BigFloat}, Cint, Cint, Int32), 'T', p.parent.n, p, p.parent.n, renew!(x), size(x, 1), size(x, 2), Base.MPFR.ROUNDING_MODE[])
             return x
         end
         function $fJ(p::TransposeFTPlan{BigFloat, FTPlan{BigFloat, 1, K}}, x::Matrix{BigFloat}) where K
             checksize(p, x)
-            xt = deepcopy.(x)
-            xc = mpfr_t.(xt)
-            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{mpfr_t}, Cint, Cint, Int32), 'T', p.parent.n, p, p.parent.n, xc, size(x, 1), size(x, 2), Base.MPFR.ROUNDING_MODE[])
-            x .= BigFloat.(xc)
+            ccall(($(string(fC)), libfasttransforms), Cvoid, (Cint, Cint, Ptr{mpfr_t}, Cint, Ptr{BigFloat}, Cint, Cint, Int32), 'T', p.parent.n, p, p.parent.n, renew!(x), size(x, 1), size(x, 2), Base.MPFR.ROUNDING_MODE[])
             return x
         end
     end

test/libfasttransformstests.jl

@@ -24,7 +24,7 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
     end
 
     α, β, γ, δ, λ, μ = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6
-    function test_1d_plans(p1, p2, x; skip::Bool=false)
+    function test_1d_plans(p1, p2, x)
         y = p1*x
         z = p2*y
         @test z ≈ x
@@ -44,51 +44,37 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
         @test z ≈ x
         P = p1*I
         Q = p2*P
-        skip ? (@test_skip Q ≈ I) : (@test Q ≈ I)
+        @test Q ≈ I
         P = p1*I
         Q = p1'P
         P = transpose(p1)*Q
         Q = transpose(p1)\P
         P = p1'\Q
         Q = p1\P
-        skip ? (@test_skip Q ≈ I) : (@test Q ≈ I)
+        @test Q ≈ I
         P = p2*I
         Q = p2'P
         P = transpose(p2)*Q
         Q = transpose(p2)\P
         P = p2'\Q
         Q = p2\P
-        skip ? (@test_skip Q ≈ I) : (@test Q ≈ I)
+        @test Q ≈ I
     end
 
-    for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
+    for T in (Float32, Float64, Complex{Float32}, Complex{Float64}, BigFloat, Complex{BigFloat})
         x = T(1)./(1:n)
         Id = Matrix{T}(I, n, n)
         for (p1, p2) in ((plan_leg2cheb(Id), plan_cheb2leg(Id)),
-                          (plan_ultra2ultra(Id, λ, μ), plan_ultra2ultra(Id, μ, λ)),
-                          (plan_jac2jac(Id, α, β, γ, δ), plan_jac2jac(Id, γ, δ, α, β)),
-                          (plan_lag2lag(Id, α, β), plan_lag2lag(Id, β, α)),
-                          (plan_jac2ultra(Id, α, β, λ), plan_ultra2jac(Id, λ, α, β)),
-                          (plan_jac2cheb(Id, α, β), plan_cheb2jac(Id, α, β)),
-                          (plan_ultra2cheb(Id, λ), plan_cheb2ultra(Id, λ)))
+                         (plan_ultra2ultra(Id, λ, μ), plan_ultra2ultra(Id, μ, λ)),
+                         (plan_jac2jac(Id, α, β, γ, δ), plan_jac2jac(Id, γ, δ, α, β)),
+                         (plan_lag2lag(Id, α, β), plan_lag2lag(Id, β, α)),
+                         (plan_jac2ultra(Id, α, β, λ), plan_ultra2jac(Id, λ, α, β)),
+                         (plan_jac2cheb(Id, α, β), plan_cheb2jac(Id, α, β)),
+                         (plan_ultra2cheb(Id, λ), plan_cheb2ultra(Id, λ)))
             test_1d_plans(p1, p2, x)
         end
     end
 
-    for T in (BigFloat, Complex{BigFloat})
-        x = T(1)./(1:n)
-        Id = Matrix{T}(I, n, n)
-        for (p1, p2) in ((plan_leg2cheb(Id), plan_cheb2leg(Id)),
-                          (plan_ultra2ultra(Id, λ, μ), plan_ultra2ultra(Id, μ, λ)),
-                          (plan_jac2jac(Id, α, β, γ, δ), plan_jac2jac(Id, γ, δ, α, β)),
-                          (plan_lag2lag(Id, α, β), plan_lag2lag(Id, β, α)),
-                          (plan_jac2ultra(Id, α, β, λ), plan_ultra2jac(Id, λ, α, β)),
-                          (plan_jac2cheb(Id, α, β), plan_cheb2jac(Id, α, β)),
-                          (plan_ultra2cheb(Id, λ), plan_cheb2ultra(Id, λ)))
-            test_1d_plans(p1, p2, x; skip=true)
-        end
-    end
-
     function test_nd_plans(p, ps, pa, A)
         B = copy(A)
         C = ps*(p*A)