Merge branch 'master' of https://github.com/MikaelSlevinsky/FastTransforms.jl

Author: dlfivefifty

.travis.yml

@@ -4,10 +4,10 @@ os:
   - linux
   - osx
 julia:
-  - 0.7
   - 1.0
   - 1.1
   - 1.2    
+  - 1.3
   - nightly
 matrix:
   allow_failures:

Project.toml

@@ -1,12 +1,12 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.5"
+version = "0.6"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
-Compat = "34da2185-b29b-5c13-b0c7-acf172513d20"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
+FastGaussQuadrature = "442a2c76-b920-505d-bb47-c5924d526838"
 HierarchicalMatrices = "7c893195-952b-5c83-bb6e-be12f22eed51"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LowRankApprox = "898213cb-b102-5a47-900c-97e73b919f73"
@@ -15,13 +15,21 @@ SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
-AbstractFFTs = "≥ 0.3.1"
-Compat = "≥ 0.17.0"
-DSP = "≥ 0.4.0"
-FFTW = "≥ 0.0.4"
-HierarchicalMatrices = "≥ 0.1.1"
-LowRankApprox = "≥ 0.1.1"
-ProgressMeter = "≥ 0.3.4"
-SpecialFunctions = "≥ 0.3.4"
-ToeplitzMatrices = "≥ 0.4.0"
-julia = "0.7, 1"
+AbstractFFTs = "0.4"
+DSP = "0.6"
+FFTW = "0.3"
+HierarchicalMatrices = "0.2"
+LowRankApprox = "0.2.3"
+ProgressMeter = "1"
+SpecialFunctions = "0.8"
+ToeplitzMatrices = "0.6"
+FastGaussQuadrature = "0.4"
+julia = "1"
+
+[extras]
+Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
+Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
+Statistics = "10745b16-79ce-11e8-11f9-7d13ad32a3b2"
+
+[targets]
+test = ["Test","Random","Statistics"]

appveyor.yml

@@ -1,9 +1,9 @@
 environment:
   matrix:
-  - julia_version: 0.7
   - julia_version: 1.0
   - julia_version: 1.1
-  - julia_version: 1.2    
+  - julia_version: 1.2
+  - julia_version: 1.3
   - julia_version: nightly
 
 platform:

examples/sphere.jl

@@ -56,14 +56,14 @@ PA = FastTransforms.plan_analysis(F);
 
 # Its spherical harmonic coefficients demonstrate that it is degree-3:
 V = zero(F);
-A_mul_B!(V, PA, F);
+mul!(V, PA, F);
 U3 = P\V
 
 # Similarly, on the tensor product grid, the Legendre polynomial P₄(z⋅y) is:
 F = [P4(z(θ,φ)⋅y) for θ in θ, φ in φ]
 
 # Its spherical harmonic coefficients demonstrate that it is exact-degree-4:
-A_mul_B!(V, PA, F);
+mul!(V, PA, F);
 U4 = P\V
 
 nrm1 = vecnorm(U4);
@@ -73,7 +73,7 @@ F = [P4(z(θ,φ)⋅x) for θ in θ, φ in φ]
 
 # It only has one nonnegligible spherical harmonic coefficient.
 # Can you spot it?
-A_mul_B!(V, PA, F);
+mul!(V, PA, F);
 U4 = P\V
 
 # That nonnegligible coefficient should be approximately √(2π/(4+1/2)),

src/FastTransforms.jl

@@ -1,45 +1,30 @@
-__precompile__()
+
 module FastTransforms
 
