Merge pull request #35 from MikaelSlevinsky/feat-triangular-harmonics

Start triangular harmonics, support experimental native multithreading

Author: MikaelSlevinsky

.gitignore

@@ -1,2 +1,2 @@
 docs/build/
-docs/site/
+docs/site/

LICENSE.md

@@ -1,6 +1,8 @@
 The FastTransforms.jl package is licensed under the MIT "Expat" License:
 
-> Copyright (c) 2017: Richard Mikael Slevinsky.
+> Copyright (c) 2016-2018: Richard Mikael Slevinsky and other contributors:
+>
+> https://github.com/MikaelSlevinsky/FastTransforms.jl/graphs/contributors
 >
 > Permission is hereby granted, free of charge, to any person obtaining
 > a copy of this software and associated documentation files (the

README.md

@@ -166,14 +166,16 @@ As with other fast transforms, `plan_sph2fourier` saves effort by caching the pr
 
    [1]  B. Alpert and V. Rokhlin. <a href="http://dx.doi.org/10.1137/0912009">A fast algorithm for the evaluation of Legendre expansions</a>, *SIAM J. Sci. Stat. Comput.*, **12**:158—179, 1991.
 
-   [2]  N. Hale and A. Townsend. <a href="http://dx.doi.org/10.1137/130932223">A fast, simple, and stable Chebyshev—Legendre transform using and asymptotic formula</a>, *SIAM J. Sci. Comput.*, **36**:A148—A167, 2014.
+   [2]  N. Hale and A. Townsend. <a href="http://dx.doi.org/10.1137/130932223">A fast, simple, and stable Chebyshev—Legendre transform using an asymptotic formula</a>, *SIAM J. Sci. Comput.*, **36**:A148—A167, 2014.
 
    [3]  J. Keiner. <a href="http://dx.doi.org/10.1137/070703065">Computing with expansions in Gegenbauer polynomials</a>, *SIAM J. Sci. Comput.*, **31**:2151—2171, 2009.
 
    [4]  D. Ruiz—Antolín and A. Townsend. <a href="https://arxiv.org/abs/1701.04492">A nonuniform fast Fourier transform based on low rank approximation</a>, arXiv:1701.04492, 2017.
 
-   [5]  R. M. Slevinsky. <a href="https://doi.org/10.1093/imanum/drw070">On the use of Hahn's asymptotic formula and stabilized recurrence for a fast, simple, and stable Chebyshev—Jacobi transform</a>, in press at *IMA J. Numer. Anal.*, 2017.
+   [5]  R. M. Slevinsky. <a href="https://doi.org/10.1093/imanum/drw070">On the use of Hahn's asymptotic formula and stabilized recurrence for a fast, simple, and stable Chebyshev—Jacobi transform</a>, *IMA J. Numer. Anal.*, **38**:102—124, 2018.
 
-   [6]  R. M. Slevinsky. <a href="https://arxiv.org/abs/1705.05448">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, arXiv:1705.05448, 2017.
+   [6]  R. M. Slevinsky. <a href="https://doi.org/10.1016/j.acha.2017.11.001">Fast and backward stable transforms between spherical harmonic expansions and bivariate Fourier series</a>, in press at *Appl. Comput. Harmon. Anal.*, 2017.
 
-   [7]  A. Townsend, M. Webb, and S. Olver. <a href="https://doi.org/10.1090/mcom/3277">Fast polynomial transforms based on Toeplitz and Hankel matrices</a>, in press at *Math. Comp.*, 2017.
+   [7]  R. M. Slevinsky, <a href="https://arxiv.org/abs/1711.07866">Conquering the pre-computation in two-dimensional harmonic polynomial transforms</a>, arXiv:1711.07866, 2017.
+
+   [8]  A. Townsend, M. Webb, and S. Olver. <a href="https://doi.org/10.1090/mcom/3277">Fast polynomial transforms based on Toeplitz and Hankel matrices</a>, in press at *Math. Comp.*, 2017.

REQUIRE

@@ -1,6 +1,6 @@
 julia 0.6
 ToeplitzMatrices 0.2
-HierarchicalMatrices 0.0.2
+HierarchicalMatrices 0.1
 LowRankApprox 0.0.2
 ProgressMeter 0.3.4
 SpecialFunctions 0.3.4

src/FastTransforms.jl

@@ -42,6 +42,10 @@ export SlowSphericalHarmonicPlan, FastSphericalHarmonicPlan, ThinSphericalHarmon
 export sph2fourier, fourier2sph, plan_sph2fourier
 export sphones, sphzeros, sphrand, sphrandn, sphevaluate
 
+export SlowTriangularHarmonicPlan
+export tri2cheb, cheb2tri, plan_tri2cheb
+export triones, trizeros, trirand, trirandn, trievaluate
+
 # Other module methods and constants:
 #export ChebyshevJacobiPlan, jac2cheb, cheb2jac
 #export sqrtpi, pochhammer, stirlingseries, stirlingremainder, Aratio, Cratio, Anαβ
@@ -51,6 +55,7 @@ export sphones, sphzeros, sphrand, sphrandn, sphevaluate
 #export fejer2, fejer_plan2, fejerweights2
 #export RecurrencePlan, forward_recurrence!, backward_recurrence
 
+include("stepthreading.jl")
 include("fftBigFloat.jl")
 include("specialfunctions.jl")
 include("clenshawcurtis.jl")
@@ -82,6 +87,7 @@ jac2jac(x...)=th_jac2jac(x...)
 
 include("hierarchical.jl")
 include("SphericalHarmonics/SphericalHarmonics.jl")
+include("TriangularHarmonics/TriangularHarmonics.jl")
 
 include("gaunt.jl")
 

src/SphericalHarmonics/fastplan.jl

@@ -8,10 +8,10 @@ struct FastSphericalHarmonicPlan{T} <: SphericalHarmonicPlan{T}
     B::Matrix{T}
 end
 
-function FastSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int; opts...)
+function FastSphericalHarmonicPlan(A::Matrix{T}, L::Int; opts...) where T
     M, N = size(A)
     n = (N+1)÷2
