add complex transforms, part deux

Closes #79

Author: MikaelSlevinsky

src/libfasttransforms.jl

@@ -436,6 +436,7 @@ for (fJ, fC, fE, K) in ((:plan_sph_synthesis, :ft_plan_sph_synthesis, :ft_execut
                     (:plan_tri_analysis, :ft_plan_tri_analysis, :ft_execute_tri_analysis, TRIANGLEANALYSIS))
     @eval begin
         $fJ(x::Matrix{T}) where T = $fJ(T, size(x, 1), size(x, 2))
+        $fJ(x::Matrix{Complex{T}}) where T <: Real = $fJ(T, size(x, 1), size(x, 2))
         function $fJ(::Type{Float64}, n::Integer, m::Integer)
             plan = ccall(($(string(fC)), libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint), n, m)
             return FTPlan{Float64, 2, $K}(plan, n, m)
@@ -451,6 +452,7 @@ for (fJ, fC, fE, K) in ((:plan_sph_synthesis, :ft_plan_sph_synthesis, :ft_execut
 end
 
 plan_tet_synthesis(x::Array{T, 3}) where T = plan_tet_synthesis(T, size(x, 1), size(x, 2), size(x, 3))
+plan_tet_synthesis(x::Array{Complex{T}, 3}) where T <: Real = plan_tet_synthesis(T, size(x, 1), size(x, 2), size(x, 3))
 
 function plan_tet_synthesis(::Type{Float64}, n::Integer, l::Integer, m::Integer)
     plan = ccall((:ft_plan_tet_synthesis, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), n, l, m)
@@ -466,6 +468,7 @@ function lmul!(p::FTPlan{Float64, 3, TETRAHEDRONSYNTHESIS}, x::Array{Float64, 3}
 end
 
 plan_tet_analysis(x::Array{T, 3}) where T = plan_tet_analysis(T, size(x, 1), size(x, 2), size(x, 3))
+plan_tet_analysis(x::Array{Complex{T}, 3}) where T <: Real = plan_tet_analysis(T, size(x, 1), size(x, 2), size(x, 3))
 
 function plan_tet_analysis(::Type{Float64}, n::Integer, l::Integer, m::Integer)
     plan = ccall((:ft_plan_tet_analysis, libfasttransforms), Ptr{ft_plan_struct}, (Cint, Cint, Cint), n, l, m)

test/libfasttransformstests.jl

@@ -5,7 +5,7 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
 @testset "libfasttransforms" begin
     n = 64
     α, β, γ, δ, λ, μ = 0.1, 0.2, 0.3, 0.4, 0.5, 0.6
-    for T in (Float32, Float64, BigFloat, Complex{Float32}, Complex{Float64}), Complex{BigFloat})
+    for T in (Float32, Float64, Complex{Float32}, Complex{Float64})
         x = T(1)./(1:n)
         Id = Matrix{T}(I, n, n)
         for (p1, p2) in ((plan_leg2cheb(Id), plan_cheb2leg(Id)),
@@ -32,7 +32,6 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
             y = p2'\z
             z = p2\y
             @test z ≈ x
-            Id = Matrix{T}(I, n, n)
             P = p1*Id
             Q = p2*P
             @test Q ≈ Id
@@ -53,8 +52,55 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
         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, λ)))
+            y = p1*x
+            z = p2*y
+            @test z ≈ x
+            y = p1*x
+            z = p1'y
+            y = transpose(p1)*z
+            z = transpose(p1)\y
+            y = p1'\z
+            z = p1\y
+            @test z ≈ x
+            y = p2*x
+            z = p2'y
+            y = transpose(p2)*z
+            z = transpose(p2)\y
+            y = p2'\z
+            z = p2\y
+            @test z ≈ x
+            P = p1*Id
+            Q = p2*P
+            @test_skip Q ≈ Id
+            P = p1*Id
+            Q = p1'P
+            P = transpose(p1)*Q
+            Q = transpose(p1)\P
+            P = p1'\Q
+            Q = p1\P
+            @test_skip Q ≈ Id
+            P = p2*Id
+            Q = p2'P
+            P = transpose(p2)*Q
+            Q = transpose(p2)\P
+            P = p2'\Q
+            Q = p2\P
+            @test_skip Q ≈ Id
+        end
+    end
+
     for T in (Float64, Complex{Float64})
-        A = sphones(T, n, 2n-1)
+        A = T <: Real ? sphones(T, n, 2n-1) : sphones(T, n, 2n-1) + im*sphones(T, n, 2n-1)
         p = plan_sph2fourier(A)
         ps = plan_sph_synthesis(A)
         pa = plan_sph_analysis(A)
@@ -63,7 +109,7 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
         A = p\(pa*C)
         @test A ≈ B
 
-        A = sphvones(T, n, 2n-1)
+        A = T <: Real ? sphvones(T, n, 2n-1) : sphvones(T, n, 2n-1) + im*sphvones(T, n, 2n-1)
         p = plan_sphv2fourier(A)
         ps = plan_sphv_synthesis(A)
         pa = plan_sphv_analysis(A)
@@ -72,7 +118,7 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
         A = p\(pa*C)
         @test A ≈ B
 
-        A = diskones(T, n, 4n-3)
+        A = T <: Real ? diskones(T, n, 4n-3) : diskones(T, n, 4n-3) + im*diskones(T, n, 4n-3)
         p = plan_disk2cxf(A)
         ps = plan_disk_synthesis(A)
         pa = plan_disk_analysis(A)
@@ -81,7 +127,7 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
         A = p\(pa*C)
         @test A ≈ B
 
-        A = triones(T, n, n)
+        A = T <: Real ? triones(T, n, n) : triones(T, n, n) + im*triones(T, n, n)
         p = plan_tri2cheb(A, α, β, γ)
         ps = plan_tri_synthesis(A)
         pa = plan_tri_analysis(A)
@@ -90,7 +136,7 @@ FastTransforms.set_num_threads(ceil(Int, Base.Sys.CPU_THREADS/2))
         A = p\(pa*C)
         @test A ≈ B
 
-        A = tetones(T, n, n, n)
+        A = T <: Real ? tetones(T, n, n, n) : tetones(T, n, n, n) + im*tetones(T, n, n, n)
         p = plan_tet2cheb(A, α, β, γ, δ)
         ps = plan_tet_synthesis(A)
         pa = plan_tet_analysis(A)