Provides diag, eig, eigfact eigvals and eigvecs for upper and lower T…

…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

base/linalg/lapack.jl

@@ -1834,48 +1834,41 @@ 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)
           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

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))
+

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