From 22117ae39b8f0b9dd164966ef54871a5a3caf46d Mon Sep 17 00:00:00 2001 From: Andreas Noack Jensen Date: Sat, 26 Apr 2014 17:02:20 +0200 Subject: [PATCH] Add arguments specifying subsets of eigenvalues in eigfact(Symmetric/Hermitian) --- base/linalg/factorization.jl | 9 ++------- base/linalg/symmetric.jl | 13 ++++++------- doc/stdlib/linalg.rst | 12 ++++++------ test/linalg1.jl | 7 ++++--- 4 files changed, 18 insertions(+), 23 deletions(-) diff --git a/base/linalg/factorization.jl b/base/linalg/factorization.jl index 4037fc2b54664..ff073e727a908 100644 --- a/base/linalg/factorization.jl +++ b/base/linalg/factorization.jl @@ -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 @@ -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) diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index 2e073c37eb0cb..eada932c25a24 100644 --- a/base/linalg/symmetric.jl +++ b/base/linalg/symmetric.jl @@ -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] diff --git a/doc/stdlib/linalg.rst b/doc/stdlib/linalg.rst index 56c7dff7657f0..8ac94cf8c839f 100644 --- a/doc/stdlib/linalg.rst +++ b/doc/stdlib/linalg.rst @@ -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) @@ -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) diff --git a/test/linalg1.jl b/test/linalg1.jl index a0d97c04754a9..0cba9ba3f2f5e 100644 --- a/test/linalg1.jl +++ b/test/linalg1.jl @@ -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 @@ -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")