-using ToeplitzMatrices, HierarchicalMatrices, LowRankApprox, ProgressMeter, Compat,
-        AbstractFFTs, SpecialFunctions
-
-if VERSION < v"0.7-"
-    using Base.FFTW
-    import Base.FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
-    import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan, plan_r2r!,
-                        REDFT00, REDFT01, REDFT10, REDFT11,
-                        RODFT00, RODFT01, RODFT10, RODFT11
-    const LAmul! = Base.A_mul_B!
-    import Base: Factorization
-    rmul!(A::AbstractArray, c::Number) = scale!(A,c)
-    lmul!(c::Number, A::AbstractArray) = scale!(c,A)
-    lmul!(A::AbstractArray, B::AbstractArray) = mul!(A,B)
-    rmul!(A::AbstractArray, B::AbstractArray) = mul!(A,B)
-    const floatmin = realmin
-else
-    using FFTW, LinearAlgebra, DSP
-    import FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
-    import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan, plan_r2r!,
-                    REDFT00, REDFT01, REDFT10, REDFT11,
-                    RODFT00, RODFT01, RODFT10, RODFT11
-    const LAmul! = LinearAlgebra.mul!
-    const libfftw = libfftw3
-    const libfftwf = libfftw3f
-    import LinearAlgebra: Factorization
-    flipdim(A,d) = reverse(A; dims=d)
-end
+using ToeplitzMatrices, HierarchicalMatrices, LowRankApprox, ProgressMeter,
+        AbstractFFTs, SpecialFunctions, FastGaussQuadrature
+
+using FFTW, LinearAlgebra, DSP
+import FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
+import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan, plan_r2r!,
+                REDFT00, REDFT01, REDFT10, REDFT11,
+                RODFT00, RODFT01, RODFT10, RODFT11
 
+const libfftw = libfftw3
+const libfftwf = libfftw3f
+import LinearAlgebra: Factorization
+flipdim(A,d) = reverse(A; dims=d)
 
 import Base: *, \, inv, size, view
-import Base: getindex, setindex!, length
-import Compat.LinearAlgebra: BlasFloat, BlasInt
+import Base: getindex, setindex!, length, axes
+import LinearAlgebra: BlasFloat, BlasInt, transpose, adjoint, lmul!, rmul!
 import HierarchicalMatrices: HierarchicalMatrix, unsafe_broadcasttimes!
-import HierarchicalMatrices: mul!, At_mul_B!, Ac_mul_B!
+import HierarchicalMatrices: mul!
 import HierarchicalMatrices: ThreadSafeVector, threadsafezeros
-import LowRankApprox: ColPerm
+import LowRankApprox: ColPerm, RowPerm
 import AbstractFFTs: Plan
-import Compat: range, transpose, adjoint, axes
+
+import FastGaussQuadrature: unweightedgausshermite
 
 export cjt, icjt, jjt, plan_cjt, plan_icjt
 export leg2cheb, cheb2leg, leg2chebu, ultra2ultra, jac2jac, plan_jac2jac
@@ -64,6 +49,8 @@ export SlowTriangularHarmonicPlan
 export tri2cheb, cheb2tri, plan_tri2cheb
 export triones, trizeros, trirand, trirandn, trievaluate
 
+export hermitepoints, weightedhermitetransform, iweightedhermitetransform
+
 # Other module methods and constants:
 #export ChebyshevJacobiPlan, jac2cheb, cheb2jac
 #export sqrtpi, pochhammer, stirlingseries, stirlingremainder, Aratio, Cratio, Anαβ
@@ -73,6 +60,8 @@ export triones, trizeros, trirand, trirandn, trievaluate
 #export plan_fejer2, fejernodes2, fejerweights2
 #export RecurrencePlan, forward_recurrence!, backward_recurrence
 
+lgamma(x) = logabsgamma(x)[1]
+
 include("stepthreading.jl")
 include("fftBigFloat.jl")
 include("specialfunctions.jl")
@@ -97,6 +86,7 @@ include("inufft.jl")
 include("cjt.jl")
 
 include("toeplitzhankel.jl")
+include("hermite.jl")
 
 #leg2cheb(x...)=th_leg2cheb(x...)
 #cheb2leg(x...)=th_cheb2leg(x...)

src/SphericalHarmonics/Butterfly.jl

