Skip to content

Commit

Permalink
Revert "Extend Eigen to keep additional information from geevx (#… (
Browse files Browse the repository at this point in the history 
#41578)

* Revert "Extend `Eigen` to keep additional information from `geevx` (#38483)"

This reverts commit 9d3a7c4.

* Remove left eigenvector handling in Float16 eigen methods.
  • Loading branch information
@andreasnoack
andreasnoack committed Jul 15, 2021
1 parent 430e5e0 commit 41ee0fa
Show file tree
Hide file tree
Showing 2 changed files with 25 additions and 105 deletions.
118 changes: 22 additions & 96 deletions stdlib/LinearAlgebra/src/eigen.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -17,7 +17,7 @@ Iterating the decomposition produces the components `F.values` and `F.vectors`.
# Examples
```jldoctest
julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}, Vector{Float64}}
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
1.0
Expand Down Expand Up @@ -47,18 +47,14 @@ julia> vals == F.values && vecs == F.vectors
true
```
"""
struct Eigen{T,V,S<:AbstractMatrix,U<:AbstractVector,R<:AbstractVector} <: Factorization{T}
struct Eigen{T,V,S<:AbstractMatrix,U<:AbstractVector} <: Factorization{T}
values::U
vectors::S
vectorsl::S
unitary::Bool
rconde::R
rcondv::R
Eigen{T,V,S,U,R}(values::AbstractVector{V}, vectors::AbstractMatrix{T}, vectorsl::AbstractMatrix{T}, unitary::Bool, rconde::R, rcondv::R) where {T,V,S,U,R} =
new(values, vectors, vectorsl, unitary, rconde, rcondv)
Eigen{T,V,S,U}(values::AbstractVector{V}, vectors::AbstractMatrix{T}) where {T,V,S,U} =
new(values, vectors)
end
Eigen(values::AbstractVector{V}, vectors::AbstractMatrix{T}, vectorsl=vectors, uni=true, rce=zeros(real(T),0), rcv=zeros(real(T), 0)) where {T,V} =
Eigen{T,V,typeof(vectors),typeof(values),typeof(rce)}(values, vectors, vectorsl, uni, rce, rcv)
Eigen(values::AbstractVector{V}, vectors::AbstractMatrix{T}) where {T,V} =
Eigen{T,V,typeof(vectors),typeof(values)}(values, vectors)

# Generalized eigenvalue problem.
"""
Expand Down Expand Up @@ -137,21 +133,10 @@ function sorteig!(λ::AbstractVector, X::AbstractMatrix, sortby::Union{Function,
if sortby !== nothing && !issorted(λ, by=sortby)
p = sortperm(λ; alg=QuickSort, by=sortby)
permute!(λ, p)
!isempty(X) && Base.permutecols!!(X, copy(p))
Base.permutecols!!(X, p)
end
return λ, X
end
function sorteig!(λ::AbstractVector, X::AbstractMatrix, sortby::Union{Function,Nothing}, Y::AbstractMatrix, rconde::AbstractVector, rcondv::AbstractVector)
if sortby !== nothing && !issorted(λ, by=sortby)
p = sortperm(λ; alg=QuickSort, by=sortby)
permute!(λ, p)
!isempty(rconde) && permute!(rconde, p)
!isempty(rcondv) && permute!(rcondv, p)
!isempty(X) && Base.permutecols!!(X, copy(p))
!isempty(Y) && X !== Y && Base.permutecols!!(Y, p)
end
return λ, X, Y, false, rconde, rcondv
end
sorteig!(λ::AbstractVector, sortby::Union{Function,Nothing}=eigsortby) = sortby === nothing ? λ : sort!(λ, by=sortby)

"""
Expand All @@ -160,32 +145,12 @@ sorteig!(λ::AbstractVector, sortby::Union{Function,Nothing}=eigsortby) = sortby
Same as [`eigen`](@ref), but saves space by overwriting the input `A` (and
`B`), instead of creating a copy.
"""
function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby, jvl::Bool=false, jvr::Bool=true, jce::Bool=false, jcv::Bool=false) where T<:BlasReal
function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T<:BlasReal
n = size(A, 2)
n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
issymmetric(A) && return eigen!(Symmetric(A), sortby=sortby)

balance = permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N')
jobvl = jvl || jce ? 'V' : 'N'
jobvr = jvr || jce ? 'V' : 'N'
sense = jce && jcv ? 'B' : jce ? 'E' : jcv ? 'V' : 'N'
A, WR, WI, VL, VR, _, _, scale, abnrm, rconde, rcondv = LAPACK.geevx!(balance, jobvl, jobvr, sense, A)
if iszero(WI)
evecr = VR
evecl = VL
evals = WR
else
evecr = complexeig(WI, VR)
evecl = complexeig(WI, VL)
evals = complex.(WR, WI)
end
rconde = jce ? inv.(rconde) : zeros(T, 0)
rcondv = jcv ? inv.(rcondv) : zeros(T, 0)
return Eigen(sorteig!(evals, evecr, sortby, evecl, rconde, rcondv)...)
end

function complexeig(WI::Vector{T}, VR::Matrix{T}) where T
n = min(size(VR)...)
A, WR, WI, VL, VR, _ = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)
iszero(WI) && return Eigen(sorteig!(WR, VR, sortby)...)
evec = zeros(Complex{T}, n, n)
j = 1
while j <= n
Expand All @@ -200,19 +165,15 @@ function complexeig(WI::Vector{T}, VR::Matrix{T}) where T
end
j += 1
end
evec
return Eigen(sorteig!(complex.(WR, WI), evec, sortby)...)
end

