Use Bidiagonal instead of BandedMatrix with bandwidth (0,1)/(1, 0) (#144)

* Bidiagonal

* Avoid [2:end]

Author: DanielVandH

Project.toml

@@ -1,7 +1,7 @@
 name = "SemiclassicalOrthogonalPolynomials"
 uuid = "291c01f3-23f6-4eb6-aeb0-063a639b53f2"
+version = "0.7.2"
 authors = ["Sheehan Olver <[email protected]>"]
-version = "0.7.1"
 
 [deps]
 ArrayLayouts = "4c555306-a7a7-4459-81d9-ec55ddd5c99a"
@@ -13,6 +13,7 @@ HypergeometricFunctions = "34004b35-14d8-5ef3-9330-4cdb6864b03a"
 InfiniteArrays = "4858937d-0d70-526a-a4dd-2d5cb5dd786c"
 InfiniteLinearAlgebra = "cde9dba0-b1de-11e9-2c62-0bab9446c55c"
 LazyArrays = "5078a376-72f3-5289-bfd5-ec5146d43c02"
+LazyBandedMatrices = "d7e5e226-e90b-4449-9968-0f923699bf6f"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 QuasiArrays = "c4ea9172-b204-11e9-377d-29865faadc5c"
 SpecialFunctions = "276daf66-3868-5448-9aa4-cd146d93841b"
@@ -27,6 +28,7 @@ HypergeometricFunctions = "0.3.4"
 InfiniteArrays = "0.15"
 InfiniteLinearAlgebra = "0.10"
 LazyArrays = "2.4"
+LazyBandedMatrices = "0.11.6"
 QuasiArrays = "0.12, 0.13"
 SpecialFunctions = "1.0, 2"
 julia = "1.10"

src/SemiclassicalOrthogonalPolynomials.jl

@@ -19,6 +19,7 @@ import ContinuumArrays: basis, Weight, @simplify, AbstractBasisLayout, BasisLayo
 import FillArrays: SquareEye
 import HypergeometricFunctions: _₂F₁general2, _₂F₁
 import SpecialFunctions: beta
+import LazyBandedMatrices
 
 export Legendre, Normalized, normalize, SemiclassicalJacobi, SemiclassicalJacobiWeight, WeightedSemiclassicalJacobi, OrthogonalPolynomialRatio
 

src/derivatives.jl

@@ -80,53 +80,46 @@ end
 function divdiff(HQ::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi}, HP::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi})
     Q = HQ.P
     P = HP.P
-    t = P.t
     a = Q.a
-    A,B,C = recurrencecoefficients(P)
-    α,β,γ = recurrencecoefficients(Q)
+    A,B,_ = recurrencecoefficients(P)
+    α,β,_ = recurrencecoefficients(Q)
     d = AccumulateAbstractVector(*, A ./ α)
