Skip to content

Commit

Permalink
Provides diag, eig, eigfact eigvals and eigvecs for upper and lower T…
Browse files Browse the repository at this point in the history 
…riangular 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
  • Loading branch information
@jiahao
jiahao committed Dec 29, 2013
1 parent 10583b8 commit 20cc7e9
Show file tree
Hide file tree
Showing 3 changed files with 83 additions and 30 deletions.
27 changes: 10 additions & 17 deletions base/linalg/lapack.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -1899,48 +1899,41 @@ 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)
ccall(($(string(trevc)),liblapack), Void,
(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
Expand Down
44 changes: 42 additions & 2 deletions base/linalg/triangular.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)]...)

Expand Down Expand Up @@ -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))
42 changes: 31 additions & 11 deletions test/linalg.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -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)
Expand Down Expand Up @@ -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)
Expand All @@ -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
Expand Down

0 comments on commit 20cc7e9

Please sign in to comment.