Skip to content

Commit

Permalink
Adds eig, eigfact, eigvals and eigvecs for Bidiagonal
Browse files Browse the repository at this point in the history
Addresses #3688
  • Loading branch information
jiahao committed Jan 7, 2014
1 parent 4fec9ba commit 5ce520c
Showing 1 changed file with 38 additions and 0 deletions.
38 changes: 38 additions & 0 deletions base/linalg/bidiag.jl
Original file line number Diff line number Diff line change
Expand Up @@ -102,6 +102,44 @@ function \{T<:BlasFloat}(M::Bidiagonal{T}, rhs::StridedVecOrMat{T})
apply(LAPACK.gtsv!, M.isupper ? (z, copy(M.dv), copy(M.ev), copy(rhs)) : (copy(M.ev), copy(M.dv), z, copy(rhs)))
end


#######################
# Eigenvalues/vectors #
#######################

eigvals{T<:Number}(M::Bidiagonal{T}) = M.isupper ? M.dv : reverse(M.dv)
function eigvecs{T<:Number}(M::Bidiagonal{T})
n = length(M.dv)
Q=zeros(T, n, n)
v=zeros(T, n)
if M.isupper
for i=1:n #index of eigenvector
v[1] = convert(T, 1.0)
for j=1:i-1 #Starting from j=i, eigenvector elements will be 0
v[j+1] = (M.dv[i]-M.dv[j])/M.ev[j] * v[j]
end
v /= norm(v)
Q[:,i] = v
end
else
for i=n:-1:1 #index of eigenvector
v[n] = convert(T, 1.0)
for j=(n-1):-1:max(1,(i-1)) #Starting from j=i-1, eigenvector elements will be 0
v[j] = (M.dv[i]-M.dv[j+1])/M.ev[j] * v[j+1]
end
v /= norm(v)
Q[:,n+1-i] = v
end
end
Q
end

eig{T<:Number}(M::Bidiagonal{T}) = eigvals(M), eigvecs(M)
eigfact{T<:Number}(M::Bidiagonal{T}) = Eigen{T,T}(eigvals(M), eigvecs(M))

###################
# Singular values #
###################
#Wrap bdsdc to compute singular values and vectors
svdvals{T<:Real}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'N', copy(M.dv), copy(M.ev))
svd {T<:Real}(M::Bidiagonal{T})=LAPACK.bdsdc!(M.isupper?'U':'L', 'I', copy(M.dv), copy(M.ev))
Expand Down

0 comments on commit 5ce520c

Please sign in to comment.