diff --git a/base/linalg/lapack.jl b/base/linalg/lapack.jl index d958d5204ec72..3637709760792 100644 --- a/base/linalg/lapack.jl +++ b/base/linalg/lapack.jl @@ -1834,20 +1834,13 @@ for (trtri, trtrs, trevc, elty) 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) @@ -1855,27 +1848,27 @@ for (trtri, trtrs, trevc, elty) 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 d94f72958f326..5ac0485e45549 100644 --- a/base/linalg/triangular.jl +++ b/base/linalg/triangular.jl @@ -16,13 +16,11 @@ function Triangular(A::Matrix) end size(A::Triangular, args...) = size(A.UL, args...) -function full(A::Triangular) - if - istril(A) return tril(A.UL) - else - return triu(A.UL) - end -end + +#Converting from Triangular to dense Matrix +convert(::Type{Matrix}, A::Triangular) = full(A) +full(A::Triangular) = istril(A) ? tril(A.UL) : triu(A.UL) + print_matrix(io::IO, A::Triangular, rows::Integer, cols::Integer) = print_matrix(io, full(A), rows, cols) istril(A::Triangular) = A.uplo == 'L' @@ -68,3 +66,41 @@ function inv{T<:BlasFloat}(A::Triangular{T}) return Ainv end inv(A::Triangular) = inv(Triangular(float(A.UL), A.uplo, A.unitdiag)) +diag(A::Triangular) = diag(A.UL) +getindex(A::Triangular,m::Int,n::Int) = getindex(A.UL, m, n) + +####################### +# 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 7230a396532f7..a484797eb1c8c 100644 --- a/test/linalg.jl +++ b/test/linalg.jl @@ -489,7 +489,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) @@ -524,7 +546,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) @@ -547,24 +569,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