function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby, jvl::Bool=false, jvr::Bool=true, jce::Bool=false, jcv::Bool=false) where T<:BlasComplex
function eigen!(A::StridedMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T<:BlasComplex
n = size(A, 2)
n == 0 && return Eigen(zeros(T, 0), zeros(T, 0, 0))
ishermitian(A) && return eigen!(Hermitian(A), sortby=sortby)
balance = permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N')
jobvl = jvl || jce ? 'V' : 'N'
jobvr = jvr || jce ? 'V' : 'N'
sense = jce && jcv ? 'B' : jce ? 'E' : jcv ? 'V' : 'N'
A, W, VL, VR, _, _, scale, abnrm, rconde, rcondv = LAPACK.geevx!(balance, jobvl, jobvr, sense, A)
return Eigen(sorteig!(W, VR, sortby, VL, rconde, rcondv)...)
eval, evec = LAPACK.geevx!(permute ? (scale ? 'B' : 'P') : (scale ? 'S' : 'N'), 'N', 'V', 'N', A)[[2,4]]
return Eigen(sorteig!(eval, evec, sortby)...)
end

"""
Expand Down Expand Up @@ -240,7 +201,7 @@ accept a `sortby` keyword.
# Examples
```jldoctest
julia> F = eigen([1.0 0.0 0.0; 0.0 3.0 0.0; 0.0 0.0 18.0])
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}, Vector{Float64}}
Eigen{Float64, Float64, Matrix{Float64}, Vector{Float64}}
values:
3-element Vector{Float64}:
1.0
Expand Down Expand Up @@ -270,21 +231,18 @@ julia> vals == F.values && vecs == F.vectors
true
```
"""
function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby, jvl::Bool=false, jvr::Bool=true, jce::Bool=false, jcv::Bool=false) where T
function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where T
AA = copy_oftype(A, eigtype(T))
isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
return eigen!(AA; permute=permute, scale=scale, sortby=sortby, jvl=jvl, jvr=jvr, jce=jce, jcv=jcv)
return eigen!(AA; permute=permute, scale=scale, sortby=sortby)
end
function eigen(A::AbstractMatrix{T}; permute::Bool=true, scale::Bool=true, sortby::Union{Function,Nothing}=eigsortby) where {T <: Union{Float16,Complex{Float16}}}
AA = copy_oftype(A, eigtype(T))
isdiag(AA) && return eigen(Diagonal(AA); permute=permute, scale=scale, sortby=sortby)
A = eigen!(AA; permute, scale, sortby)
values = convert(AbstractVector{isreal(A.values) ? Float16 : Complex{Float16}}, A.values)
vectors = convert(AbstractMatrix{isreal(A.vectors) ? Float16 : Complex{Float16}}, A.vectors)
vectorsl = convert(AbstractMatrix{isreal(A.vectors) ? Float16 : Complex{Float16}}, A.vectorsl)
rconde = convert(Vector{Float16}, A.rconde)
rcondv = convert(Vector{Float16}, A.rcondv)
return Eigen(values, vectors, vectorsl, A.unitary, rconde, rcondv)
return Eigen(values, vectors)
end
eigen(x::Number) = Eigen([x], fill(one(x), 1, 1))

