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

Add balancing options for eigenvector calculations #5428

Merged
merged 1 commit into from
Jan 25, 2014
Merged

Add balancing options for eigenvector calculations #5428

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
2 changes: 2 additions & 0 deletions NEWS.md
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -108,6 +108,8 @@ Library improvements

* `LinAlg` (linear algebra) improvements

* Balancing options for eigenvector calculations

* Sparse linear algebra

* Faster sparse `kron` ([#4958]).
Expand Down
20 changes: 10 additions & 10 deletions base/linalg/factorization.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -440,12 +440,12 @@ function getindex(A::Union(Eigen,GeneralizedEigen), d::Symbol)
throw(KeyError(d))
end

function eigfact!{T<:BlasReal}(A::StridedMatrix{T})
function eigfact!{T<:BlasReal}(A::StridedMatrix{T}; balance::Symbol=:balance)
n = size(A, 2)
n==0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
issym(A) && return eigfact!(Symmetric(A))

WR, WI, VL, VR = LAPACK.geev!('N', 'V', A)
A, WR, WI, VL, VR, _ = LAPACK.geevx!(balance == :balance ? 'B' : (balance == :diagonal ? 'S' : (balance == :permute ? 'P' : (balance == :nobalance ? 'N' : throw(ArgumentError("balance must be either ':balance', ':diagonal', ':permute' or ':nobalance'"))))), 'N', 'V', 'N', A)
all(WI .== 0.) && return Eigen(WR, VR)
evec = zeros(Complex{T}, n, n)
j = 1
Expand All @@ -462,24 +462,24 @@ function eigfact!{T<:BlasReal}(A::StridedMatrix{T})
return Eigen(complex(WR, WI), evec)
end

function eigfact!{T<:BlasComplex}(A::StridedMatrix{T})
function eigfact!{T<:BlasComplex}(A::StridedMatrix{T}; balance::Symbol=:balance)
n = size(A, 2)
n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
ishermitian(A) && return eigfact!(Hermitian(A))
return Eigen(LAPACK.geev!('N', 'V', A)[[1,3]]...)
return Eigen(LAPACK.geevx!(balance == :balance ? 'B' : (balance == :diagonal ? 'S' : (balance == :permute ? 'P' : (balance == :nobalance ? 'N' : throw(ArgumentError("balance must be either ':balance', 'diagonal', 'permute' or 'nobalance'"))))), 'N', 'V', 'N', A)[[2,4]]...)
end
eigfact!(A::StridedMatrix) = eigfact!(float(A))
eigfact{T<:BlasFloat}(x::StridedMatrix{T}) = eigfact!(copy(x))
eigfact(A::StridedMatrix) = eigfact!(float(A))
eigfact!(A::StridedMatrix; balance::Symbol=:balance) = eigfact!(float(A), balance=balance)
eigfact{T<:BlasFloat}(x::StridedMatrix{T}; balance::Symbol=:balance) = eigfact!(copy(x), balance=balance)
eigfact(A::StridedMatrix; balance::Symbol=:balance) = eigfact!(float(A), balance=balance)
eigfact(x::Number) = Eigen([x], fill(one(x), 1, 1))

function eig(A::Union(Number, AbstractMatrix))
F = eigfact(A)
function eig(A::Union(Number, AbstractMatrix), args...)
F = eigfact(A, args...)
F[:values], F[:vectors]
end

#Calculates eigenvectors
eigvecs(A::Union(Number, AbstractMatrix)) = eigfact(A)[:vectors]
eigvecs(A::Union(Number, AbstractMatrix); balance::Symbol=:balance) = eigfact(A, balance=balance)[:vectors]

function eigvals!{T<:BlasReal}(A::StridedMatrix{T})
issym(A) && return eigvals!(Symmetric(A))
Expand Down
121 changes: 115 additions & 6 deletions base/linalg/lapack.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -1167,10 +1167,65 @@ for (geev, gesvd, gesdd, ggsvd, elty, relty) in
end
end
end
for (ggev, elty) in
((:dggev_,:Float64),
(:sggev_,:Float32))
## Expert driver and generalized eigenvlua problem
for (geevx, ggev, elty) in
((:dgeevx_,:dggev_,:Float64),
(:sgeevx_,:sggev_,:Float32))
@eval begin
# SUBROUTINE DGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, WR, WI,
# VL, LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM,
# RCONDE, RCONDV, WORK, LWORK, IWORK, INFO )
#
# .. Scalar Arguments ..
# CHARACTER BALANC, JOBVL, JOBVR, SENSE
# INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
# DOUBLE PRECISION ABNRM
# ..
# .. Array Arguments ..
# INTEGER IWORK( * )
# DOUBLE PRECISION A( LDA, * ), RCONDE( * ), RCONDV( * ),
# $ SCALE( * ), VL( LDVL, * ), VR( LDVR, * ),
# $ WI( * ), WORK( * ), WR( * )
function geevx!(balanc::Char, jobvl::Char, jobvr::Char, sense::Char, A::StridedMatrix{$elty})
n = chksquare(A)
lda = max(1,stride(A,2))
wr = Array($elty, n)
wi = Array($elty, n)
ldvl = jobvl == 'V' ? n : (jobvl == 'N' ? 0 : throw(ArgumentError("jobvl must be 'V' or 'N'")))
VL = Array($elty, ldvl, n)
ldvr = jobvr == 'V' ? n : (jobvr == 'N' ? 0 : throw(ArgumentError("jobvr must be 'V' or 'N'")))
VR = Array($elty, ldvr, n)
ilo = Array(BlasInt, 1)
ihi = Array(BlasInt, 1)
scale = Array($elty, n)
abnrm = Array($elty, 1)
rconde = Array($elty, n)
rcondv = Array($elty, n)
work = Array($elty, 1)
lwork::BlasInt = -1
iwork = Array(BlasInt, sense == 'N' || sense == 'E' ? 0 : (sense == 'V' || sense == 'B' ? 2n-2 : throw(ArgumentError("argument sense must be 'N', 'E', 'V' or 'B'"))))
info = Array(BlasInt, 1)
for i = 1:2
ccall(($(string(geevx)),Base.liblapack_name), Void,
(Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8},
Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
Ptr{$elty}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$elty},
Ptr{$elty}, Ptr{$elty}, Ptr{$elty}, Ptr{$elty},
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{BlasInt}),
&balanc, &jobvl, &jobvr, &sense,
&n, A, &lda, wr,
wi, VL, &max(1,ldvl), VR,
&max(1,ldvr), ilo, ihi, scale,
abnrm, rconde, rcondv, work,
&lwork, iwork, info)
lwork = convert(BlasInt, work[1])
work = Array($elty, lwork)
end
@lapackerror
A, wr, wi, VL, VR, ilo[1], ihi[1], scale, abnrm[1], rconde, rcondv
end

