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.

Sign up for GitHub

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

Jump to bottom

Add arguments specifying subsets of eigenvalues in eigfact(Symmetric/Hermitian) #6652

Merged
merged 1 commit into from
Apr 27, 2014
Merged

Add arguments specifying subsets of eigenvalues in eigfact(Symmetric/Hermitian) #6652

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
1 commit
Select commit
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
9 changes: 2 additions & 7 deletions base/linalg/factorization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -548,8 +548,8 @@ eigfact(x::Number) = Eigen([x], fill(one(x), 1, 1))
# F = eigfact(A, permute=permute, scale=scale)
# F[:values], F[:vectors]
# end
function eig(A::Union(Number, AbstractMatrix); kwargs...)
F = eigfact(A, kwargs...)
function eig(A::Union(Number, AbstractMatrix), args...; kwargs...)
F = eigfact(A, args..., kwargs...)
F[:values], F[:vectors]
end
#Calculates eigenvectors
Expand Down Expand Up @@ -612,11 +612,6 @@ end
eigfact{T<:BlasFloat}(A::AbstractMatrix{T}, B::StridedMatrix{T}) = eigfact!(copy(A),copy(B))
eigfact{TA,TB}(A::AbstractMatrix{TA}, B::AbstractMatrix{TB}) = (S = promote_type(Float32,typeof(one(TA)/norm(one(TA))),TB); eigfact!(S != TA ? convert(Matrix{S},A) : copy(A), S != TB ? convert(Matrix{S},B) : copy(B)))

function eig(A::AbstractMatrix, B::AbstractMatrix)
F = eigfact(A, B)
F[:values], F[:vectors]
end

function eigvals!{T<:BlasReal}(A::StridedMatrix{T}, B::StridedMatrix{T})
issym(A) && issym(B) && return eigvals!(Symmetric(A), Symmetric(B))
alphar, alphai, beta, vl, vr = LAPACK.ggev!('N', 'N', A, B)
Expand Down
13 changes: 6 additions & 7 deletions base/linalg/symmetric.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -37,14 +37,13 @@ ctranspose(A::Hermitian) = A
factorize(A::HermOrSym) = bkfact(A.S, symbol(A.uplo), issym(A))
\(A::HermOrSym, B::StridedVecOrMat) = \(bkfact(A.S, symbol(A.uplo), issym(A)), B)

eigfact!{T<:BlasReal}(A::Symmetric{T}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.S, 0.0, 0.0, 0, 0, -1.0)...)
eigfact!{T<:BlasComplex}(A::Hermitian{T}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.S, 0.0, 0.0, 0, 0, -1.0)...)
eigfact{T<:BlasFloat}(A::HermOrSym{T}) = eigfact!(copy(A))
eigfact!{T<:BlasFloat}(A::HermOrSym{T}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.S, 0.0, 0.0, 0, 0, -1.0)...)
eigfact!{T<:BlasFloat}(A::HermOrSym{T}, il::Int, ih::Int) = Eigen(LAPACK.syevr!('V', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)...)
eigfact!{T<:BlasFloat}(A::HermOrSym{T}, vl::Real, vh::Real) = Eigen(LAPACK.syevr!('V', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)...)
eigfact{T}(A::HermOrSym{T}) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? eigfact!(convert(AbstractMatrix{S}, A)) : eigfact!(copy(A)))
eigvals!{T<:BlasReal}(A::Symmetric{T}, il::Int=1, ih::Int=size(A,1)) = LAPACK.syevr!('N', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)[1]
eigvals!{T<:BlasReal}(A::Symmetric{T}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)[1]
eigvals!{T<:BlasComplex}(A::Hermitian{T}, il::Int=1, ih::Int=size(A,1)) = LAPACK.syevr!('N', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)[1]
eigvals!{T<:BlasComplex}(A::Hermitian{T}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)[1]
eigfact{T,U<:Number}(A::HermOrSym{T},l::U,h::U) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? eigfact!(convert(AbstractMatrix{S}, A),l,h) : eigfact!(copy(A),l,h))
eigvals!{T<:BlasFloat}(A::HermOrSym{T}, il::Int=1, ih::Int=size(A,1)) = LAPACK.syevr!('N', 'I', A.uplo, A.S, 0.0, 0.0, il, ih, -1.0)[1]
eigvals!{T<:BlasFloat}(A::HermOrSym{T}, vl::Real, vh::Real) = LAPACK.syevr!('N', 'V', A.uplo, A.S, vl, vh, 0, 0, -1.0)[1]
eigvals{T<:BlasFloat}(A::HermOrSym{T},l::Real=1,h::Real=size(A,1)) = eigvals!(copy(A),l,h)
eigvals{T}(A::HermOrSym{T},l::Real=1,h::Real=size(A,1)) = (S = promote_type(Float32,typeof(one(T)/norm(one(T)))); S != T ? eigvals!(convert(AbstractMatrix{S}, A, l, h)) : eigvals!(copy(A), l, h))
eigmax(A::HermOrSym) = eigvals(A, size(A, 1), size(A, 1))[1]
Expand Down
12 changes: 6 additions & 6 deletions doc/stdlib/linalg.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -179,17 +179,17 @@ Linear algebra functions in Julia are largely implemented by calling functions f

