Fix of svd to allow for zero dimensional matrices. sdd has been remov…

…ed and svd calls gesdd. Thin svd method has been added. Change eig to allow for vals only method for efficiency. Added eigvals function for convenience. Updated exported names for the changes in svd and eig. Change norm and rank to allow for zero dim matrices.

base/export.jl

@@ -634,6 +634,7 @@ export
     diagmm!,
     dot,
     eig,
+    eigvals,
     expm,
     eye,
     factors,
@@ -661,7 +662,6 @@ export
     randsym,
     rank,
     rref,
-    sdd,
     solve,
     svd,
     svdvals,

base/linalg.jl

@@ -314,7 +314,9 @@ diag(A::AbstractVector) = error("Perhaps you meant to use diagm().")
 #diagm{T}(v::Union(AbstractVector{T},AbstractMatrix{T}))
 
 function norm(x::AbstractVector, p::Number)
-    if p == Inf
+    if length(x) == 0
+        return 0.0
+    elseif p == Inf
         return max(abs(x))
     elseif p == -Inf
         return min(abs(x))
@@ -326,14 +328,17 @@ end
 norm(x::AbstractVector) = sqrt(real(dot(x,x)))
 
 function norm(A::AbstractMatrix, p)
-    if size(A,1) == 1 || size(A,2) == 1
+    m, n = size(A)
+    if m == 0 || n == 0
+        return 0.0
+    elseif m == 1 || n == 1
         return norm(reshape(A, numel(A)), p)
     elseif p == 1
         return max(sum(abs(A),1))
     elseif p == 2
         return max(svd(A)[2])
     elseif p == Inf
-        max(sum(abs(A),2))
+        return max(sum(abs(A),2))
     elseif p == "fro"
         return sqrt(sum(diag(A'*A)))
     else
@@ -344,6 +349,8 @@ end
 norm(A::AbstractMatrix) = norm(A, 2)
 rank(A::AbstractMatrix, tol::Real) = sum(svdvals(A) .> tol)
 function rank(A::AbstractMatrix)
+    m,n = size(A)
+    if m == 0 || n == 0; return 0; end
     sv = svdvals(A)
     sum(sv .> max(size(A,1),size(A,2))*eps(sv[1]))
 end

base/linalg_lapack.jl

@@ -1,41 +1,58 @@
 # linear algebra functions that use the Lapack module
 
-eig{T<:Integer}(x::StridedMatrix{T}) = eig(float64(x))
-
-function eig{T<:LapackScalar}(A::StridedMatrix{T})
-    if ishermitian(A) return Lapack.syev!('V','U',copy(A)) end
-                                        # Only compute right eigenvectors
-    if iscomplex(A) return Lapack.geev!('N','V',copy(A))[2:3] end
-    VL, WR, WI, VR = Lapack.geev!('N','V',copy(A))
-    if all(WI .== 0.) return WR, VR end
+function eig{T<:LapackScalar}(A::StridedMatrix{T}, vecs::Bool)
     n = size(A, 2)
-    evec = complex(zeros(T, n, n))
-    j = 1
-    while j <= n
-        if WI[j] == 0.0
-            evec[:,j] = VR[:,j]
-        else
-            evec[:,j]   = VR[:,j] + im*VR[:,j+1]
-            evec[:,j+1] = VR[:,j] - im*VR[:,j+1]
+    if n == 0; return vecs ? (zeros(T, 0), zeros(T, 0, 0)) : zeros(T, 0, 0); end
+
+    if ishermitian(A); return Lapack.syev!(vecs ? 'V' : 'N', 'U', copy(A)); end
+
+    if iscomplex(A)
+        W, VR = Lapack.geev!('N', vecs ? 'V' : 'N', copy(A))[2:3]
+        if vecs; return W, VR; end
+        return W
+    end
+
+    VL, WR, WI, VR = Lapack.geev!('N', vecs ? 'V' : 'N', copy(A))
+    if all(WI .== 0.) 
+        if vecs; return WR, VR; end
+        return WR
+    end
+    if vecs    
+        evec = complex(zeros(T, n, n))
+        j = 1
+        while j <= n
+            if WI[j] == 0.0
+                evec[:,j] = VR[:,j]
+            else
+                evec[:,j]   = VR[:,j] + im*VR[:,j+1]
+                evec[:,j+1] = VR[:,j] - im*VR[:,j+1]
+                j += 1
+            end
             j += 1
         end
-        j += 1
+        return complex(WR, WI), evec
     end
-    complex(WR, WI), evec
+    complex(WR, WI)
 end
 
-sdd!{T<:LapackScalar}(A::StridedMatrix{T},vecs::Char) = Lapack.gesdd!(vecs, copy(A))
-sdd{T<:LapackScalar}(A::StridedMatrix{T},vecs::Char) = sdd!(copy(A), vecs)
-sdd{T<:Real}(x::StridedMatrix{T},vecs::Char) = sdd(float64(x),vecs)
-sdd(A) = sdd(A, 'A')
+eig{T<:Integer}(x::StridedMatrix{T}, vecs::Bool) = eig(float64(x), vecs)
+eig(x::StridedMatrix) = eig(x, true)
+eigvals(x::StridedMatrix) = eig(x, false)
 
-function svd{T<:LapackScalar}(A::StridedMatrix{T},vecs::Bool)
-    Lapack.gesvd!(vecs ? 'A' : 'N', vecs ? 'A' : 'N', copy(A))
+function svd{T<:LapackScalar}(A::StridedMatrix{T},vecs::Bool,thin::Bool)
+    m,n = size(A)
+    if m == 0 || n == 0
+        if vecs; return (eye(m, thin ? n : m), zeros(0), eye(n,n)); end
+        return (zeros(T, 0, 0), zeros(T, 0), zeros(T, 0, 0))
+    end
+    if vecs; return Lapack.gesdd!(thin ? 'O' : 'A', copy(A)); end
+    Lapack.gesdd!('N', copy(A))
 end
 
-svd{T<:Integer}(x::StridedMatrix{T},vecs) = svd(float64(x),vecs)
-svd(A) = svd(A,true)
-svdvals(A) = svd(A,false)[2]
+svd{T<:Integer}(x::StridedMatrix{T},vecs,thin) = svd(float64(x),vecs,thin)
+svd(A::StridedMatrix) = svd(A,true,false)
+svd(A::StridedMatrix, thin::Bool) = svd(A,true,thin)
+svdvals(A) = svd(A,false,true)[2]
 
 function (\){T<:LapackScalar}(A::StridedMatrix{T}, B::StridedVecOrMat{T})
     Acopy = copy(A)