Implement Olver's algorithm (alternative that avoids overflow)

.github/workflows/ci.yml

@@ -49,4 +49,5 @@ jobs:
       - uses: julia-actions/julia-processcoverage@v1
       - uses: codecov/codecov-action@v4
         with:
+          token: ${{ secrets.CODECOV_TOKEN }}
           file: lcov.info

Project.toml

@@ -6,7 +6,6 @@ version = "0.0.1"
 [deps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
-LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 
 [weakdeps]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
@@ -25,8 +24,9 @@ julia = "1.6"
 [extras]
 FillArrays = "1a297f60-69ca-5386-bcde-b61e274b549b"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
+LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 
 [targets]
-test = ["FillArrays", "LazyArrays", "SpecialFunctions", "Test"]
+test = ["FillArrays", "LazyArrays", "LinearAlgebra", "SpecialFunctions", "Test"]

ext/RecurrenceRelationshipsFillArraysExt.jl

@@ -1,6 +1,5 @@
 module RecurrenceRelationshipsFillArraysExt
 using RecurrenceRelationships, FillArrays
-using RecurrenceRelationships.LinearAlgebra
 using FillArrays: AbstractFill, getindex_value
 
 import RecurrenceRelationships: _forwardrecurrence_next, _clenshaw_next

ext/RecurrenceRelationshipsLazyArraysExt.jl

@@ -1,6 +1,5 @@
 module RecurrenceRelationshipsLazyArraysExt
 using RecurrenceRelationships, LazyArrays
-using RecurrenceRelationships.LinearAlgebra
 using LazyArrays.FillArrays
 
 import RecurrenceRelationships: _forwardrecurrence_next, _clenshaw_next

src/RecurrenceRelationships.jl

@@ -1,9 +1,9 @@
 module RecurrenceRelationships
-using LinearAlgebra
-export forwardrecurrence, forwardrecurrence!, clenshaw, clenshaw!
+export forwardrecurrence, forwardrecurrence!, clenshaw, clenshaw!, olver, olver!
 
 include("forward.jl")
 include("clenshaw.jl")
+include("olver.jl")
 
 if !isdefined(Base, :get_extension)
     include("../ext/RecurrenceRelationshipsFillArraysExt.jl")

src/olver.jl

@@ -0,0 +1,26 @@
+function olver!(d::AbstractVector{T}, r, a, b, c; atol=eps(T)) where T
+    resize!(r, max(length(r), 10))
+
+    r[1] = b[1]/c[1]
+    M = zero(T)
+    # d[1] is left unmodified
+    for k = 2:length(d)
+        r[k] = b[k] - c[k-1]a[k-1]/r[k-1]
+        d[k] -= a[k-1]d[k-1]/r[k-1]
+        M = max(M,one(T)) / abs(r[k])
+        if M * abs(d[k]) ≤ atol
+            resize!(d, k)
+            break
+        end
+    end
+
+    # backsubstitution
+    N = length(d)
+    d[N] /= -r[N]
+    for k = length(d)-1:-1:1
+        d[k] = (d[k]-c[k]*d[k+1])/r[k]
+    end
+    d
+end
+
+olver(d::AbstractVector{T}, a, b, c; atol=eps(T)) where T = olver!(Base.copymutable(d), similar(d), a, b, c; atol)

test/runtests.jl

@@ -1,6 +1,7 @@
-using RecurrenceRelationships, FillArrays, LazyArrays, Test
+using RecurrenceRelationships, LinearAlgebra, Test
+using FillArrays, LazyArrays
 
-@testset "clenshaw" begin
+@testset "forward/clenshaw" begin
     @testset "Chebyshev T" begin
         for elty in (Float64, Float32)
             c = [1,2,3]
@@ -136,3 +137,16 @@ using RecurrenceRelationships, FillArrays, LazyArrays, Test
         @test clenshaw(c, A, B, C, 0.1) == clenshaw(c, A, Vector(B), C, 0.1)
     end
 end
+
+@testset "olver" begin
+    N = 1000
+    x = 0.1
+    a,b,c = ones(N-1), -range(2; step=2, length=N)/x, ones(N-1)
+    j = olver([1; zeros(N-1)], a,b,c)
+    T = SymTridiagonal(Vector(b), c)
+    if VERSION ≥ v"1.10" # NoPivot doesn't exist in v1.6
+        L, U = lu(T, NoPivot())
+        n = length(j)
+        @test U[1:n,1:n] \ (L[1:n,1:n] \ [1; zeros(n-1)]) ≈ j
+    end
+end