# SUBROUTINE DGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHAR, ALPHAI,
# $ BETA, VL, LDVL, VR, LDVR, WORK, LWORK, INFO )
# * .. Scalar Arguments ..
Expand Down Expand Up @@ -1219,10 +1274,63 @@ for (ggev, elty) in
end
end
end
for (ggev, elty, relty) in
((:zggev_,:Complex128,:Float64),
(:cggev_,:Complex64,:Float32))
for (geevx, ggev, elty, relty) in
((:zgeevx_,:zggev_,:Complex128,:Float64),
(:cgeevx_,:cggev_,:Complex64,:Float32))
@eval begin
# SUBROUTINE ZGEEVX( BALANC, JOBVL, JOBVR, SENSE, N, A, LDA, W, VL,
# LDVL, VR, LDVR, ILO, IHI, SCALE, ABNRM, RCONDE,
# RCONDV, WORK, LWORK, RWORK, INFO )
#
# .. Scalar Arguments ..
# CHARACTER BALANC, JOBVL, JOBVR, SENSE
# INTEGER IHI, ILO, INFO, LDA, LDVL, LDVR, LWORK, N
# DOUBLE PRECISION ABNRM
# ..
# .. Array Arguments ..
# DOUBLE PRECISION RCONDE( * ), RCONDV( * ), RWORK( * ),
# $ SCALE( * )
# COMPLEX*16 A( LDA, * ), VL( LDVL, * ), VR( LDVR, * ),
# $ W( * ), WORK( * )
function geevx!(balanc::Char, jobvl::Char, jobvr::Char, sense::Char, A::StridedMatrix{$elty})
n = chksquare(A)
lda = max(1,stride(A,2))
w = Array($elty, n)
ldvl = jobvl == 'V' ? n : (jobvl == 'N' ? 0 : throw(ArgumentError("jobvl must be 'V' or 'N'")))
VL = Array($elty, ldvl, n)
ldvr = jobvr == 'V' ? n : (jobvr == 'N' ? 0 : throw(ArgumentError("jobvr must be 'V' or 'N'")))
VR = Array($elty, ldvr, n)
ilo = Array(BlasInt, 1)
ihi = Array(BlasInt, 1)
scale = Array($relty, n)
abnrm = Array($relty, 1)
rconde = Array($relty, n)
rcondv = Array($relty, n)
work = Array($elty, 1)
lwork::BlasInt = -1
rwork = Array($relty, 2n)
info = Array(BlasInt, 1)
for i = 1:2
ccall(($(string(geevx)),Base.liblapack_name), Void,
(Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8}, Ptr{Uint8},
Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty},
Ptr{$elty}, Ptr{BlasInt}, Ptr{$elty}, Ptr{BlasInt},
Ptr{BlasInt}, Ptr{BlasInt}, Ptr{$relty}, Ptr{$relty},
Ptr{$relty}, Ptr{$relty}, Ptr{$elty}, Ptr{BlasInt},
Ptr{$relty}, Ptr{BlasInt}),
&balanc, &jobvl, &jobvr, &sense,
&n, A, &lda, w,
VL, &max(1,ldvl), VR, &max(1,ldvr),
ilo, ihi, scale, abnrm,
rconde, rcondv, work, &lwork,
rwork, info)
lwork = convert(BlasInt, work[1])
work = Array($elty, lwork)
end
@lapackerror
A, w, VL, VR, ilo[1], ihi[1], scale, abnrm[1], rconde, rcondv
end