-    RP = RotationPlan(T, n-1)
+    RP = RotationPlan(T, M-1)
     p1 = plan_normleg2cheb(A)
     p2 = plan_normleg12cheb2(A)
     p1inv = plan_cheb2normleg(A)
@@ -77,7 +77,7 @@ function Base.At_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix)
         At_mul_B_col_J!(Y, BF[J-1], B, 2J)
         At_mul_B_col_J!(Y, BF[J-1], B, 2J+1)
     end
-    zero_spurious_modes!(Y)
+    sph_zero_spurious_modes!(Y)
 end
 
 Base.Ac_mul_B!(Y::Matrix, FP::FastSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, FP, X)

src/SphericalHarmonics/slowplan.jl

@@ -2,7 +2,7 @@ import Base.LinAlg: Givens, AbstractRotation
 
 ### These three A_mul_B! should be in Base, but for the time being they do not add methods to Base.A_mul_B!; they add methods to the internal A_mul_B!.
 
-function A_mul_B!{T<:Real}(G::Givens{T}, A::AbstractVecOrMat)
+function A_mul_B!(G::Givens{T}, A::AbstractVecOrMat) where T<:Real
     m, n = size(A, 1), size(A, 2)
     if G.i2 > m
         throw(DimensionMismatch("column indices for rotation are outside the matrix"))
@@ -28,7 +28,7 @@ function A_mul_B!(A::AbstractMatrix, G::Givens)
     return A
 end
 
-function A_mul_B!{T<:Real}(A::AbstractMatrix, G::Givens{T})
+function A_mul_B!(A::AbstractMatrix, G::Givens{T}) where T<:Real
     m, n = size(A, 1), size(A, 2)
     if G.i2 > n
         throw(DimensionMismatch("column indices for rotation are outside the matrix"))
@@ -45,7 +45,7 @@ struct Pnmp2toPlm{T} <: AbstractRotation{T}
     rotations::Vector{Givens{T}}
 end
 
-function Pnmp2toPlm{T}(::Type{T}, n::Int, m::Int)
+function Pnmp2toPlm(::Type{T}, n::Int, m::Int) where T
     G = Vector{Givens{T}}(n)
     @inbounds for ℓ = 1:n
         c = sqrt(T((2m+2)*(2ℓ+2m+3))/T((ℓ+2m+2)*(ℓ+2m+3)))
@@ -75,16 +75,77 @@ end
 
 struct RotationPlan{T} <: AbstractRotation{T}
     layers::Vector{Pnmp2toPlm{T}}
