Merge pull request #44 from MikaelSlevinsky/dl/0.7

WIP Julia v0.7-beta support

Author: dlfivefifty

.travis.yml

@@ -8,7 +8,7 @@ julia:
   - nightly
 matrix:
   allow_failures:
-    - julia: nightly
+    - julia: nightly  
 notifications:
   email: false
 after_success:

REQUIRE

@@ -1,9 +1,10 @@
 julia 0.6
-ToeplitzMatrices 0.2
+ToeplitzMatrices 0.4
 HierarchicalMatrices 0.1.1
 LowRankApprox 0.1.1
 ProgressMeter 0.3.4
 SpecialFunctions 0.3.4
 Compat 0.17
 AbstractFFTs 0.3.1
-FFTW 0.0.4
+FFTW 0.0.4
+DSP 0.4

src/FastTransforms.jl

@@ -8,21 +8,33 @@ if VERSION < v"0.7-"
     using Base.FFTW
     import Base.FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
     import Base.FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan
+    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)
 else
-    using FFTW
+    using FFTW, LinearAlgebra, DSP
     import FFTW: r2rFFTWPlan, unsafe_execute!, fftwSingle, fftwDouble, fftwNumber
-    import FFTW: libfftw, libfftwf, PlanPtr, r2rFFTWPlan
+    import FFTW: libfftw3, libfftw3f, PlanPtr, r2rFFTWPlan
+    const LAmul! = LinearAlgebra.mul!
+    const libfftw = libfftw3
+    const libfftwf = libfftw3f
+    import LinearAlgebra: Factorization
+    flipdim(A,d) = reverse(A; dims=d)
 end
 
+
 import Base: *, \, size, view
-import Base: getindex, setindex!, Factorization, length
-import Base.LinAlg: BlasFloat, BlasInt
+import Base: getindex, setindex!, length
+import Compat.LinearAlgebra: BlasFloat, BlasInt
 import HierarchicalMatrices: HierarchicalMatrix, unsafe_broadcasttimes!
-import HierarchicalMatrices: A_mul_B!, At_mul_B!, Ac_mul_B!
+import HierarchicalMatrices: mul!, At_mul_B!, Ac_mul_B!
 import HierarchicalMatrices: ThreadSafeVector, threadsafezeros
 import LowRankApprox: ColPerm
 import AbstractFFTs: Plan
-
+import Compat: range, transpose, adjoint
 
 export cjt, icjt, jjt, plan_cjt, plan_icjt
 export leg2cheb, cheb2leg, leg2chebu, ultra2ultra, jac2jac

src/PaduaTransform.jl

@@ -1,7 +1,8 @@
-doc"""
+
+# lex indicates if its lexigraphical (i.e., x, y) or reverse (y, x)
+"""
 Pre-plan an Inverse Padua Transform.
 """
-# lex indicates if its lexigraphical (i.e., x, y) or reverse (y, x)
 struct IPaduaTransformPlan{lex,IDCTPLAN,T}
     cfsmat::Matrix{T}
     idctplan::IDCTPLAN
@@ -10,37 +11,38 @@ end
 IPaduaTransformPlan(cfsmat::Matrix{T},idctplan,::Type{Val{lex}}) where {T,lex} =
     IPaduaTransformPlan{lex,typeof(idctplan),T}(cfsmat,idctplan)
 
-doc"""
+"""
 Pre-plan an Inverse Padua Transform.
 """
 function plan_ipaduatransform!(::Type{T},N::Integer,lex) where T
     n=Int(cld(-3+sqrt(1+8N),2))
     if N ≠ div((n+1)*(n+2),2)
         error("Padua transforms can only be applied to vectors of length (n+1)*(n+2)/2.")
     end
-    IPaduaTransformPlan(Array{T}(n+2,n+1),FFTW.plan_r2r!(Array{T}(n+2,n+1),FFTW.REDFT00),lex)
+    IPaduaTransformPlan(Array{T}(undef,n+2,n+1),FFTW.plan_r2r!(Array{T}(undef,n+2,n+1),FFTW.REDFT00),lex)
 end
 
 
 plan_ipaduatransform!(::Type{T},N::Integer) where {T} = plan_ipaduatransform!(T,N,Val{true})
 plan_ipaduatransform!(v::AbstractVector{T},lex...) where {T} = plan_ipaduatransform!(eltype(v),length(v),lex...)
 
+
 function *(P::IPaduaTransformPlan,v::AbstractVector{T}) where T
     cfsmat=trianglecfsmat(P,v)
     n,m=size(cfsmat)
-    scale!(view(cfsmat,:,2:m-1),0.5)
-    scale!(view(cfsmat,2:n-1,:),0.5)
+    rmul!(view(cfsmat,:,2:m-1),0.5)
+    rmul!(view(cfsmat,2:n-1,:),0.5)
     tensorvals=P.idctplan*cfsmat
     paduavec!(v,P,tensorvals)
 end
 
 ipaduatransform!(v::AbstractVector,lex...) = plan_ipaduatransform!(v,lex...)*v
