Return vector instead of matrix for eigenvalues and singular values.

Close issue #177

j/linalg_arpack.j

@@ -230,7 +230,7 @@ function eigs{T}(A::AbstractMatrix{T}, k::Int, evtype::ASCIIString)
     if (info[1] != 0); error("Error in ARPACK eupd"); end
 
     if iscomplex(A) || isrealsymA
-        return (diagm(d), v[1:n, 1:nev])
+        return (d, v[1:n, 1:nev])
     else
         evec = complex(zeros(T, n, nev+1), zeros(T, n, nev+1))
         for j=1:nev
@@ -242,7 +242,7 @@ function eigs{T}(A::AbstractMatrix{T}, k::Int, evtype::ASCIIString)
                 j += 1
             end
         end
-        return (diagm(complex(dr[1:nev],di[1:nev])), evec[1:n, 1:nev])
+        return (complex(dr[1:nev],di[1:nev]), evec[1:n, 1:nev])
     end
 
 end
@@ -307,6 +307,6 @@ function svds{T}(A::AbstractMatrix{T}, k::Int)
     v = v[1:n, 1:nev]
     u = A*v*diagm(1./d)
 
-    return (u, diagm(d), v.')
+    return (u, d, v.')
 
 end

j/linalg_lapack.j

@@ -363,7 +363,7 @@ function eig{T}(A::Matrix{T})
             info = jl_lapack_syev(jobz, uplo, n, EV, stride(A,2), W, work, lwork)
         end
 
-        if info == 0; return (diagm(W), EV); end
+        if info == 0; return (W, EV); end
         error("Error in LAPACK syev/heev");
 
     else # Non-symmetric case
@@ -404,7 +404,7 @@ function eig{T}(A::Matrix{T})
         if info != 0; error("Error in LAPACK geev"); end
 
         if iscomplex(A)
-            return (diagm(W), VR)
+            return (W, VR)
         else
             evec = complex(zeros(T, n, n), zeros(T, n, n))
             for j=1:n
@@ -416,7 +416,7 @@ function eig{T}(A::Matrix{T})
                     j += 1
                 end
             end
-            return (diagm(complex(WR, WI)), evec)
+            return (complex(WR, WI), evec)
         end
 
     end # symmetric / non-symmetric case
@@ -506,10 +506,7 @@ function svd{T}(A::Matrix{T})
         info = jl_lapack_gesvd(jobu, jobvt, m, n, X, stride(A,2), S, U, m, VT, n, work, lwork)
     end
 
-    SIGMA = zeros(T, m, n)
-    for i=1:k; SIGMA[i,i] = S[i]; end
-
-    if info == 0; return (U, SIGMA, VT); end
+    if info == 0; return (U, S, VT); end
     error("Error in LAPACK gesvd");
 end
 

test/tests.j

@@ -766,8 +766,10 @@ r = chol(asym)
 @assert sum(q*r[:,p] - a) < 1e-8
 (d,v) = eig(asym)
 @assert sum(asym*v[:,1]-d[1]*v[:,1]) < 1e-8
+(d,v) = eig(a)
+@assert abs(sum(a*v[:,1]-d[1]*v[:,1])) < 1e-8
 (u,s,vt) = svd(a)
-@assert sum(u*s*vt - a) < 1e-8
+@assert sum(u*diagm(s)*vt - a) < 1e-8
 x = a \ b
 @assert sum(a*x-b) < 1e-8
 x = triu(a) \ b
@@ -777,10 +779,10 @@ x = tril(a) \ b
 
 # arpack
 (d,v) = eigs(asym, 3)
-@assert sum(asym*v[:,1]-d[1,1]*v[:,1]) < 1e-8
+@assert sum(asym*v[:,1]-d[1]*v[:,1]) < 1e-8
 
 (d,v) = eigs(a,3)
-@assert abs(sum(a*v[:,2]-d[2,2]*v[:,2])) < 1e-8
+@assert abs(sum(a*v[:,2]-d[2]*v[:,2])) < 1e-8
 
 # hash table
 h = HashTable()

test/unittests.j

@@ -200,7 +200,7 @@ a = rand(n,n)
 if WORD_SIZE==64
     # TODO: hangs on 32-bit
     (d,v) = eigs(a, 3)
-    @assert abs(sum(a*v[:,1]-d[1,1]*v[:,1])) < 1e-8
+    @assert abs(sum(a*v[:,1]-d[1]*v[:,1])) < 1e-8
 end
 
 # hash table