-    v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
-    v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))
-
-    _BandedMatrix(
-        Vcat(
-            ((a:∞) .* v2 .- ((a+1):∞) .* Vcat(1,v1[2:end] .* d))',
-            (((a+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
+    v1 = MulAddAccumulate(Vcat(0,α[2:∞] ./ α), β)
+    v2 = MulAddAccumulate(Vcat(0,A[2:∞] ./ α), Vcat(B[1], B[2:end] .* d))
+    p = ((a+1):∞) .* v2 .- ((a+2):∞) .* v1 .* d
+    q = ((a+1):∞) .* Vcat(1,d)
+    return LazyBandedMatrices.Bidiagonal(q, p, :U)
 end
 
 function divdiff(HQ::HalfWeighted{:b,<:Any,<:SemiclassicalJacobi}, HP::HalfWeighted{:b,<:Any,<:SemiclassicalJacobi})
     Q = HQ.P
     P = HP.P
-    t = P.t
     b = Q.b
-    A,B,C = recurrencecoefficients(P)
-    α,β,γ = recurrencecoefficients(Q)
+    A,B,_ = recurrencecoefficients(P)
+    α,β,_ = recurrencecoefficients(Q)
     d = AccumulateAbstractVector(*, A ./ α)
     d2 = AccumulateAbstractVector(*, A ./ Vcat(1,α))
-    v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
-    v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))
-
-    _BandedMatrix(
-        Vcat(
-            (-(b:∞) .* v2 .+ ((b+1):∞) .* Vcat(1,v1[2:end] .* d) .+ Vcat(0,(1:∞) .* d2))',
-            (-((b+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
+    v1 = MulAddAccumulate(Vcat(0,α[2:∞] ./ α), β)
+    v2 = MulAddAccumulate(Vcat(0,A[2:∞] ./ α), Vcat(B[1], B[2:end] .* d))
+    p = -((b+1):∞) .* v2 .+ ((b+2):∞) .* v1 .* d .+ (1:∞) .* d2
+    q = -((b+1):∞) .* Vcat(1,d)
+    return LazyBandedMatrices.Bidiagonal(q, p, :U)
 end
 
 function divdiff(HQ::HalfWeighted{:c,<:Any,<:SemiclassicalJacobi}, HP::HalfWeighted{:c,<:Any,<:SemiclassicalJacobi})
     Q = HQ.P
     P = HP.P
     t = P.t
     c = Q.c
-    A,B,C = recurrencecoefficients(P)
-    α,β,γ = recurrencecoefficients(Q)
+    A,B,_ = recurrencecoefficients(P)
+    α,β,_ = recurrencecoefficients(Q)
     d = AccumulateAbstractVector(*, A ./ α)
     d2 = AccumulateAbstractVector(*, A ./ Vcat(1,α))
-    v1 = MulAddAccumulate(Vcat(0,0,α[2:∞] ./ α), Vcat(0,β))
-    v2 = MulAddAccumulate(Vcat(0,0,A[2:∞] ./ α), Vcat(0,B[1], B[2:end] .* d))
-    _BandedMatrix(
-        Vcat(
-            (-(c:∞) .* v2 .+ ((c+1):∞) .* Vcat(1,v1[2:end] .* d) .+ Vcat(0,(t:t:∞) .* d2))',
-            (-((c+1):∞) .* Vcat(1,d))'), ℵ₀, 0,1)
+    v1 = MulAddAccumulate(Vcat(0,α[2:∞] ./ α), β)
+    v2 = MulAddAccumulate(Vcat(0,A[2:∞] ./ α), Vcat(B[1], B[2:end] .* d))
+    p = -((c+1):∞) .* v2 .+ ((c+2):∞) .* v1 .* d .+ (t:t:∞) .* d2
+    q = -((c+1):∞) .* Vcat(1,d)
+    return LazyBandedMatrices.Bidiagonal(q, p, :U)
 end
 
 function diff(HP::HalfWeighted{:a,<:Any,<:SemiclassicalJacobi}; dims=1)
@@ -165,8 +158,9 @@ function divdiff(HQ::HalfWeighted{:ab}, HP::HalfWeighted{:ab})
     e = AccumulateAbstractVector(*, Vcat(1,A ./ α))
     f = MulAddAccumulate(Vcat(0,0,A[2:end] ./ α[2:end]), Vcat(0, (B./ α) .* e))
     g = cumsum(β ./ α)
-    _BandedMatrix(Vcat((((a+1):∞) .* e .- ((b+a+1):∞).*f .+ ((a+b+2):∞) .* e .* g )',
-            (-((a+b+2):∞)  .* d)'),ℵ₀,1,0)
+    p = ((a+1):∞) .* e .- ((b+a+1):∞).*f .+ ((a+b+2):∞) .* e .* g 
+    q = -((a+b+2):∞)  .* d
+    return LazyBandedMatrices.Bidiagonal(p, q, :L)
 end
 
 
@@ -182,8 +176,9 @@ function divdiff(HQ::HalfWeighted{:bc}, HP::HalfWeighted{:bc})
     e = AccumulateAbstractVector(*, Vcat(1,A ./ α))
     f = MulAddAccumulate(Vcat(0,0,A[2:end] ./ α[2:end]), Vcat(0, (B./ α) .* e))
     g = cumsum(β ./ α)
-    _BandedMatrix(Vcat((-((t+1)* (0:∞) .+ (t*(b+1) + c+1)) .* e .+ ((c+b+1):∞).*f .- ((b+c+2):∞) .* e .* g )',
-            (((b+c+2):∞)  .* d)'),ℵ₀,1,0)
+    p = -((t+1)* (0:∞) .+ (t*(b+1) + c+1)) .* e .+ ((c+b+1):∞).*f .- ((b+c+2):∞) .* e .* g 
+    q = ((b+c+2):∞)  .* d
+    return LazyBandedMatrices.Bidiagonal(p, q, :L)
 end
 
 function divdiff(HQ::HalfWeighted{:ac}, HP::HalfWeighted{:ac})
@@ -198,8 +193,9 @@ function divdiff(HQ::HalfWeighted{:ac}, HP::HalfWeighted{:ac})
     e = AccumulateAbstractVector(*, Vcat(1,A ./ α))
     f = MulAddAccumulate(Vcat(0,0,A[2:end] ./ α[2:end]), Vcat(0, (B./ α) .* e))
     g = cumsum(β ./ α)
-    _BandedMatrix(Vcat((t* ((a+1):∞) .* e .- ((c+a+1):∞).*f .+ ((a+c+2):∞) .* e .* g )',
-            (-((a+c+2):∞)  .* d)'),ℵ₀,1,0)
+    p = t* ((a+1):∞) .* e .- ((c+a+1):∞).*f .+ ((a+c+2):∞) .* e .* g 
+    q = -((a+c+2):∞)  .* d
+    return LazyBandedMatrices.Bidiagonal(p, q, :L)
 end