Speed up log kernel for functions (#33)

* Speed up log kernel for functions

* real-valued log

* add test for sub-Legendre

Project.toml

@@ -1,7 +1,7 @@
 name = "SingularIntegrals"
 uuid = "d7440221-8b5e-42fc-909c-0567823f424a"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.2.2"
+version = "0.2.3"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -19,16 +19,16 @@ QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 
 [compat]
-ArrayLayouts = "1.0.11"
+ArrayLayouts = "1.4"
 BandedMatrices = "0.17.17"
 ClassicalOrthogonalPolynomials = "0.11"
-ContinuumArrays = "0.14, 0.15, 0.16"
+ContinuumArrays = "0.16"
 FastTransforms = "0.15"
 FillArrays = "1"
 HypergeometricFunctions = "0.3.4"
-InfiniteArrays = "0.12.11, 0.13"
-LazyArrays = "1"
-LazyBandedMatrices = "0.8.5, 0.9"
+InfiniteArrays = "0.13"
+LazyArrays = "1.8"
+LazyBandedMatrices = "0.9"
 QuasiArrays = "0.11"
 SpecialFunctions = "1, 2"
 julia = "1.9"

src/SingularIntegrals.jl

@@ -1,11 +1,11 @@
 module SingularIntegrals
 using ClassicalOrthogonalPolynomials, ContinuumArrays, QuasiArrays, LazyArrays, LazyBandedMatrices, FillArrays, BandedMatrices, LinearAlgebra, SpecialFunctions, HypergeometricFunctions, InfiniteArrays
-using ContinuumArrays: @simplify, Weight, AbstractAffineQuasiVector, inbounds_getindex, broadcastbasis, MappedBasisLayouts, MemoryLayout, MappedWeightLayout, ExpansionLayout, demap, basismap, AbstractBasisLayout
+using ContinuumArrays: @simplify, Weight, AbstractAffineQuasiVector, inbounds_getindex, broadcastbasis, MappedBasisLayouts, MemoryLayout, MappedWeightLayout, ExpansionLayout, demap, basismap, AbstractBasisLayout, SubBasisLayout
 using QuasiArrays: AbstractQuasiMatrix, BroadcastQuasiMatrix, LazyQuasiArrayStyle, AbstractQuasiVecOrMat
 import ClassicalOrthogonalPolynomials: AbstractJacobiWeight, WeightedBasis, jacobimatrix, orthogonalityweight, recurrencecoefficients, _p0, Clenshaw, chop, initiateforwardrecurrence, MappedOPLayouts, unweighted
 using LazyBandedMatrices: Tridiagonal, SymTridiagonal, subdiagonaldata, supdiagonaldata, diagonaldata, ApplyLayout
 import LazyArrays: AbstractCachedMatrix, AbstractCachedArray, paddeddata, arguments, resizedata!, cache_filldata!, zero!, cacheddata
-import Base: *, +, -, /, \, Slice, axes, getindex, sum, ==, oneto, size, broadcasted, copy, tail
+import Base: *, +, -, /, \, Slice, axes, getindex, sum, ==, oneto, size, broadcasted, copy, tail, view
 import LinearAlgebra: dot
 using BandedMatrices: _BandedMatrix
 using FastTransforms: _forwardrecurrence!, _forwardrecurrence_next

src/logkernel.jl

@@ -52,6 +52,7 @@ for lk in (:logkernel, :complexlogkernel)
         end
 
         $lk_layout(::ExpansionLayout, A, dims...) = $lk_layout(ApplyLayout{typeof(*)}(), A, dims...)
+        $lk_layout(::SubBasisLayout, A, dims...) = $lk(parent(A), dims...)[:, parentindices(A)[2]]
     end
 end
 
@@ -89,7 +90,7 @@ end
 ####
 
 
-function logkernel(wP::Weighted{T,<:ChebyshevU}, z) where T
+function logkernel(wP::Weighted{T,<:ChebyshevU}, z::Number) where T
     if z in axes(wP,1)
         Tn = Vcat(convert(T,π)*log(2*one(T)), convert(T,π)*ChebyshevT{T}()[z,2:end]./oneto(∞))
         return transpose((Tn[3:end]-Tn[1:end])/2)
@@ -106,15 +107,57 @@ function logkernel(wP::Weighted{T,<:ChebyshevU}, z) where T
 
 end
 
-function complexlogkernel(P::Legendre{T}, z) where T
+function complexlogkernel(P::Legendre, z::Number)
+    T = promote_type(eltype(P), typeof(z))
     r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T)
     r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1)
     r2 = z*r1 + 2*one(T)/3
     transpose(RecurrenceArray(z, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), [r0,r1,r2]))
 end
 
