diff --git a/base/export.jl b/base/export.jl index 8e2279ddf05cd..060ef03779ade 100644 --- a/base/export.jl +++ b/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, diff --git a/base/linalg.jl b/base/linalg.jl index 771e4b3263495..8ee32e2b0e0ae 100644 --- a/base/linalg.jl +++ b/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 diff --git a/base/linalg_lapack.jl b/base/linalg_lapack.jl index 3c4c49a10ab0d..bd24799b655d5 100644 --- a/base/linalg_lapack.jl +++ b/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)