``sqrtm`` uses a polyalgorithm, computing the matrix square root using Schur factorizations (:func:`schurfact`) unless it detects the matrix to be Hermitian or real symmetric, in which case it computes the matrix square root from an eigendecomposition (:func:`eigfact`). In the latter situation for positive definite matrices, the matrix square root has ``Real`` elements, otherwise it has ``Complex`` elements.

.. function:: eig(A,[permute=true,][scale=true]) -> D, V
.. function:: eig(A,[il,][iu,][vl,][vu,][permute=true,][scale=true]) -> D, V

Wrapper around ``eigfact`` extracting all parts the factorization to a tuple. Direct use of ``eigfact`` is therefore generally more efficient. Computes eigenvalues and eigenvectors of ``A``. See :func:`eigfact` for details on the ``permute`` and ``scale`` keyword arguments.
Wrapper around ``eigfact`` extracting all parts the factorization to a tuple. Direct use of ``eigfact`` is therefore generally more efficient. Computes eigenvalues and eigenvectors of ``A``. See :func:`eigfact` for details on adiitional arguments.

.. function:: eig(A, B) -> D, V

Wrapper around ``eigfact`` extracting all parts the factorization to a tuple. Direct use of ``eigfact`` is therefore generally more efficient. Computes generalized eigenvalues and vectors of ``A`` with respect to ``B``.

.. function:: eigvals(A)
.. function:: eigvals(A,[il,][iu,][vl,][vu])

Returns the eigenvalues of ``A``.
Returns the eigenvalues of ``A``. If ``A`` is ``Symmetric`` or ``Hermitian``, it is possible to calculate only a subset of the eigenvalues by specifying either a pair ``il`` and `iu`` for the lower and upper indices of the sorted eigenvalues or a pair ``vl`` and ``vu`` for the lower and upper boundaries of the eigenvalues.

.. function:: eigmax(A)

Expand All @@ -206,11 +206,11 @@ Linear algebra functions in Julia are largely implemented by calling functions f

For ``SymTridiagonal`` matrices, if the optional vector of eigenvalues ``eigvals`` is specified, returns the specific corresponding eigenvectors.

.. function:: eigfact(A,[permute=true,][scale=true])
.. function:: eigfact(A,[il,][iu,][vl,][vu,][permute=true,][scale=true])

Compute the eigenvalue decomposition of ``A`` and return an ``Eigen`` object. If ``F`` is the factorization object, the eigenvalues can be accessed with ``F[:values]`` and the eigenvectors with ``F[:vectors]``. The following functions are available for ``Eigen`` objects: ``inv``, ``det``.

For general non-symmetric matrices it is possible to specify how the matrix is balanced before the eigenvector calculation. The option ``permute=true`` permutes the matrix to become closer to upper triangular, and ``scale=true`` scales the matrix by its diagonal elements to make rows and columns more equal in norm. The default is ``true`` for both options.
If ``A`` is ``Symmetric`` or ``Hermitian`` it is possible to calculate only a subset of the eigenvalues by specifying either a pair ``il`` and `iu`` for the lower and upper indices of the sorted eigenvalues or a pair ``vl`` and ``vu`` for the lower and upper boundaries of the eigenvalues. For general non-symmetric matrices it is possible to specify how the matrix is balanced before the eigenvector calculation. The option ``permute=true`` permutes the matrix to become closer to upper triangular, and ``scale=true`` scales the matrix by its diagonal elements to make rows and columns more equal in norm. The default is ``true`` for both options.

.. function:: eigfact(A, B)

Expand Down
7 changes: 4 additions & 3 deletions test/linalg1.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -19,8 +19,8 @@ areal = randn(n,n)/2
aimg = randn(n,n)/2
breal = randn(n,2)/2
bimg = randn(n,2)/2
for eltya in (Float32, Float64, Complex32, Complex64, Complex128, BigFloat, Int)
for eltyb in (Float32, Float64, Complex32, Complex64, Complex128, Int)
for eltya in (Float32, Float64, Complex64, Complex128, BigFloat, Int)
for eltyb in (Float32, Float64, Complex64, Complex128, Int)
a = eltya == Int ? rand(1:7, n, n) : convert(Matrix{eltya}, eltya <: Complex ? complex(areal, aimg) : areal)
b = eltyb == Int ? rand(1:5, n, 2) : convert(Matrix{eltyb}, eltyb <: Complex ? complex(breal, bimg) : breal)
asym = a'+a # symmetric indefinite
Expand Down Expand Up @@ -131,8 +131,9 @@ debug && println("symmetric eigen-decomposition")
if eltya != BigFloat && eltyb != BigFloat # Revisit when implemented in julia
d,v = eig(asym)
@test_approx_eq asym*v[:,1] d[1]*v[:,1]
@test_approx_eq v*scale(d,v') asym
@test_approx_eq v*Diagonal(d)*v' asym
@test isequal(eigvals(asym[1]), eigvals(asym[1:1,1:1]))
@test_approx_eq abs(eigfact(Hermitian(asym), 1,2)[:vectors]'v[:,1:2]) eye(eltya, 2)
end

debug && println("non-symmetric eigen decomposition")
Expand Down