From 20cc7e9fb2cc3aa64171ec4d8f1ab620f033e52a Mon Sep 17 00:00:00 2001 From: Jiahao Chen Date: Sat, 19 Oct 2013 01:30:47 -0400 Subject: [PATCH] Provides diag, eig, eigfact eigvals and eigvecs for upper and lower Triangular matrices Addresses #3688 - Also provides simple wrappers around the full routines for svd, svdfact, svdfact!, svdvals, svdvecs for Triangular - Adds convert routine from Triangular to dense Matrix - Clean up bidiagonal Matrix tests --- base/linalg/lapack.jl | 27 +++++++++--------------- base/linalg/triangular.jl | 44 +++++++++++++++++++++++++++++++++++++-- test/linalg.jl | 42 +++++++++++++++++++++++++++---------- 3 files changed, 83 insertions(+), 30 deletions(-) diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index 87a2e9f1bb85e..5a898ccd7dc3e 100644 --- a/base/linalg/lapack.jl +++ b/base/linalg/lapack.jl @@ -1899,20 +1899,13 @@ for (trcon, elty, relty) in # DOUBLE PRECISION T( LDT, * ), VL( LDVL, * ), VR( LDVR, * ), #$ WORK( * ) function trevc!(side::BlasChar, howmny::BlasChar, - select::StridedMatrix{Bool}, A::StridedMatrix{$elty}, + select::Vector{Bool}, A::StridedMatrix{$elty}, VL::StridedMatrix{$elty}, VR::StridedMatrix{$elty}) chkstride1(A) chksquare(A) - if !(side in ['R', 'L', 'B']) error("Unsupported value of side") end - if !(howmny in ['A', 'B', 'S']) error("Unsupported value of howmny") end ldt, n = size(A) - if n!=length(select) error("Wrong size of select array") end - if ldt < max(1,n) error("Wrong dimension 1 of A") end ldvl, mm = size(VL) - if ldvl < side=='A' ? 1 : n error("Wrong dimension 1 of VL") end - ldvr, mm2= size(VR) - if mm!=mm2 error("VL and VR have different dimension 2s") end - if ldvr < side=='A' ? 1 : n error("Wrong dimension 1 of VR") end + ldvr, mm = size(VR) m = Array(BlasInt, 1) work = Array($elty, 3n) info = Array(BlasInt, 1) @@ -1920,27 +1913,27 @@ for (trcon, elty, relty) in (Ptr{BlasChar}, Ptr{BlasChar}, Ptr{Bool}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}), - &side, &howmny, &select, &n, &A, &ldt, &VL, &ldvl, &VR, &ldvr, &mm, - &m, &work, &info + &side, &howmny, select, &n, A, &ldt, VL, &ldvl, VR, &ldvr, &mm, + m, work, info ) if info[1] < 0 throw(LAPACKException(info[1])) end #Decide what exactly to return if howmny=='S' #compute selected eigenvectors if side=='L' #left eigenvectors only - return select, VL[:,1:m] + return select, VL[:,1:m[1]] elseif side=='R' #right eigenvectors only - return select, VR[:,1:m] + return select, VR[:,1:m[1]] else #side=='B' #both eigenvectors - return select, VL[:,1:m], VR[:,1:m] + return select, VL[:,1:m[1]], VR[:,1:m[1]] end else #compute all eigenvectors if side=='L' #left eigenvectors only - return VL[:,1:m] + return VL[:,1:m[1]] elseif side=='R' #right eigenvectors only - return VR[:,1:m] + return VR[:,1:m[1]] else #side=='B' #both eigenvectors - return VL[:,1:m], VR[:,1:m] + return VL[:,1:m[1]], VR[:,1:m[1]] end end end diff --git a/base/linalg/triangular.jl b/base/linalg/triangular.jl index 2f25b31188856..fa0f749d3c281 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -57,16 +57,20 @@ end # Generic routines # #################### +size(A::Triangular, args...) = size(A.UL, args...) + +print_matrix(io::IO, A::Triangular, rows::Integer, cols::Integer) = print_matrix(io, full(A), rows, cols) + +convert(::Type{Matrix}, A::Triangular) = full(A) full(A::Triangular) = (istril(A) ? tril! : triu!)(A.UL) getindex{T}(A::Triangular{T}, i::Integer, j::Integer) = i == j ? A.UL[i,j] : ((A.uplo == 'U') == (i < j) ? getindex(A.UL, i, j) : zero(T)) istril(A::Triangular) = A.uplo == 'L' || istriu(A.UL) istriu(A::Triangular) = A.uplo == 'U' || istril(A.UL) -size(A::Triangular, args...) = size(A.UL, args...) - transpose(A::Triangular) = Triangular(A.UL, A.uplo=='U':'L':'U', A.unitdiag) ctranspose(A::Triangular) = conj(transpose(A)) +diag(A::Triangular) = diag(A.UL) #Generic solver using naive substitution function naivesub!(A::Triangular, b::AbstractVector, x::AbstractVector=b) @@ -98,6 +102,7 @@ function naivesub!(A::Triangular, b::AbstractVector, x::AbstractVector=b) end x end + \{T<:Number}(A::Triangular{T}, b::AbstractVector{T}) = naivesub!(A, b, similar(b)) \{T<:Number}(A::Triangular{T}, B::AbstractMatrix{T}) = hcat([naivesub!(A, B[:,i], similar(B[:,i])) for i=1:size(B,2)]...) @@ -148,3 +153,38 @@ function eigvecs{T}(A::Triangular{T}) end eigfact(A::Triangular) = Eigen(eigvals(A), eigvecs(A)) + +####################### +# Eigenvalues/vectors # +####################### + +eigvals(A::Triangular) = A.uplo=='U' ? diag(A) : reverse(diag(A)) +function eigvecs{T<:BlasFloat}(A::Triangular{T}) + V = LAPACK.trevc!('R', 'A', Array(Bool,1), A.uplo=='U' ? A.UL : transpose(A.UL), + Array(T,size(A)), Array(T, size(A))) + if A.uplo=='L' #This is the transpose of the Schur form + #The eigenvectors must be transformed + VV = inv(Triangular(transpose(V))) + N = size(V,2) + for i=1:N #Reorder eigenvectors to follow LAPACK convention + V[:,i]=VV[:,N+1-i] + end + end + #Need to normalize + for i=1:size(V,2) + V[:,i] /= norm(V[:,i]) + end + V +end +eig(M::Triangular) = eigvals(M), eigvecs(M) +eigfact(M::Triangular) = Eigen(eigvals(M), eigvecs(M)) + +############################# +# Singular values / vectors # +############################# + +svd(M::Triangular) = svd(full(M)) +svdfact(M::Triangular) = svdfact(full(M)) +svdfact!(M::Triangular) = svdfact!(full(M)) +svdvals(M::Triangular) = svdvals(full(M)) +svdvecs(M::Triangular) = svdvecs(full(M)) diff --git a/test/linalg.jl b/test/linalg.jl index 36f46fe76bc1e..195843f753cc3 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -601,7 +601,29 @@ function test_approx_eq_vecs(a, b) end end -#LAPACK tests for symmetric tridiagonal matrices +############################## +# Tests for special matrices # +############################## + +#Triangular matrices +n=5 +for elty in (Float32, Float64) + AUfull = convert(Matrix{elty}, sum([diagm(randn(n-i), i) for i=0:n-1])) + for Afull in {AUfull, AUfull'} #Test upper and lower triangular + A = Triangular(Afull) + #Test eigenvalues/vectors against dense routines + d1, d2 = eigvals(A), eigvals(Afull) + v1, v2 = eigvecs(A), eigvecs(Afull) + @test_approx_eq d1 d2 + test_approx_eq_vecs(v1, v2) + #Test spectral decomposition + #XXX This is not the correct error bound + @test_approx_eq_eps 0 norm(v1 * diagm(d1) * inv(v1) - Afull) sqrt(eps(elty)) + @test_approx_eq_eps 0 norm(v2 * diagm(d2) * inv(v2) - Afull) sqrt(eps(elty)) + end +end + +#SymTridiagonal (symmetric tridiagonal) matrices n=5 Ainit = randn(n) Binit = randn(n-1) @@ -630,7 +652,7 @@ for elty in (Float32, Float64) end -#Test bidiagonal matrices and their SVDs +#Bidiagonal matrices dv = randn(n) ev = randn(n-1) for elty in (Float32, Float64, Complex64, Complex128) @@ -653,24 +675,22 @@ for elty in (Float32, Float64, Complex64, Complex128) @test ctranspose(ctranspose(T)) == T if (elty <: Real) - #XXX If I run either of these tests separately, by themselves, things are OK. - # Enabling BOTH tests results in segfault. - # Where is the memory corruption??? - - @test_approx_eq svdvals(full(T)) svdvals(T) - u1, d1, v1 = svd(full(T)) + Tfull = full(T) + #Test singular values/vectors + @test_approx_eq svdvals(Tfull) svdvals(T) + u1, d1, v1 = svd(Tfull) u2, d2, v2 = svd(T) @test_approx_eq d1 d2 test_approx_eq_vecs(u1, u2) test_approx_eq_vecs(v1, v2) #Test eigenvalues/vectors - d1, v1 = eig(full(T)) + d1, v1 = eig(Tfull) d2, v2 = eigvals(T), eigvecs(T) @test_approx_eq d1 d2 test_approx_eq_vecs(v1, v2) - @test_approx_eq_eps 0 norm(v1 * diagm(d1) * inv(v1) - full(T)) 1e-14 - @test_approx_eq_eps 0 norm(v2 * diagm(d2) * inv(v2) - full(T)) 1e-14 + @test_approx_eq_eps 0 norm(v1 * diagm(d1) * inv(v1) - Tfull) eps(elty)*n*(n+1) + @test_approx_eq_eps 0 norm(v2 * diagm(d2) * inv(v2) - Tfull) eps(elty)*n*(n+1) end end end