+function complexlogkernel(P::Legendre, zs::AbstractVector)
+    T = promote_type(eltype(P), eltype(zs))
+    m = length(zs)
+    data = Matrix{T}(undef, 3, m)
+    for j = 1:m
+        z = zs[j]
+        r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T)
+        r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1)
+        r2 = z*r1 + 2*one(T)/3
+        data[1,j] = r0; data[2,j] = r1; data[3,j] = r2;
+    end
+
+    transpose(RecurrenceArray(zs, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), data))
+end
+
+
 logkernel(P::Legendre, z) = real.(complexlogkernel(P, complex(z)))
 
+function logkernel(P::Legendre, x::Real)
+    T = promote_type(eltype(P), typeof(x))
+    z = complex(x)
+    r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T)
+    r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1)
+    r2 = z*r1 + 2*one(T)/3
+    transpose(RecurrenceArray(x, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), [real(r0),real(r1),real(r2)]))
+end
+
+function logkernel(P::Legendre, x::AbstractVector{<:Real})
+    T = promote_type(eltype(P), eltype(x))
+    m = length(x)
+    data = Matrix{T}(undef, 3, m)
+    for j = 1:m
+        z = complex(x[j])
+        r0 = (1 + z)log(1 + z) - (z-1)log(z-1) - 2one(T)
+        r1 = (z+1)*r0/2 + 1 - (z+1)log(z+1)
+        r2 = z*r1 + 2*one(T)/3
+        data[1,j] = real(r0); data[2,j] = real(r1); data[3,j] = real(r2);
+    end
+
+    transpose(RecurrenceArray(x, ((one(real(T)):2:∞)./(2:∞), Zeros{real(T)}(∞), (-one(real(T)):∞)./(2:∞)), data))
+end
+
+
 ###
 # Maps
 ###

src/recurrence.jl

@@ -154,10 +154,15 @@ end
 
 _getindex_resize_iffinite!(A, kr, jr, _) = layout_getindex(A, kr, jr)
 
-@inline getindex(A::RecurrenceMatrix, kr::AbstractUnitRange, jr::AbstractUnitRange) = _getindex_resize_iffinite!(A, kr, jr, last(kr))
-@inline getindex(A::RecurrenceMatrix, kr::AbstractVector, jr::AbstractVector) = _getindex_resize_iffinite!(A, kr, jr, last(kr))
+@inline getindex(A::RecurrenceMatrix, kr::AbstractUnitRange, jr::AbstractUnitRange) = _getindex_resize_iffinite!(A, kr, jr, maximum(kr))
+@inline getindex(A::RecurrenceMatrix, kr::AbstractVector, jr::AbstractVector) = _getindex_resize_iffinite!(A, kr, jr, maximum(kr))
 @inline getindex(A::RecurrenceMatrix, k::Integer, jr::AbstractVector) = _getindex_resize_iffinite!(A, k, jr, k)
 @inline getindex(A::RecurrenceMatrix, k::Integer, jr::AbstractUnitRange) = _getindex_resize_iffinite!(A, k, jr, k)
 @inline getindex(A::RecurrenceMatrix, k::Integer, ::Colon) = _getindex_resize_iffinite!(A, k, :, k)
-@inline getindex(A::RecurrenceMatrix, kr::AbstractVector, ::Colon) = _getindex_resize_iffinite!(A, kr, :, last(kr))
-@inline getindex(A::RecurrenceMatrix, kr::AbstractUnitRange, ::Colon) = _getindex_resize_iffinite!(A, kr, :, last(kr))
+@inline getindex(A::RecurrenceMatrix, kr::AbstractVector, ::Colon) = _getindex_resize_iffinite!(A, kr, :, maximum(kr))
+@inline getindex(A::RecurrenceMatrix, kr::AbstractUnitRange, ::Colon) = _getindex_resize_iffinite!(A, kr, :, maximum(kr))
+
+function view(A::RecurrenceVector, kr::AbstractVector)
+    resizedata!(A, maximum(kr))
+    view(A.data, kr)
+end

test/test_logkernel.jl

@@ -140,4 +140,18 @@ end
         L = log.(abs.(x .- x'))
         @test T[0.2,:]'*((T\L*Weighted(T)) * (T\exp.(x))) ≈ -2.9976362326874373 # Mathematica
     end
+
+    @testset "Legendre" begin
+        P = Legendre()
+        f = expand(P, exp)
+        @test logkernel(f, 0.1) ≈ logkernel(f, complex(0.1)) ≈ -2.3204982810441956
+        @test logkernel(f, [0.1,1.2,-0.5]) ≈ logkernel(f, ComplexF64[0.1,1.2,-0.5])
+
+        @test logkernel(f, 0.1 + 0.2im) ≈ -1.6570185704416018
+    end
+
+    @testset "sub-Legendre" begin
+        P = Legendre()
+        @test logkernel(P[:,1:10],0.1)' == logkernel(P, 0.1)'[1:10]
+    end
 end