+    snm::Vector{T}
+    cnm::Vector{T}
 end
 
-function RotationPlan{T}(::Type{T}, n::Int)
+function RotationPlan(::Type{T}, n::Int) where T
     layers = Vector{Pnmp2toPlm{T}}(n-1)
     @inbounds for m = 0:n-2
         layers[m+1] = Pnmp2toPlm(T, n-1-m, m)
     end
-    RotationPlan(layers)
+    n = n+1
+    snm = zeros(T, (n*(n+1))÷2)
+    cnm = zeros(T, (n*(n+1))÷2)
+    @inbounds for l = 0:n-1
+        for m = 0:n-l-1
+            nums = T((l+1)*(l+2))
+            numc = T((2*m+2)*(2*l+2*m+5))
+            den = T((l+2*m+3)*(l+2*m+4))
+            snm[l+(m*(2*n+1-m))÷2+1] = sqrt(nums/den)
+            cnm[l+(m*(2*n+1-m))÷2+1] = sqrt(numc/den)
+        end
+    end
+    RotationPlan(layers, snm, cnm)
+end
+
+function Base.A_mul_B!(P::RotationPlan, A::AbstractMatrix)
+    N, M = size(A)
+    snm = P.snm
+    cnm = P.cnm
+    @stepthreads for m = M÷2:-1:2
+        @inbounds for j = m:-2:2
+            for l = N-j:-1:1
+                s = snm[l+(j-2)*(2*N+3-j)÷2]
+                c = cnm[l+(j-2)*(2*N+3-j)÷2]
+                a1 = A[l+N*(2*m-1)]
+                a2 = A[l+2+N*(2*m-1)]
+                a3 = A[l+N*(2*m)]
+                a4 = A[l+2+N*(2*m)]
+                A[l+N*(2*m-1)] = c*a1 + s*a2
+                A[l+2+N*(2*m-1)] = c*a2 - s*a1
+                A[l+N*(2*m)] = c*a3 + s*a4
+                A[l+2+N*(2*m)] = c*a4 - s*a3
+            end
+        end
+    end
+    A
+end
+
+function Base.At_mul_B!(P::RotationPlan, A::AbstractMatrix)
+    N, M = size(A)
+    snm = P.snm
+    cnm = P.cnm
+    @stepthreads for m = M÷2:-1:2
+        @inbounds for j = reverse(m:-2:2)
+            for l = 1:N-j
+                s = snm[l+(j-2)*(2*N+3-j)÷2]
+                c = cnm[l+(j-2)*(2*N+3-j)÷2]
+                a1 = A[l+N*(2*m-1)]
+                a2 = A[l+2+N*(2*m-1)]
+                a3 = A[l+N*(2*m)]
+                a4 = A[l+2+N*(2*m)]
+                A[l+N*(2*m-1)] = c*a1 - s*a2
+                A[l+2+N*(2*m-1)] = c*a2 + s*a1
+                A[l+N*(2*m)] = c*a3 - s*a4
+                A[l+2+N*(2*m)] = c*a4 + s*a3
+            end
+        end
+    end
+    A
 end
 
+#=
 function Base.A_mul_B!(P::RotationPlan, A::AbstractMatrix)
     M, N = size(A)
     @inbounds for m = N÷2-2:-1:0
@@ -122,6 +183,7 @@ function Base.At_mul_B!(P::RotationPlan, A::AbstractMatrix)
     end
     A
 end
+=#
 
 Base.Ac_mul_B!(P::RotationPlan, A::AbstractMatrix) = At_mul_B!(P, A)
 
@@ -135,10 +197,9 @@ struct SlowSphericalHarmonicPlan{T} <: SphericalHarmonicPlan{T}
     B::Matrix{T}
 end
 
-function SlowSphericalHarmonicPlan{T}(A::Matrix{T})
+function SlowSphericalHarmonicPlan(A::Matrix{T}) where T
     M, N = size(A)
-    n = (N+1)÷2
-    RP = RotationPlan(T, n-1)
+    RP = RotationPlan(T, M-1)
     p1 = plan_normleg2cheb(A)
     p2 = plan_normleg12cheb2(A)
     p1inv = plan_cheb2normleg(A)
@@ -177,7 +238,7 @@ function Base.At_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix)
         A_mul_B_col_J!!(Y, p1inv, B, J)
         A_mul_B_col_J!!(Y, p1inv, B, J+1)
     end
