Support clenshaw! with any DenseColumnMajor blas vector (#113)

* Support clenshaw! with any DenseColumnMajor blas vector

* Fix out of bounds error

* reactivate coverage

* v0.9.3

* Allow strided coefficients

* fix c-call

* Test Zeros diagonal special case

* Update clenshawtests.jl

Co-authored-by: Mikael Slevinsky <Richard.Slevinsky@umanitoba.ca>

Author: dlfivefifty

.github/workflows/ci.yml

@@ -41,9 +41,10 @@ jobs:
             ${{ runner.os }}-
       - uses: julia-actions/julia-buildpkg@latest
       - uses: julia-actions/julia-runtest@latest
-      - uses: julia-actions/julia-uploadcodecov@latest
-        env:
-          CODECOV_TOKEN: ${{ secrets.CODECOV_TOKEN }}
+      - uses: julia-actions/julia-processcoverage@v1
+      - uses: codecov/codecov-action@v1
+        with:
+          file: lcov.info
   docs:
     name: Documentation
     runs-on: ubuntu-latest

Project.toml

@@ -1,9 +1,10 @@
 name = "FastTransforms"
 uuid = "057dd010-8810-581a-b7be-e3fc3b93f78c"
-version = "0.9.2"
+version = "0.9.3"
 
 [deps]
 AbstractFFTs = "621f4979-c628-5d54-868e-fcf4e3e8185c"
+ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
 BinaryProvider = "b99e7846-7c00-51b0-8f62-c81ae34c0232"
 DSP = "717857b8-e6f2-59f4-9121-6e50c889abd2"
 FFTW = "7a1cc6ca-52ef-59f5-83cd-3a7055c09341"
@@ -19,6 +20,7 @@ ToeplitzMatrices = "c751599d-da0a-543b-9d20-d0a503d91d24"
 
 [compat]
 AbstractFFTs = "0.4, 0.5"
+ArrayLayouts = "0.3.7"
 BinaryProvider = "0.5"
 DSP = "0.6"
 FFTW = "1"

README.md

@@ -1,6 +1,6 @@
 # FastTransforms.jl
 
-[![Build Status](https://github.com/JuliaApproximation/FastTransforms.jl/workflows/CI/badge.svg)](https://github.com/JuliaApproximation/FastTransforms.jl/actions?query=workflow%3ACI) [![Travis](https://travis-ci.org/JuliaApproximation/FastTransforms.jl.svg?branch=master)](https://travis-ci.org/JuliaApproximation/FastTransforms.jl) [![codecov](https://codecov.io/gh/JuliaApproximation/FastTransforms.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaApproximation/FastTransforms.jl) [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaApproximation.github.io/FastTransforms.jl/stable) [![](https://img.shields.io/badge/docs-latest-blue.svg)](https://JuliaApproximation.github.io/FastTransforms.jl/latest)
+[![Build Status](https://github.com/JuliaApproximation/FastTransforms.jl/workflows/CI/badge.svg)](https://github.com/JuliaApproximation/FastTransforms.jl/actions?query=workflow%3ACI) [![Travis](https://travis-ci.org/JuliaApproximation/FastTransforms.jl.svg?branch=master)](https://travis-ci.org/JuliaApproximation/FastTransforms.jl) [![codecov](https://codecov.io/gh/JuliaApproximation/FastTransforms.jl/branch/master/graph/badge.svg)](https://codecov.io/gh/JuliaApproximation/FastTransforms.jl) [![](https://img.shields.io/badge/docs-stable-blue.svg)](https://JuliaApproximation.github.io/FastTransforms.jl/stable) [![](https://img.shields.io/badge/docs-dev-blue.svg)](https://JuliaApproximation.github.io/FastTransforms.jl/dev)
 
 `FastTransforms.jl` allows the user to conveniently work with orthogonal polynomials with degrees well into the millions.
 

src/FastTransforms.jl

@@ -1,7 +1,7 @@
 module FastTransforms
 
 using FastGaussQuadrature, LinearAlgebra
-using Reexport, SpecialFunctions, ToeplitzMatrices, FillArrays
+using Reexport, SpecialFunctions, ToeplitzMatrices, FillArrays, ArrayLayouts
 
 import DSP
 
@@ -48,6 +48,8 @@ export plan_leg2cheb, plan_cheb2leg, plan_ultra2ultra, plan_jac2jac,
        plan_tet2cheb, plan_tet_synthesis, plan_tet_analysis,
        plan_spinsph2fourier, plan_spinsph_synthesis, plan_spinsph_analysis
 
+include("clenshaw.jl")
+
 include("libfasttransforms.jl")
 
 export plan_nufft, plan_nufft1, plan_nufft2, plan_nufft3, plan_inufft1, plan_inufft2
@@ -96,6 +98,6 @@ lgamma(x) = logabsgamma(x)[1]
 
 include("specialfunctions.jl")
 
-include("clenshaw.jl")
+
 
 end # module

src/clenshaw.jl

@@ -68,10 +68,14 @@ function clenshaw!(c::AbstractVector, A::AbstractVector, B::AbstractVector, C::A
 end
 
 
-@inline _clenshaw_next(n, A, B, C, x, c, bn1, bn2) = muladd(muladd(A[n],x,B[n]), bn1, muladd(-C[n+1],bn2,c[n]))
-@inline _clenshaw_next(n, A, ::Zeros, C, x, c, bn1, bn2) = muladd(A[n]*x, bn1, muladd(-C[n+1],bn2,c[n]))
+Base.@propagate_inbounds _clenshaw_next(n, A, B, C, x, c, bn1, bn2) = muladd(muladd(A[n],x,B[n]), bn1, muladd(-C[n+1],bn2,c[n]))
+Base.@propagate_inbounds _clenshaw_next(n, A, ::Zeros, C, x, c, bn1, bn2) = muladd(A[n]*x, bn1, muladd(-C[n+1],bn2,c[n]))
 # Chebyshev U
-@inline _clenshaw_next(n, A::AbstractFill, ::Zeros, C::Ones, x, c, bn1, bn2) = muladd(getindex_value(A)*x, bn1, -bn2+c[n])
+Base.@propagate_inbounds _clenshaw_next(n, A::AbstractFill, ::Zeros, C::Ones, x, c, bn1, bn2) = muladd(getindex_value(A)*x, bn1, -bn2+c[n])
+
+# allow special casing first arg, for ChebyshevT in OrthogonalPolynomialsQuasi
+Base.@propagate_inbounds _clenshaw_first(A, B, C, x, c, bn1, bn2) = muladd(muladd(A[1],x,B[1]), bn1, muladd(-C[2],bn2,c[1]))
+
 
 """
     clenshaw(c, A, B, C, x)
@@ -90,9 +94,11 @@ function clenshaw(c::AbstractVector, A::AbstractVector, B::AbstractVector, C::Ab
     @inbounds begin
         bn2 = zero(T)
         bn1 = convert(T,c[N])
-        for n = N-1:-1:1
+        N == 1 && return bn1
+        for n = N-1:-1:2
             bn1,bn2 = _clenshaw_next(n, A, B, C, x, c, bn1, bn2),bn1
         end
+        bn1 = _clenshaw_first(A, B, C, x, c, bn1, bn2)
     end
     bn1
 end
@@ -120,7 +126,9 @@ clenshaw!(c::AbstractVector, x::AbstractVector) = clenshaw!(c, x, x)
 evaluates the first-kind Chebyshev (T) expansion with coefficients `c` at points `x`,
 overwriting `f` with the results.
 """
-function clenshaw!(c::AbstractVector, x::AbstractVector, f::AbstractVector)
+clenshaw!(c::AbstractVector, x::AbstractVector, f::AbstractVector) = _clenshaw!(MemoryLayout(c), MemoryLayout(x), MemoryLayout(f), c, x, f)
+
+function _clenshaw!(_, _, _, c::AbstractVector, x::AbstractVector, f::AbstractVector)
     f .= clenshaw.(Ref(c), x)
 end
 

src/libfasttransforms.jl

@@ -73,15 +73,19 @@ function check_clenshaw_points(x, ϕ₀, f)
     length(x) == length(ϕ₀) == length(f) || throw(ArgumentError("Dimensions must match"))
 end
 
-function clenshaw!(c::Vector{Float64}, x::Vector{Float64}, f::Vector{Float64})
-    @assert length(x) == length(f)
-    ccall((:ft_clenshaw, libfasttransforms), Cvoid, (Cint, Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Float64}), length(c), c, 1, length(x), x, f)
+function check_clenshaw_points(x, f)
+    length(x) == length(f) || throw(ArgumentError("Dimensions must match"))
+end
+
+function _clenshaw!(::AbstractStridedLayout, ::AbstractColumnMajor, ::AbstractColumnMajor, c::AbstractVector{Float64}, x::AbstractVector{Float64}, f::AbstractVector{Float64})
+    @boundscheck check_clenshaw_points(x, f)
+    ccall((:ft_clenshaw, libfasttransforms), Cvoid, (Cint, Ptr{Float64}, Cint, Cint, Ptr{Float64}, Ptr{Float64}), length(c), c, stride(c,1), length(x), x, f)
     f
 end
 
-function clenshaw!(c::Vector{Float32}, x::Vector{Float32}, f::Vector{Float32})
-    @assert length(x) == length(f)
-    ccall((:ft_clenshawf, libfasttransforms), Cvoid, (Cint, Ptr{Float32}, Cint, Cint, Ptr{Float32}, Ptr{Float32}), length(c), c, 1, length(x), x, f)
+function _clenshaw!(::AbstractStridedLayout, ::AbstractColumnMajor, ::AbstractColumnMajor, c::AbstractVector{Float32}, x::AbstractVector{Float32}, f::AbstractVector{Float32})
+    @boundscheck check_clenshaw_points(x, f)
+    ccall((:ft_clenshawf, libfasttransforms), Cvoid, (Cint, Ptr{Float32}, Cint, Cint, Ptr{Float32}, Ptr{Float32}), length(c), c, stride(c,1), length(x), x, f)
     f
 end
 

test/clenshawtests.jl

@@ -3,20 +3,34 @@ import FastTransforms: clenshaw, clenshaw!, forwardrecurrence!, forwardrecurrenc
 
 @testset "clenshaw" begin
     @testset "Chebyshev T" begin
-        c = [1,2,3]
-        cf = float(c)
-        @test @inferred(clenshaw(c,1)) ≡ 1 + 2 + 3
-        @test @inferred(clenshaw(c,0)) ≡ 1 + 0 - 3
-        @test @inferred(clenshaw(c,0.1)) == 1 + 2*0.1 + 3*cos(2acos(0.1))
-        @test @inferred(clenshaw(c,[-1,0,1])) == clenshaw!(c,[-1,0,1]) == [2,-2,6]
-        @test clenshaw(c,[-1,0,1]) isa Vector{Int}
-        @test @inferred(clenshaw(Float64[],1)) ≡ 0.0
+        for elty in (Float64, Float32)
+            c = [1,2,3]
+            cf = elty.(c)
+            @test @inferred(clenshaw(c,1)) ≡ 1 + 2 + 3
+            @test @inferred(clenshaw(c,0)) ≡ 1 + 0 - 3
+            @test @inferred(clenshaw(c,0.1)) == 1 + 2*0.1 + 3*cos(2acos(0.1))
+            @test @inferred(clenshaw(c,[-1,0,1])) == clenshaw!(c,[-1,0,1]) == [2,-2,6]
+            @test clenshaw(c,[-1,0,1]) isa Vector{Int}
+            @test @inferred(clenshaw(elty[],1)) ≡ zero(elty)
 
-        x = [1,0,0.1]
-        @test @inferred(clenshaw(c,x)) ≈ @inferred(clenshaw!(c,copy(x))) ≈ 
-            @inferred(clenshaw!(c,x,similar(x))) ≈
-            @inferred(clenshaw(cf,x)) ≈ @inferred(clenshaw!(cf,copy(x))) ≈ 
-            @inferred(clenshaw!(cf,x,similar(x))) ≈ [6,-2,-1.74]
+            x = elty[1,0,0.1]
+            @test @inferred(clenshaw(c,x)) ≈ @inferred(clenshaw!(c,copy(x))) ≈ 
+                @inferred(clenshaw!(c,x,similar(x))) ≈
+                @inferred(clenshaw(cf,x)) ≈ @inferred(clenshaw!(cf,copy(x))) ≈ 
+                @inferred(clenshaw!(cf,x,similar(x))) ≈ elty[6,-2,-1.74]
+
+            @testset "Strided" begin
+                cv = view(cf,:)
+                xv = view(x,:)
+                @test clenshaw!(cv, xv, similar(xv)) == clenshaw!(cf,x,similar(x))
+
+                cv2 = view(cf,1:2:3)
+                @test clenshaw!(cv2, xv, similar(xv)) == clenshaw([1,3], x)
+
+                # modifies x and xv
+                @test clenshaw!(cv2, xv) == xv == x == clenshaw([1,3], elty[1,0,0.1])
+            end
+        end
     end
 
     @testset "Chebyshev U" begin
@@ -101,6 +115,13 @@ import FastTransforms: clenshaw, clenshaw!, forwardrecurrence!, forwardrecurrenc
         @test v_f isa Vector{Float64}
 
         j = 3
-        clenshaw([zeros(Int,j-1); 1; zeros(Int,N-j)], A, B, C, 1) == v_i[j]
+        @test clenshaw([zeros(Int,j-1); 1; zeros(Int,N-j)], A, B, C, 1) == v_i[j]
+    end
+
+    @testset "Zeros diagonal" begin
+        N = 10; A = randn(N); B = Zeros{Int}(N); C = randn(N+1)
+        @test forwardrecurrence(N, A, B, C, 0.1) == forwardrecurrence(N, A, Vector(B), C, 0.1)
+        c = randn(N)
+        @test clenshaw(c, A, B, C, 0.1) == clenshaw(c, A, Vector(B), C, 0.1)
     end
 end