Expand Down Expand Up @@ -470,19 +428,7 @@ function eigmin(A::Union{Number, AbstractMatrix};
minimum(v)
end

"""
spectral(f, F::Eigen)
Construct a matrix from an eigen-decomposition `F` by applying the function to
the spectrum (diagonal) of `F`.
"""
function spectral(f, A::Eigen)
d = Diagonal(f.(A.values))
v = A.vectors
vd = v * d
A.unitary ? vd * v' : vd / v
end
inv(A::Eigen) = spectral(inv, A)
inv(A::Eigen) = A.vectors * inv(Diagonal(A.values)) / A.vectors
det(A::Eigen) = prod(A.values)

# Generalized eigenproblem
Expand Down Expand Up @@ -672,28 +618,8 @@ function show(io::IO, mime::MIME{Symbol("text/plain")}, F::Union{Eigen,Generaliz
summary(io, F); println(io)
println(io, "values:")
show(io, mime, F.values)
if !isdefined(F, :vectorsl) || (!isempty(F.vectors) && (F.vectors === F.vectorsl || isempty(F.vectorsl)))
println(io, "\nvectors:")
show(io, mime, F.vectors)
else
if !isempty(F.vectors)
println(io, "\nright vectors:")
show(io, mime, F.vectors)
end
if !isempty(F.vectorsl)
println(io, "\nleft vectors:")
show(io, mime, F.vectorsl)
end
end
if isdefined(F, :rconde) && !isempty(F.rconde)
println(io, "\ncondition values:")
show(io, mime, F.rconde)
end
if isdefined(F, :rcondv) && !isempty(F.rcondv)
println(io, "\ncondition vectors:")
show(io, mime, F.rcondv)
end
nothing
println(io, "\nvectors:")
show(io, mime, F.vectors)
end

function Base.hash(F::Eigen, h::UInt)
Expand All @@ -709,7 +635,7 @@ end
# Conversion methods

## Can we determine the source/result is Real? This is not stored in the type Eigen
AbstractMatrix(F::Eigen) = spectral(identity, F)
AbstractMatrix(F::Eigen) = F.vectors * Diagonal(F.values) / F.vectors
AbstractArray(F::Eigen) = AbstractMatrix(F)
Matrix(F::Eigen) = Array(AbstractArray(F))
Array(F::Eigen) = Matrix(F)
12 changes: 3 additions & 9 deletions stdlib/LinearAlgebra/test/eigen.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -186,27 +186,21 @@ end
D32 = eigen(Float32.(C))
F = eigen(complex(C))
F32 = eigen(complex(Float32.(C)))
@test B isa Eigen{Float16, Float16, Matrix{Float16}, Vector{Float16}, Vector{Float16}}
@test B isa Eigen{Float16, Float16, Matrix{Float16}, Vector{Float16}}
@test B.values isa Vector{Float16}
@test B.vectors isa Matrix{Float16}
@test B.vectorsl isa Matrix{Float16}
@test B.values ≈ B32.values
@test B.vectors ≈ B32.vectors
@test B.vectorsl ≈ B32.vectorsl
@test D isa Eigen{ComplexF16, ComplexF16, Matrix{ComplexF16}, Vector{ComplexF16}, Vector{Float16}}
@test D isa Eigen{ComplexF16, ComplexF16, Matrix{ComplexF16}, Vector{ComplexF16}}
@test D.values isa Vector{ComplexF16}
@test D.vectors isa Matrix{ComplexF16}
@test D.vectorsl isa Matrix{ComplexF16}
@test D.values ≈ D32.values
@test D.vectors ≈ D32.vectors
@test D.vectorsl ≈ D32.vectorsl
@test F isa Eigen{ComplexF16, ComplexF16, Matrix{ComplexF16}, Vector{ComplexF16}, Vector{Float16}}
@test F isa Eigen{ComplexF16, ComplexF16, Matrix{ComplexF16}, Vector{ComplexF16}}
@test F.values isa Vector{ComplexF16}
@test F.vectors isa Matrix{ComplexF16}
@test F.vectorsl isa Matrix{ComplexF16}
@test F.values ≈ F32.values
@test F.vectors ≈ F32.vectors
@test F.vectorsl ≈ F32.vectorsl
end

end # module TestEigen

0 comments on commit 41ee0fa

Please sign in to comment.