Add LDLt factorization for Strang (#78)

* Add factorize for Strang

* Show Strang LDL'

* Limit factorize to >= v1.8

docs/lit/examples/1-overview.jl

@@ -23,6 +23,7 @@ Toeplitz, Hankel, and circulant matrices.
 
 using SpecialMatrices
 using Polynomials
+using LinearAlgebra: factorize, Diagonal
 
 
 # ## [`Cauchy` matrix](http://en.wikipedia.org/wiki/Cauchy_matrix)
@@ -111,7 +112,14 @@ Riemann(5)
 A special symmetric, tridiagonal, Toeplitz matrix named after Gilbert Strang.
 =#
 
-Strang(5)
+S = Strang(5)
+
+# The Strang matrix has a special ``L D L'`` factorization:
+F = factorize(S)
+
+# Here is a verification:
+F.L * Diagonal(F.D) * F.L'
+
 
 # ## [`Vandermonde` matrix](http://en.wikipedia.org/wiki/Vandermonde_matrix)
 

src/strang.jl

@@ -1,5 +1,7 @@
 export Strang
 
+using LinearAlgebra: LDLt, SymTridiagonal
+
 """
     Strang([T=Int,] n::Int)
 
@@ -61,3 +63,29 @@ end
     y[2:end] .-= @view x[1:end-1]
     return y
 end
+
+
+# LU factors for future reference:
+#=
+SparseArrays: spdiagm
+   L = spdiagm(0 => ones(n), -1 => -(1:n-1) ./ (2:n))
+   U = spdiagm(0 => (2:n+1) ./ (1:n), 1 => -ones(n-1))
+=#
+
+if VERSION >= v"1.8"
+#=
+Strang is a special case of SymTridiagonal,
+and the default factorization of that class is LDLt
+https://docs.julialang.org/en/v1/stdlib/LinearAlgebra/#LinearAlgebra.factorize
+so we use that here,
+even though it does not fully exploit the sparse structure
+of the LL' and LDL' factorizations of a Strang matrix.
+=#
+function LinearAlgebra.factorize(A::Strang)
+    n = A.n
+    ev = -(1:n-1) ./ (2:n) # L = spdiagm(0 => ones(n), -1 => ev)
+    dv = (2:n+1) ./ (1:n) # D = Diagonal(dv)
+    S = SymTridiagonal(dv, ev)
+    return LDLt(S)
+end
+end

test/strang.jl

@@ -2,6 +2,7 @@
 
 using SpecialMatrices: Strang
 using LinearAlgebra: diagm, issymmetric, ishermitian, isposdef
+using LinearAlgebra: factorize, SymTridiagonal, Diagonal
 using Test: @test, @testset, @test_throws, @inferred
 
 @testset "strang" begin
@@ -30,3 +31,22 @@ using Test: @test, @testset, @test_throws, @inferred
     @test ishermitian(A)
     @test isposdef(A)
 end
+
+
+if VERSION >= v"1.8"
+
+@testset "strang-factor" begin
+    # factorize:
+    n = 5
+    A = Strang(5)
+    M = Matrix(A)
+    T = SymTridiagonal(M)
+    F = factorize(T)
+    G = factorize(A)
+    H = G.L * Diagonal(G.D) * G.L'
+    @test H ≈ M
+    @test F.L ≈ G.L
+    @test F.D ≈ G.D
+end
+
+end