From 8cbba67b7b67a728b1cb1a773f7d82916e6bc451 Mon Sep 17 00:00:00 2001 From: MichelJuillard Date: Sat, 3 Dec 2022 13:46:58 +0100 Subject: [PATCH] fix bug in general schur eigenvalues computation --- src/schur.jl | 20 +++++++++++--------- test/schur_test.jl | 20 ++++++++++++-------- 2 files changed, 23 insertions(+), 17 deletions(-) diff --git a/src/schur.jl b/src/schur.jl index 9f496ae..a2c02b8 100644 --- a/src/schur.jl +++ b/src/schur.jl @@ -363,27 +363,29 @@ for (gges, elty) in ((:dgges_, :Float64), end chklapackerror(info[]) @inbounds for i in axes(A, 1) - ws.eigen_values[i] = complex(ws.αr[i], ws.αi[i]) + ws.eigen_values[i] = complex(ws.αr[i], ws.αi[i])/ws.β[i] end - return A, B, ws.eigen_values, ws.β, ws.vsl, ws.vsr + return A, B, complex.(ws.αr, ws.αi), ws.β, ws.vsl, ws.vsr end end end """ - gges!(ws, jobvsl, jobvsr, A, B; select=nothing, resize=true) -> (A, B, ws.eigen_values, ws.β, ws.vsl, ws.vsr) + gges!(ws, jobvsl, jobvsr, A, B; select=nothing, resize=true) -> (A, B, ws.α, ws.β, ws.vsl, ws.vsr) Computes the generalized eigenvalues, generalized Schur form, left Schur vectors (`jobsvl = V`), or right Schur vectors (`jobvsr = V`) of `A` and `B`, using preallocated [`GeneralizedSchurWs`](@ref) workspace `ws`. If `ws` is not of the right size, and `resize==true` it will be resized appropriately. -It is possible to specify `select`, a function used to sort the eigenvalues during the decomposition. -The function should have the signature `f(αr::T, αi::T, β::T) -> Bool`, where -`αr` and `αi` are the real and imaginary parts of the eigenvalue, `β` the factor, and `T == eltype(A). - -The generalized eigenvalues are returned in `ws.eigen_values` and `ws.β`. The left Schur -vectors are returned in `ws.vsl` and the right Schur vectors are returned in `ws.vsr`. +It is possible to specify `select`, a function used to sort the eigenvalues to the top left corner of the Schur form. +The function should have the signature `f(αr::T, αi::T, β::T) -> Bool` where `T == eltype(A)`. +An eigenvalue `(αr[j]+αi[j])/β[j]` is selected if `f(αr[j],αi[j],β[j])` is true, +i.e. if either one of a complex conjugate pair of eigenvalues is selected, +then both complex eigenvalues are selected. +The generalized eigenvalues components are returned in `ws.α` and `ws.β` where `ws.α` is a complex vector and `ẁs.β`, a real vector. +The generalized eigenvalues (`ws.α./ws.β`) are returned in `ws.eigen_values`, a complex vector. +The left Schur vectors are returned in `ws.vsl` and the right Schur vectors are returned in `ws.vsr`. """ gges!(ws::GeneralizedSchurWs, jobvsl::AbstractChar, jobvsr::AbstractChar, A::AbstractMatrix, B::AbstractMatrix) diff --git a/test/schur_test.jl b/test/schur_test.jl index e7c3519..ef842a1 100644 --- a/test/schur_test.jl +++ b/test/schur_test.jl @@ -53,24 +53,26 @@ end B0 = randn(n, n) ws = GeneralizedSchurWs(copy(A0)) - A1, B1, eig1, β1, vsl1, vsr1 = LAPACK.gges!('V', 'V', copy(A0), copy(B0)) - A2, B2, eig2, β2, vsl2, vsr2 = LAPACK.gges!(ws, 'V', 'V', copy(A0), copy(B0)) + A1, B1, α1, β1, vsl1, vsr1 = LAPACK.gges!('V', 'V', copy(A0), copy(B0)) + A2, B2, α2, β2, vsl2, vsr2 = LAPACK.gges!(ws, 'V', 'V', copy(A0), copy(B0)) @test isapprox(A1, A2) @test isapprox(B1, B2) - @test isapprox(eig1, eig2) + @test isapprox(α1, α2) @test isapprox(β1, β2) @test isapprox(vsl1, vsl2) @test isapprox(vsr1, vsr2) # Using Workspace, factorize! ws = Workspace(LAPACK.gges!, copy(A0)) - A2, B2, eig2, β2, vsl2, vsr2 = factorize!(ws, 'V', 'V', copy(A0), copy(B0)) + A2, B2, α2, β2, vsl2, vsr2 = factorize!(ws, 'V', 'V', copy(A0), copy(B0)) @test isapprox(A1, A2) @test isapprox(B1, B2) - @test isapprox(eig1, eig2) + @test isapprox(α1, α2) @test isapprox(β1, β2) @test isapprox(vsl1, vsl2) @test isapprox(vsr1, vsr2) + @test isapprox(sort(abs.(eigen(A0, B0).values)), sort(abs.(ws.eigen_values))) + show(devnull, "text/plain", ws) for div in (-1,1) @@ -91,14 +93,16 @@ end 0.86268 -0.212549 -0.211994] ws = GeneralizedSchurWs(copy(A0)) - A0, B0, eig, β, vsl, vsr = LAPACK.gges!(ws, 'V', 'V', copy(A0), + A0, B0, α, β, vsl, vsr = LAPACK.gges!(ws, 'V', 'V', copy(A0), copy(B0); select = (ar, ai, b) -> ar^2 + ai^2 < FastLapackInterface.SCHUR_CRITERIUM * b^2) @test ws.sdim[] == 1 - @test sign(real(eig[1])) == -1 - @test sign(real(eig[2])) == 1 + @test (real(α[1])/β[1])^2 < FastLapackInterface.SCHUR_CRITERIUM + @test (real(α[2])/β[2])^2 > FastLapackInterface.SCHUR_CRITERIUM + @test (real(α[3])/β[3])^2 > FastLapackInterface.SCHUR_CRITERIUM + @test isapprox(sort(abs.(eigen(A0, B0).values)), sort(abs.(ws.eigen_values))) show(devnull, "text/plain", ws) end end