Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Various changes to svd, eig, norm and rank #1312

Merged
merged 1 commit into from
Sep 29, 2012
Merged
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion base/export.jl
Original file line number Diff line number Diff line change
Expand Up @@ -634,6 +634,7 @@ export
diagmm!,
dot,
eig,
eigvals,
expm,
eye,
factors,
Expand Down Expand Up @@ -661,7 +662,6 @@ export
randsym,
rank,
rref,
sdd,
solve,
svd,
svdvals,
Expand Down
13 changes: 10 additions & 3 deletions base/linalg.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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))
Expand All @@ -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
Expand All @@ -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
Expand Down
71 changes: 44 additions & 27 deletions base/linalg_lapack.jl
Original file line number Diff line number Diff line change
@@ -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)
Expand Down