From a0fd83b9d3f74313993549d00f1d053556a11e3c Mon Sep 17 00:00:00 2001 From: Jiahao Chen Date: Sat, 11 Jan 2014 01:17:09 -0500 Subject: [PATCH] Provide inv and det of Tridiagonal and SymTridiagonal matrices --- base/linalg/tridiag.jl | 111 +++++++++++++++++++++++++++++++++++++++++ test/linalg.jl | 34 +++++++++++++ 2 files changed, 145 insertions(+) diff --git a/base/linalg/tridiag.jl b/base/linalg/tridiag.jl index 8d84655705590..5899a7afb5693 100644 --- a/base/linalg/tridiag.jl +++ b/base/linalg/tridiag.jl @@ -89,6 +89,65 @@ eigmin(m::SymTridiagonal) = eigvals(m, 1, 1)[1] eigvecs(m::SymTridiagonal) = eig(m)[2] eigvecs{Eigenvalue<:Real}(m::SymTridiagonal, eigvals::Vector{Eigenvalue}) = LAPACK.stein!(m.dv, m.ev, eigvals) +################### +# Generic methods # +################### + +#Needed for inv() +type ZeroOffsetVector + data::Vector +end +getindex (a::ZeroOffsetVector, i) = a.data[i+1] +setindex!(a::ZeroOffsetVector, x, i) = a.data[i+1]=x + +#Implements the inverse using principal minors +#Ref: +# R. Usmani, "Inversion of a tridiagonal Jacobi matrix", +# Linear Algebra and its Applications 212-213 (1994), pp.413-414 +# doi:10.1016/0024-3795(94)90414-6 +function inv{T}(A::SymTridiagonal{T}) + n = size(A, 1) + θ = ZeroOffsetVector(zeros(T, n+1)) #principal minors of A + θ[0] = 1 + n>=1 && (θ[1] = A.dv[1]) + for i=2:n + θ[i] = A.dv[i]*θ[i-1]-A.ev[i-1]^2*θ[i-2] + end + φ = zeros(T, n+1) + φ[n+1] = 1 + n>=1 && (φ[n] = A.dv[n]) + for i=n-1:-1:1 + φ[i] = A.dv[i]*φ[i+1]-A.ev[i]^2*φ[i+2] + end + α = Array(T, n, n) + for i=1:n, j=1:n + sign = (i+j)%2==0 ? (+) : (-) + if ij + α[i,j]=(sign)(prod(A.ev[j:i-1]))*θ[j-1]*φ[i+1]/θ[n] + end + end + α +end + +#Implements the determinant using principal minors +# R. Usmani, "Inversion of a tridiagonal Jacobi matrix", +# Linear Algebra and its Applications 212-213 (1994), pp.413-414 +# doi:10.1016/0024-3795(94)90414-6 +function det{T}(A::SymTridiagonal{T}) + n = size(A, 1) + θa = one(T) + n==0 && return θa + θb = A.dv[1] + for i=2:n + θb, θa = A.dv[i]*θb-A.ev[i-1]^2*θa, θb + end + return θb +end + ## Tridiagonal matrices ## type Tridiagonal{T} <: AbstractMatrix{T} dl::Vector{T} # sub-diagonal @@ -175,6 +234,10 @@ function diag{T}(M::Tridiagonal{T}, n::Integer=0) end end +################### +# Generic methods # +################### + +(A::Tridiagonal, B::Tridiagonal) = Tridiagonal(A.dl+B.dl, A.d+B.d, A.du+B.du) -(A::Tridiagonal, B::Tridiagonal) = Tridiagonal(A.dl-B.dl, A.d-B.d, A.du+B.du) *(A::Tridiagonal, B::Number) = Tridiagonal(A.dl*B, A.d*B, A.du*B) @@ -185,6 +248,54 @@ end ==(A::Tridiagonal, B::SymTridiagonal) = (A.dl==A.du==B.ev) && (A.d==B.dv) ==(A::SymTridiagonal, B::SymTridiagonal) = B==A +#Implements the inverse using principal minors +#Ref: +# R. Usmani, "Inversion of a tridiagonal Jacobi matrix", +# Linear Algebra and its Applications 212-213 (1994), pp.413-414 +# doi:10.1016/0024-3795(94)90414-6 +function inv{T}(A::Tridiagonal{T}) + n = size(A, 1) + θ = ZeroOffsetVector(zeros(T, n+1)) #principal minors of A + θ[0] = 1 + n>=1 && (θ[1] = A.d[1]) + for i=2:n + θ[i] = A.d[i]*θ[i-1]-A.dl[i-1]*A.du[i-1]*θ[i-2] + end + φ = zeros(T, n+1) + φ[n+1] = 1 + n>=1 && (φ[n] = A.d[n]) + for i=n-1:-1:1 + φ[i] = A.d[i]*φ[i+1]-A.du[i]*A.dl[i]*φ[i+2] + end + α = Array(T, n, n) + for i=1:n, j=1:n + sign = (i+j)%2==0 ? (+) : (-) + if ij + α[i,j]=(sign)(prod(A.dl[j:i-1]))*θ[j-1]*φ[i+1]/θ[n] + end + end + α +end + +#Implements the determinant using principal minors +# R. Usmani, "Inversion of a tridiagonal Jacobi matrix", +# Linear Algebra and its Applications 212-213 (1994), pp.413-414 +# doi:10.1016/0024-3795(94)90414-6 +function det{T}(A::Tridiagonal{T}) + n = size(A, 1) + θa = one(T) + n==0 && return θa + θb = A.d[1] + for i=2:n + θb, θa = A.d[i]*θb-A.dl[i-1]*A.du[i-1]*θa, θb + end + return θb +end + # Elementary operations that mix Tridiagonal and SymTridiagonal matrices convert(::Type{Tridiagonal}, A::SymTridiagonal) = Tridiagonal(A.ev, A.dv, A.ev) +(A::Tridiagonal, B::SymTridiagonal) = Tridiagonal(A.dl+B.ev, A.d+B.dv, A.du+B.ev) diff --git a/test/linalg.jl b/test/linalg.jl index 466c531b4bdd5..be60400354dde 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -603,7 +603,41 @@ for relty in (Float16, Float32, Float64, BigFloat), elty in (relty, Complex{relt end end +#Tridiagonal matrices +for relty in (Float16, Float32, Float64), elty in (relty, Complex{relty}) + a = convert(Vector{elty}, randn(n-1)) + b = convert(Vector{elty}, randn(n)) + c = convert(Vector{elty}, randn(n-1)) + if elty <: Complex + a += im*convert(Vector{elty}, randn(n-1)) + b += im*convert(Vector{elty}, randn(n)) + c += im*convert(Vector{elty}, randn(n-1)) + end + + A=Tridiagonal(a, b, c) + fA=(elty<:Complex?complex128:float64)(full(A)) + for func in (det, inv) + @test_approx_eq_eps func(A) func(fA) n^2*sqrt(eps(relty)) + end +end + #SymTridiagonal (symmetric tridiagonal) matrices +for relty in (Float16, Float32, Float64), elty in (relty, )#XXX Complex{relty}) doesn't work + a = convert(Vector{elty}, randn(n)) + b = convert(Vector{elty}, randn(n-1)) + if elty <: Complex + relty==Float16 && continue + a += im*convert(Vector{elty}, randn(n)) + b += im*convert(Vector{elty}, randn(n-1)) + end + + A=SymTridiagonal(a, b) + fA=(elty<:Complex?complex128:float64)(full(A)) + for func in (det, inv) + @test_approx_eq_eps func(A) func(fA) n^2*sqrt(eps(relty)) + end +end + Ainit = randn(n) Binit = randn(n-1) for elty in (Float32, Float64)