@@ -110,10 +110,9 @@ function Butterfly(A::AbstractMatrix{T}, L::Int; isorthogonal::Bool = false, opt
     Butterfly(columns, factors, permutations, indices, threadsafezeros(T, kk), threadsafezeros(T, kk), threadsafezeros(T, kk), threadsafezeros(T, kk))
 end
 
-if VERSION ≥ v"0.7-"
-    LinearAlgebra.adjoint(B::Butterfly) = Adjoint(B)
-    LinearAlgebra.transpose(B::Butterfly) = Transpose(B)
-end
+LinearAlgebra.adjoint(B::Butterfly) = Adjoint(B)
+LinearAlgebra.transpose(B::Butterfly) = Transpose(B)
+
 
 function sumkmax(indices::Vector{Vector{Int}})
     ret = 0
@@ -173,66 +172,61 @@ function rowperm!(fwd::Bool, y::AbstractVector, x::AbstractVector, p::Vector{Int
 end
 
 ## ColumnPermutation
-mul!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = rowperm!(false, B, A.p, jstart)
-At_mul_B!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = rowperm!(true, B, A.p, jstart)
-Ac_mul_B!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = At_mul_B!(A, B, jstart)
+lmul!(A::ColPerm, B::AbstractVecOrMat, jstart::Int) = rowperm!(false, B, A.p, jstart)
+lmul!(At::RowPerm, B::AbstractVecOrMat, jstart::Int) = rowperm!(true, B, transpose(At).p, jstart)
 
 mul!(y::AbstractVector, A::ColPerm, x::AbstractVector, jstart::Int) = rowperm!(false, y, x, A.p, jstart)
-At_mul_B!(y::AbstractVector, A::ColPerm, x::AbstractVector, jstart::Int) = rowperm!(true, y, x, A.p, jstart)
-Ac_mul_B!(y::AbstractVector, A::ColPerm, x::AbstractVector, jstart::Int) = At_mul_B!(y, x, A, jstart)
+mul!(y::AbstractVector, At::RowPerm, x::AbstractVector, jstart::Int) = rowperm!(true, y, x, transpose(At).p, jstart)
 
-# Fast mul!, At_mul_B!, and Ac_mul_B! for an ID. These overwrite the output.
+# Fast mul! for an ID. These overwrite the output.
 
 
 function mul!(y::AbstractVecOrMat{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where {T}
     k, n = size(A)
-    At_mul_B!(P, x, jstart)
+    lmul!(transpose(P), x, jstart)
     copyto!(y, istart, x, jstart, k)
     mul!(y, A.T, x, istart, jstart+k)
-    mul!(P, x, jstart)
+    lmul!(P, x, jstart)
     y
 end
 
 function mul!(y::AbstractVector{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVector{T}, temp::AbstractVector{T}, istart::Int, jstart::Int) where {T}
     k, n = size(A)
-    At_mul_B!(temp, P, x, jstart)
+    mul!(temp, transpose(P), x, jstart)
     copyto!(y, istart, temp, jstart, k)
     mul!(y, A.T, temp, istart, jstart+k)
     y
 end
 
-### mul!, At_mul_B!, and  Ac_mul_B! for a Butterfly factorization.
+### mul! for a Butterfly factorization.
 mul!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)
 
-for f! in (:At_mul_B!, :Ac_mul_B!)
+for Trans in (:Transpose, :Adjoint)
     @eval begin
-        function $f!(y::AbstractVecOrMat{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where {T}
+        function mul!(y::AbstractVecOrMat{T}, Ac::$Trans{T,IDPackedV{T}}, P::ColumnPermutation, x::AbstractVecOrMat{T}, istart::Int, jstart::Int) where {T}
+            A = parent(Ac)
             k, n = size(A)
             copyto!(y, istart, x, jstart, k)
-            $f!(y, A.T, x, istart+k, jstart)
-            mul!(P, y, istart)
+            mul!(y, $Trans(A.T), x, istart+k, jstart)
+            lmul!(P, y, istart)
             y
         end
 
-        function $f!(y::AbstractVector{T}, A::IDPackedV{T}, P::ColumnPermutation, x::AbstractVector{T}, temp::AbstractVector{T}, istart::Int, jstart::Int) where {T}
+        function mul!(y::AbstractVector{T}, Ac::$Trans{T,IDPackedV{T}}, P::ColumnPermutation, x::AbstractVector{T}, temp::AbstractVector{T}, istart::Int, jstart::Int) where {T}
+            A = parent(Ac)
             k, n = size(A)
             copyto!(temp, istart, x, jstart, k)
-            $f!(temp, A.T, x, istart+k, jstart)
+            mul!(temp, $Trans(A.T), x, istart+k, jstart)
             mul!(y, P, temp, istart)
             y
         end
     end
 end
 
-if VERSION < v"0.7-"
-    Base.A_mul_B!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)
-    Base.At_mul_B!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = At_mul_B_col_J!(u, B, b, 1)
-    Base.Ac_mul_B!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = Ac_mul_B_col_J!(u, B, b, 1)
-else
-    LinearAlgebra.mul!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)
-    LinearAlgebra.mul!(u::Vector{T}, Bt::Transpose{T,Butterfly{T}}, b::Vector{T}) where T = At_mul_B_col_J!(u, parent(Bt), b, 1)
-    LinearAlgebra.mul!(u::Vector{T}, Bc::Adjoint{T,Butterfly{T}}, b::Vector{T}) where T = Ac_mul_B_col_J!(u, parent(Bc), b, 1)
-end
+LinearAlgebra.mul!(u::Vector{T}, B::Butterfly{T}, b::Vector{T}) where T = mul_col_J!(u, B, b, 1)
+LinearAlgebra.mul!(u::Vector{T}, Bt::Transpose{T,Butterfly{T}}, b::Vector{T}) where T = mul_col_J!(u, Bt, b, 1)
+LinearAlgebra.mul!(u::Vector{T}, Bc::Adjoint{T,Butterfly{T}}, b::Vector{T}) where T = mul_col_J!(u, Bc, b, 1)
+
 
 
 function mul_col_J!(u::VecOrMat{T}, B::Butterfly{T}, b::VecOrMat{T}, J::Int) where T
@@ -290,10 +284,10 @@ function mul_col_J!(u::VecOrMat{T}, B::Butterfly{T}, b::VecOrMat{T}, J::Int) whe
     u
 end
 
-for f! in (:At_mul_B!,:Ac_mul_B!)
-    f_col_J! = Meta.parse(string(f!)[1:end-1]*"_col_J!")
+for Trans in (:Transpose,:Adjoint)
     @eval begin
-        function $f_col_J!(u::VecOrMat{T}, B::Butterfly{T}, b::VecOrMat{T}, J::Int) where T
+        function mul_col_J!(u::VecOrMat{T}, Bc::$Trans{T,Butterfly{T}}, b::VecOrMat{T}, J::Int) where T
+            B = parent(Bc)
             L = length(B.factors) - 1
             tL = 2^L
 
@@ -315,7 +309,7 @@ for f! in (:At_mul_B!,:Ac_mul_B!)
             for i = 1:tL
                 ml = mu+1
                 mu += size(columns[i], 1)
-                $f!(temp1, columns[i], b, inds[i], ml+COLSHIFT)
+                mul!(temp1, $Trans(columns[i]), b, inds[i], ml+COLSHIFT)
             end
 
             ii, jj = tL, 1
@@ -329,7 +323,7 @@ for f! in (:At_mul_B!,:Ac_mul_B!)
                     shft = 2jj*div(ctr,2jj)
                     fill!(temp4, zero(T))
                     for j = 1:jj
-                        $f!(temp3, factors[j+ctr], permutations[j+ctr], temp1, temp4, indsout[2j+shft-1], indsin[j+ctr])
+                        mul!(temp3, $Trans(factors[j+ctr]), permutations[j+ctr], temp1, temp4, indsout[2j+shft-1], indsin[j+ctr])
                         addtemp3totemp2!(temp2, temp3, indsout[2j+shft-1], indsout[2j+shft+1]-1)
                     end
                     ctr += jj
@@ -346,7 +340,7 @@ for f! in (:At_mul_B!,:Ac_mul_B!)
             for j = 1:tL
                 nl = nu+1
                 nu += size(factors[j], 2)
-                $f!(u, factors[j], permutations[j], temp1, nl+COLSHIFT, inds[j])
+                mul!(u, $Trans(factors[j]), permutations[j], temp1, nl+COLSHIFT, inds[j])
             end
 
             u

src/SphericalHarmonics/SphericalHarmonics.jl

@@ -4,14 +4,9 @@ function *(P::SphericalHarmonicPlan, X::AbstractMatrix)
     mul!(zero(X), P, X)
 end
 
-if VERSION < v"0.7-"
-    function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
-        At_mul_B!(zero(X), P, X)
-    end
-else
-    function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
-        mul!(zero(X), transpose(P), X)
-    end
+
+function \(P::SphericalHarmonicPlan, X::AbstractMatrix)
+    mul!(zero(X), transpose(P), X)
 end
 
 include("sphfunctions.jl")