-    zero_spurious_modes!(At_mul_B!(RP, Y))
+    sph_zero_spurious_modes!(At_mul_B!(RP, Y))
 end
 
 Base.Ac_mul_B!(Y::Matrix, SP::SlowSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, SP, X)

src/SphericalHarmonics/sphfunctions.jl

@@ -1,4 +1,4 @@
-function zero_spurious_modes!(A::AbstractMatrix)
+function sph_zero_spurious_modes!(A::AbstractMatrix)
     M, N = size(A)
     n = N÷2
     @inbounds for j = 1:n
@@ -10,7 +10,7 @@ function zero_spurious_modes!(A::AbstractMatrix)
     A
 end
 
-function sphrand{T}(::Type{T}, m::Int, n::Int)
+function sphrand(::Type{T}, m::Int, n::Int) where T
     A = zeros(T, m, 2n-1)
     for i = 1:m
         A[i,1] = rand(T)
@@ -24,7 +24,7 @@ function sphrand{T}(::Type{T}, m::Int, n::Int)
     A
 end
 
-function sphrandn{T}(::Type{T}, m::Int, n::Int)
+function sphrandn(::Type{T}, m::Int, n::Int) where T
     A = zeros(T, m, 2n-1)
     for i = 1:m
         A[i,1] = randn(T)
@@ -38,7 +38,7 @@ function sphrandn{T}(::Type{T}, m::Int, n::Int)
     A
 end
 
-function sphones{T}(::Type{T}, m::Int, n::Int)
+function sphones(::Type{T}, m::Int, n::Int) where T
     A = zeros(T, m, 2n-1)
     for i = 1:m
         A[i,1] = one(T)
@@ -52,7 +52,7 @@ function sphones{T}(::Type{T}, m::Int, n::Int)
     A
 end
 
-sphzeros{T}(::Type{T}, m::Int, n::Int) = zeros(T, m, 2n-1)
+sphzeros(::Type{T}, m::Int, n::Int) where T = zeros(T, m, 2n-1)
 
 function normalizecolumns!(A::AbstractMatrix)
     m, n = size(A)

src/SphericalHarmonics/synthesisanalysis.jl

@@ -5,7 +5,7 @@ struct SynthesisPlan{T, P1, P2}
     temp::Vector{T}
 end
 
-function plan_synthesis{T<:fftwNumber}(A::Matrix{T})
+function plan_synthesis(A::Matrix{T}) where T<:fftwNumber
     m, n = size(A)
     x = FFTW.FakeArray(T, m)
     y = FFTW.FakeArray(T, n)
@@ -22,7 +22,7 @@ struct AnalysisPlan{T, P1, P2}
     temp::Vector{T}
 end
 
-function plan_analysis{T<:fftwNumber}(A::Matrix{T})
+function plan_analysis(A::Matrix{T}) where T<:fftwNumber
     m, n = size(A)
     x = FFTW.FakeArray(T, m)
     y = FFTW.FakeArray(T, n)
@@ -32,7 +32,7 @@ function plan_analysis{T<:fftwNumber}(A::Matrix{T})
     AnalysisPlan(planθ, planφ, C, zeros(T, n))
 end
 
-function Base.A_mul_B!{T}(Y::Matrix{T}, P::SynthesisPlan{T}, X::Matrix{T})
+function Base.A_mul_B!(Y::Matrix{T}, P::SynthesisPlan{T}, X::Matrix{T}) where T
     M, N = size(X)
 
     # Column synthesis
@@ -73,7 +73,7 @@ function Base.A_mul_B!{T}(Y::Matrix{T}, P::SynthesisPlan{T}, X::Matrix{T})
     Y
 end
 
-function Base.A_mul_B!{T}(Y::Matrix{T}, P::AnalysisPlan{T}, X::Matrix{T})
+function Base.A_mul_B!(Y::Matrix{T}, P::AnalysisPlan{T}, X::Matrix{T}) where T
     M, N = size(X)
 
     # Row analysis
@@ -112,7 +112,7 @@ end
 
 
 
-function row_analysis!{T}(P, C, vals::Vector{T})
+function row_analysis!(P, C, vals::Vector{T}) where T
     n = length(vals)
     cfs = scale!(two(T)/n,P*vals)
     cfs[1] *= half(T)
