diff --git a/docs/lit/examples/1-overview.jl b/docs/lit/examples/1-overview.jl index f429472..e0c19f8 100644 --- a/docs/lit/examples/1-overview.jl +++ b/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) diff --git a/src/strang.jl b/src/strang.jl index 228e9f8..a423be5 100644 --- a/src/strang.jl +++ b/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 diff --git a/test/strang.jl b/test/strang.jl index ccc142c..78d42f1 100644 --- a/test/strang.jl +++ b/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