# SUBROUTINE ZGGEV( JOBVL, JOBVR, N, A, LDA, B, LDB, ALPHA, BETA,
# $ VL, LDVL, VR, LDVR, WORK, LWORK, RWORK, INFO )
# * .. Scalar Arguments ..
Expand Down Expand Up @@ -1271,6 +1379,7 @@ for (ggev, elty, relty) in
end
end
end

# One step incremental condition estimation of max/min singular values
for (laic1, elty) in
((:dlaic1_,:Float64),
Expand Down
10 changes: 5 additions & 5 deletions doc/stdlib/linalg.rst
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -97,9 +97,9 @@ Linear algebra functions in Julia are largely implemented by calling functions f

Compute the matrix square root of ``A``. If ``B = sqrtm(A)``, then ``B*B == A`` within roundoff error.

.. function:: eig(A) -> D, V
.. function:: eig(A,[balance=:balance]) -> D, V

Compute eigenvalues and eigenvectors of A
Compute eigenvalues and eigenvectors of A. See ``eigfact`` for details on the ``balance`` keyword argument.

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

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

Returns the smallest eigenvalue of ``A``.

.. function:: eigvecs(A, [eigvals])
.. function:: eigvecs(A, [eigvals,][balance=:balance])

Returns the eigenvectors of ``A``.

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

.. function:: eigfact(A)
.. function:: eigfact(A,[balance=:balance])

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``.
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. Possible values are ``:balance``(default), ``:permute``,``:diagonal`` and ``:nobalance``.

.. function:: eigfact(A, B)

Expand Down
13 changes: 13 additions & 0 deletions test/linalg.jl
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -847,6 +847,19 @@ Hinv = Rational{BigInt}[(-1)^(i+j)*(i+j-1)*binomial(nHilbert+i-1,nHilbert-j)*bin
with_bigfloat_precision(2^10) do
@test norm(float64(inv(float(H)) - float(Hinv))) < 1e-100
end

# Test balancing in eigenvector calculations
for elty in (Float32, Float64, Complex64, Complex128)
A = convert(Matrix{elty}, [ 3.0 -2.0 -0.9 2*eps(real(one(elty)));
-2.0 4.0 1.0 -eps(real(one(elty)));
-eps(real(one(elty)))/4 eps(real(one(elty)))/2 -1.0 0;
-0.5 -0.5 0.1 1.0])
F = eigfact(A,balance=:nobalance)
@test_approx_eq F[:vectors]*Diagonal(F[:values])/F[:vectors] A
F = eigfact(A)
@test norm(F[:vectors]*Diagonal(F[:values])/F[:vectors] - A) > 0.01
end

## Issue related tests
# issue 1447
let
Expand Down