@@ -123,7 +123,7 @@ function row_analysis!{T}(P, C, vals::Vector{T})
     negateeven!(reverseeven!(A_mul_B!(C, cfs)))
 end
 
-function row_synthesis!{T}(P, C, cfs::Vector{T})
+function row_synthesis!(P, C, cfs::Vector{T}) where T
     n = length(cfs)
     Ac_mul_B!(C, reverseeven!(negateeven!(cfs)))
     if iseven(n)
@@ -171,17 +171,17 @@ function negateeven!(x::Vector)
     x
 end
 
-function A_mul_B_col_J!{T}(Y::Matrix{T}, P::r2rFFTWPlan{T}, X::Matrix{T}, J::Int)
+function A_mul_B_col_J!(Y::Matrix{T}, P::r2rFFTWPlan{T}, X::Matrix{T}, J::Int) where T
     unsafe_execute_col_J!(P, X, Y, J)
     return Y
 end
 
-function unsafe_execute_col_J!{T<:fftwDouble}(plan::r2rFFTWPlan{T}, X::Matrix{T}, Y::Matrix{T}, J::Int)
+function unsafe_execute_col_J!(plan::r2rFFTWPlan{T}, X::Matrix{T}, Y::Matrix{T}, J::Int) where T<:fftwDouble
     M = size(X, 1)
     ccall((:fftw_execute_r2r, libfftw), Void, (PlanPtr, Ptr{T}, Ptr{T}), plan, pointer(X, M*(J-1)+1), pointer(Y, M*(J-1)+1))
 end
 
-function unsafe_execute_col_J!{T<:fftwSingle}(plan::r2rFFTWPlan{T}, X::Matrix{T}, Y::Matrix{T}, J::Int)
+function unsafe_execute_col_J!(plan::r2rFFTWPlan{T}, X::Matrix{T}, Y::Matrix{T}, J::Int) where T<:fftwSingle
     M = size(X, 1)
     ccall((:fftwf_execute_r2r, libfftwf), Void, (PlanPtr, Ptr{T}, Ptr{T}), plan, pointer(X, M*(J-1)+1), pointer(Y, M*(J-1)+1))
 end

src/SphericalHarmonics/thinplan.jl

@@ -12,10 +12,10 @@ struct ThinSphericalHarmonicPlan{T} <: SphericalHarmonicPlan{T}
     B::Matrix{T}
 end
 
-function ThinSphericalHarmonicPlan{T}(A::Matrix{T}, L::Int; opts...)
+function ThinSphericalHarmonicPlan(A::Matrix{T}, L::Int; opts...) where T
     M, N = size(A)
     n = (N+1)÷2
-    RP = RotationPlan(T, n-1)
+    RP = RotationPlan(T, M-1)
     p1 = plan_normleg2cheb(A)
     p2 = plan_normleg12cheb2(A)
     p1inv = plan_cheb2normleg(A)
@@ -148,7 +148,7 @@ function Base.At_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix)
         end
     end
 
-    zero_spurious_modes!(Y)
+    sph_zero_spurious_modes!(Y)
 end
 
 Base.Ac_mul_B!(Y::Matrix, TP::ThinSphericalHarmonicPlan, X::Matrix) = At_mul_B!(Y, TP, X)

src/TriangularHarmonics/TriangularHarmonics.jl

@@ -0,0 +1,24 @@
+@compat abstract type TriangularHarmonicPlan{T} end
+
+function *(P::TriangularHarmonicPlan, X::AbstractMatrix)
+    A_mul_B!(zero(X), P, X)
+end
+
+function \(P::TriangularHarmonicPlan, X::AbstractMatrix)
+    At_mul_B!(zero(X), P, X)
+end
+
+include("trifunctions.jl")
+include("slowplan.jl")
+
+function plan_tri2cheb(A::AbstractMatrix, α, β, γ; opts...)
+    M, N = size(A)
+    if M ≤ 1023
+        SlowTriangularHarmonicPlan(A, α, β, γ)
+    else
+        ThinTriangularHarmonicPlan(A, α, β, γ; opts...)
+    end
+end
+
+tri2cheb(A::AbstractMatrix, α, β, γ; opts...) = plan_tri2cheb(A, α, β, γ; opts...)*A
+cheb2tri(A::AbstractMatrix, α, β, γ; opts...) = plan_tri2cheb(A, α, β, γ; opts...)\A