-doc"""
+"""
 Inverse Padua Transform maps the 2D Chebyshev coefficients to the values of the interpolation polynomial at the Padua points.
 """
 ipaduatransform(v::AbstractVector,lex...) = plan_ipaduatransform!(v,lex...)*copy(v)
 
-doc"""
+"""
 Creates ``(n+2)x(n+1)`` Chebyshev coefficient matrix from triangle coefficients.
 """
 function trianglecfsmat(P::IPaduaTransformPlan{true},cfs::AbstractVector)
@@ -81,7 +83,7 @@ function trianglecfsmat(P::IPaduaTransformPlan{false},cfs::AbstractVector)
     return cfsmat
 end
 
-doc"""
+"""
 Vectorizes the function values at the Padua points.
 """
 function paduavec!(v,P::IPaduaTransformPlan,padmat::Matrix)
@@ -100,7 +102,7 @@ function paduavec!(v,P::IPaduaTransformPlan,padmat::Matrix)
     return v
 end
 
-doc"""
+"""
 Pre-plan a Padua Transform.
 """
 struct PaduaTransformPlan{lex,DCTPLAN,T}
@@ -111,15 +113,15 @@ end
 PaduaTransformPlan(vals::Matrix{T},dctplan,::Type{Val{lex}}) where {T,lex} =
     PaduaTransformPlan{lex,typeof(dctplan),T}(vals,dctplan)
 
-doc"""
+"""
 Pre-plan a Padua Transform.
 """
 function plan_paduatransform!(::Type{T},N::Integer,lex) where T
     n=Int(cld(-3+sqrt(1+8N),2))
     if N ≠ ((n+1)*(n+2))÷2
         error("Padua transforms can only be applied to vectors of length (n+1)*(n+2)/2.")
     end
-    PaduaTransformPlan(Array{T}(n+2,n+1),FFTW.plan_r2r!(Array{T}(n+2,n+1),FFTW.REDFT00),lex)
+    PaduaTransformPlan(Array{T}(undef,n+2,n+1),FFTW.plan_r2r!(Array{T}(undef,n+2,n+1),FFTW.REDFT00),lex)
 end
 
 plan_paduatransform!(::Type{T},N::Integer) where {T} = plan_paduatransform!(T,N,Val{true})
@@ -131,21 +133,21 @@ function *(P::PaduaTransformPlan,v::AbstractVector{T}) where T
     vals=paduavalsmat(P,v)
     tensorcfs=P.dctplan*vals
     m,l=size(tensorcfs)
-    scale!(tensorcfs,T(2)/(n*(n+1)))
-    scale!(view(tensorcfs,1,:),0.5)
-    scale!(view(tensorcfs,:,1),0.5)
-    scale!(view(tensorcfs,m,:),0.5)
-    scale!(view(tensorcfs,:,l),0.5)
+    rmul!(tensorcfs,T(2)/(n*(n+1)))
+    rmul!(view(tensorcfs,1,:),0.5)
+    rmul!(view(tensorcfs,:,1),0.5)
+    rmul!(view(tensorcfs,m,:),0.5)
+    rmul!(view(tensorcfs,:,l),0.5)
     trianglecfsvec!(v,P,tensorcfs)
 end
 
 paduatransform!(v::AbstractVector,lex...) = plan_paduatransform!(v,lex...)*v
-doc"""
+"""
 Padua Transform maps from interpolant values at the Padua points to the 2D Chebyshev coefficients.
 """
 paduatransform(v::AbstractVector,lex...) = plan_paduatransform!(v,lex...)*copy(v)
 
-doc"""
+"""
 Creates ``(n+2)x(n+1)`` matrix of interpolant values on the tensor grid at the ``(n+1)(n+2)/2`` Padua points.
 """
 function paduavalsmat(P::PaduaTransformPlan,v::AbstractVector)
@@ -165,7 +167,7 @@ function paduavalsmat(P::PaduaTransformPlan,v::AbstractVector)
     return vals
 end
 
-doc"""
+"""
 Creates length ``(n+1)(n+2)/2`` vector from matrix of triangle Chebyshev coefficients.
 """
 function trianglecfsvec!(v,P::PaduaTransformPlan{true},cfs::Matrix)
@@ -194,12 +196,12 @@ function trianglecfsvec!(v,P::PaduaTransformPlan{false},cfs::Matrix)
     return v
 end
 
-doc"""
+"""
 Returns coordinates of the ``(n+1)(n+2)/2`` Padua points.
 """
-function paduapoints(::Type{T},n::Integer) where T
+function paduapoints(::Type{T}, n::Integer) where T
     N=div((n+1)*(n+2),2)
-    MM=Matrix{T}(N,2)
+    MM=Matrix{T}(undef,N,2)
     m=0
     delta=0
     NN=fld(n+2,2)