From 563dee9ea0b7069fd36ceea720ef4926e656d148 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Wed, 26 Nov 2025 17:00:16 -0500 Subject: [PATCH 1/4] use syevr meta-algorithm for SymTridiagonal --- src/tridiag.jl | 57 +++++++++++++++++++++++++++++++++++++++++++------ test/tridiag.jl | 11 ++++++++++ 2 files changed, 62 insertions(+), 6 deletions(-) diff --git a/src/tridiag.jl b/src/tridiag.jl index 2357ecc1..1de09294 100644 --- a/src/tridiag.jl +++ b/src/tridiag.jl @@ -278,29 +278,74 @@ end ldiv!(A::SymTridiagonal, B::AbstractVecOrMat; shift::Number=false) = ldiv!(ldlt(A, shift=shift), B) rdiv!(B::AbstractVecOrMat, A::SymTridiagonal; shift::Number=false) = rdiv!(B, ldlt(A, shift=shift)) -eigen!(A::SymTridiagonal{<:BlasReal,<:StridedVector}) = Eigen(LAPACK.stegr!('V', A.dv, A.ev)...) +# tridiagonal eigensolver meta-algorithm from LAPACK.syevr! for alg==RobustRepresentations() +# - if all eigenvalues are desired, call stev == sterf (eigvals) or stegr == stemr (eigen) +# - otherwise, and also if an error occurs, fall back to stebz and (if eigvecs wanted) stein +function syevr_tri_eigen(range::AbstractChar, dv::AbstractVector{T}, ev::AbstractVector{T}, vl::Real, vu::Real, il::Integer, iu::Integer) where {T<:BlasReal} + if range == 'A' || (range == 'I' && il == 1 && iu == length(dv)) + try + # need to copy dv, ev so that they are available for fallbacks, below + return Eigen(LAPACK.stegr!('V', range, copymutable(dv), copymutable(ev), vl, vu, il, iu)...) + catch ex + ex isa LAPACKException || rethrow() + end + end + # note that these functions do not actually modify dv, ev, despite the ! + values, iblock, isplit = LAPACK.stebz!(range, 'B', vl, vu, il, iu, -1.0, dv, ev) + vectors = LAPACK.stein!(dv, ev, values, iblock, isplit) + return Eigen(values, vectors) +end +function syevr_tri_eigvals(range::AbstractChar, dv::AbstractVector{T}, ev::AbstractVector{T}, vl::Real, vu::Real, il::Integer, iu::Integer) where {T<:BlasReal} + if range == 'A' || (range == 'I' && il == 1 && iu == length(dv)) + try + # need to copy dv, ev so that they are available for fallbacks, below + return LAPACK.stev!('N', copymutable(dv), copymutable(ev))[1] + catch ex + ex isa LAPACKException || rethrow() + end + end + # note that this function does not actually modify dv, ev, despite the ! + values, iblock, isplit = LAPACK.stebz!(range, 'B', vl, vu, il, iu, -1.0, dv, ev) + return values +end +syevr_tri_eigen(dv::AbstractVector{T}, ev::AbstractVector{T}) where {T<:BlasReal} = + syevr_tri_eigen('A', dv, ev, 0.0, 0.0, 0, 0) +syevr_tri_eigvals(dv::AbstractVector{T}, ev::AbstractVector{T}) where {T<:BlasReal} = + syevr_tri_eigvals('A', dv, ev, 0.0, 0.0, 0, 0) + +eigen!(A::SymTridiagonal{<:BlasReal,<:StridedVector}) = syevr_tri_eigen(A.dv, A.ev) +eigen(A::SymTridiagonal{<:BlasReal,<:StridedVector}) = syevr_tri_eigen(A.dv, A.ev) eigen(A::SymTridiagonal{T}) where T = eigen!(copymutable_oftype(A, eigtype(T))) eigen!(A::SymTridiagonal{<:BlasReal,<:StridedVector}, irange::UnitRange) = - Eigen(LAPACK.stegr!('V', 'I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop)...) + syevr_tri_eigen('I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop) +eigen(A::SymTridiagonal{<:BlasReal,<:StridedVector}, irange::UnitRange) = + syevr_tri_eigen('I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop) eigen(A::SymTridiagonal{T}, irange::UnitRange) where T = eigen!(copymutable_oftype(A, eigtype(T)), irange) eigen!(A::SymTridiagonal{<:BlasReal,<:StridedVector}, vl::Real, vu::Real) = - Eigen(LAPACK.stegr!('V', 'V', A.dv, A.ev, vl, vu, 0, 0)...) + syevr_tri_eigen('V', A.dv, A.ev, vl, vu, 0, 0) +eigen(A::SymTridiagonal{<:BlasReal,<:StridedVector}, vl::Real, vu::Real) = + syevr_tri_eigen('V', A.dv, A.ev, vl, vu, 0, 0) eigen(A::SymTridiagonal{T}, vl::Real, vu::Real) where T = eigen!(copymutable_oftype(A, eigtype(T)), vl, vu) -eigvals!(A::SymTridiagonal{<:BlasReal,<:StridedVector}) = LAPACK.stev!('N', A.dv, A.ev)[1] +eigvals!(A::SymTridiagonal{<:BlasReal,<:StridedVector}) = syevr_tri_eigvals(A.dv, A.ev) +eigvals(A::SymTridiagonal{<:BlasReal,<:StridedVector}) = syevr_tri_eigvals(A.dv, A.ev) eigvals(A::SymTridiagonal{T}) where T = eigvals!(copymutable_oftype(A, eigtype(T))) eigvals!(A::SymTridiagonal{<:BlasReal,<:StridedVector}, irange::UnitRange) = - LAPACK.stegr!('N', 'I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop)[1] + syevr_tri_eigvals('I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop) +eigvals(A::SymTridiagonal{<:BlasReal,<:StridedVector}, irange::UnitRange) = + syevr_tri_eigvals('I', A.dv, A.ev, 0.0, 0.0, irange.start, irange.stop) eigvals(A::SymTridiagonal{T}, irange::UnitRange) where T = eigvals!(copymutable_oftype(A, eigtype(T)), irange) eigvals!(A::SymTridiagonal{<:BlasReal,<:StridedVector}, vl::Real, vu::Real) = - LAPACK.stegr!('N', 'V', A.dv, A.ev, vl, vu, 0, 0)[1] + syevr_tri_eigvals('V', A.dv, A.ev, vl, vu, 0, 0) +eigvals(A::SymTridiagonal{<:BlasReal,<:StridedVector}, vl::Real, vu::Real) = + syevr_tri_eigvals('V', A.dv, A.ev, vl, vu, 0, 0) eigvals(A::SymTridiagonal{T}, vl::Real, vu::Real) where T = eigvals!(copymutable_oftype(A, eigtype(T)), vl, vu) diff --git a/test/tridiag.jl b/test/tridiag.jl index 7a4d78b9..7678ee98 100644 --- a/test/tridiag.jl +++ b/test/tridiag.jl @@ -1245,4 +1245,15 @@ end end end +@testset "issue #1491" begin + # example where dstemr fails and we need to fall back to dstebz + dstein + mat = SymTridiagonal([1.9081958235744878e-17, 2.4936649967166602e-17, 0.07337398524939442, -0.037808419232506454, -0.019817057861955256, 0.02058130717549595, 0.0954455291925313, -0.14462978163628976, -0.03603900926534216, 0.10574604437162191, 0.0032290029813105, -0.060091988785585734, 0.15705171399550666, -0.14949686924737385, -0.07781372420855004, 0.03553277340671564, -0.016676626898563265, 0.041446049523777104, 0.10027425296079334, -0.10113668506693083, -0.1489473358740452, 0.1254312751044567, 0.34898062637027333, -0.33073386964350204, 0.2590576535805837, -0.25350583802618526, 0.43822601856881893, -0.4307646088223485, 0.779796594155107, -0.8224956491946511, 0.9089301791662714, -0.8619090557448831, 0.2318453071607881, -0.23320912386320142, 1.3140530211788515, -1.3139256900796155, -0.08390903751033216, 0.08392205661902721, 0.22240033452483454, -0.22246103504490297, 0.029954811490418787, -0.027812830389636778, 0.5263757772376259, -0.5365156423835895, 0.39737424822653156, -0.37936748640222284, 0.9333689384817099, -0.9703849367869426, 1.4044213669711443, -1.3377550567558987, 0.13172978466921414, -0.13519350037871353, 0.4999843222780838, -0.5567430191155477, 1.581258524837276, -1.5341564150622489, 0.8586746358146186, -0.8739294540224967, -0.17581567809216825, 0.15338653603265703, 0.780052785514151, -0.7342698870670235, -0.3173535579981542, 1.3364975241689332, -2.007075402840063, 1.0579643498230522, 0.31139215377118973, -0.39386762936936015, 1.8984071979199968, -1.920544791889951, 0.49358795540233863, -0.4613352087464951, 2.455277289025621, -2.4875372089056578, 0.9570609479497718, -0.9570638796447994, -1.1779359208299722, 1.163074811976796, -1.9301179244054274, 1.9450298141233324, -1.1378646951262403, 0.6350150906475753, 0.08802515616268328, 0.44163280306776015, -0.32773792909920113, 0.32952024534152274, 1.4740728008471815, -1.5038699938845244, 0.6847048634935279, -0.6874625497241057, 1.8751566802154134, -1.8547286700155103, 1.2640098482403956, -1.2174383863708917], [4.0208830900619335, 2.571546145080502e-15, 3.6768867672571894, 2.1275700268385056, 2.6591525597707077, 1.5630051292577454, 3.169596676339297, 1.4476618215762946, 3.1868304175740785, 1.424991114261858, 3.4173704890401693, 1.081122203108007, 3.8388891688820155, 0.7188503200238101, 3.628959037595643, 1.2527365479664625, 3.655577511598652, 0.5795107784623137, 3.9154918603711404, 0.7305602983115349, 4.204984775263665, 0.4667273034707952, 4.086056779436833, 0.7605098745710153, 3.765995021766351, 0.684123305522116, 4.180477278972855, 0.22109649225044853, 4.1892495574511095, 0.5187290599147804, 4.267105890594963, 0.38956213309842774, 4.395771702058166, 0.012968327575103515, 4.370023869671557, 9.368149211179955e-5, 1.7478905204265782, 0.012088046413693718, 2.5400351007261417, 0.02355280591484052, 2.966532136162213, 0.1157714633509318, 3.4610631638791225, 0.21831668966320836, 3.6404645575284627, 0.2911241355881879, 4.444637765978322, 0.46127544435962, 3.9217404428871023, 0.4673958058890463, 3.6185829343499414, 0.6185854819553929, 4.063721954940776, 0.4953177806242663, 4.017151353166678, 0.4035551805802303, 4.318876887238373, 0.17260927157837835, 4.41629817034136, 0.251496875950196, 4.017922104317157, 0.47982184361329844, 3.1370471883575424, 2.4334652498110247, 1.9981994742982454, 0.8095271305424269, 4.40837268034087, 0.18738580231835025, 4.147189662985785, 0.02664936509207341, 2.521497611392047, 0.4752292595225756, 3.795983764313679, 0.004636199088478572, 2.8022242947838416, 0.0004261020805263317, 3.2721199558188827, 0.16661903081014354, 3.1196340292912215, 0.03715487568175004, 2.6328978561492336, 1.5366796979500048, 2.1884275437040994, 0.7509912208995596, 4.056762760327159, 0.3470101947911033, 4.267342895676449, 0.4303372808790385, 4.051346888184243, 0.14621115143865943, 4.2045273832895145, 0.2540956144073191, 4.064260328749705]) + λ = [-4.683759859765017, -4.683759859765017, -4.680622745270949, -4.680066202102797, -4.680066202102788, -4.642387544636312, -4.642387544636312, -4.62850079279404, -4.5645761799521845, -4.564576179952184, -4.563095222112828, -4.563095222112824, -4.4746776824461, -4.474677682446098, -4.372795904675404, -4.372795904675403, -4.32974284978971, -4.329742849789709, -4.190416629535501, -4.160130115632619, -4.160130115632617, -4.158746679334386, -4.020883090061934, -4.0208830900619335, -4.0208830900619335, -3.992984903781014, -3.992984903781012, -3.992659934264744, -3.699861925653249, -3.699861925653247, -3.699861925580689, -3.4509279319955604, -3.4509279319955577, -3.4509279319623327, -3.4082484759685445, -3.4082484759685436, -3.4082484739031704, -2.9611516414045025, -2.9611516414045007, -2.9611516414045003, -2.5495001812022124, -2.549500181202207, -2.5495001812022027, -1.7498663309027835, -1.7498663308995552, -1.749866330899554, -1.749866330899551, 1.7498663308995548, 1.7498663308995575, 1.749866330899561, 1.7498663309055233, 2.549500181202202, 2.5495001812022067, 2.5495001812022093, 2.961151641404497, 2.9611516414045007, 2.9611516414045043, 3.408248467695482, 3.4082484759685476, 3.408248475968552, 3.450927931921516, 3.45092793199556, 3.4509279319955635, 3.6998619253791105, 3.6998619256532437, 3.6998619256532463, 3.99283345851548, 3.9929849037810126, 3.9929849037810143, 4.0208830900619335, 4.0208830900619335, 4.020883090061935, 4.159323562240434, 4.160130115632611, 4.160130115632616, 4.245716444792351, 4.32974284978971, 4.329742849789711, 4.372795904675407, 4.37279590467541, 4.474677682446099, 4.474677682446103, 4.563095222112824, 4.563095222112828, 4.564576179952181, 4.564576179952186, 4.6363803854054515, 4.642387544636308, 4.642387544636311, 4.679778864396294, 4.6800662021027755, 4.680066202102799, 4.6837598597650185, 4.683759859765019] + densemat = Symmetric(Matrix(mat)) + F = eigen(mat) + Fdense = eigen(densemat) + @test λ ≈ eigvals(mat) ≈ eigvals(densemat) ≈ F.values ≈ Fdense.values + @test F.vectors ≈ Fdense.vectors +end + end # module TestTridiagonal From 3c939aeb387387f400ba438ae166edb0a9c14d36 Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 27 Nov 2025 23:49:07 -0500 Subject: [PATCH 2/4] Update src/tridiag.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Mateus Araújo --- src/tridiag.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tridiag.jl b/src/tridiag.jl index 1de09294..17076f0f 100644 --- a/src/tridiag.jl +++ b/src/tridiag.jl @@ -291,7 +291,7 @@ function syevr_tri_eigen(range::AbstractChar, dv::AbstractVector{T}, ev::Abstrac end end # note that these functions do not actually modify dv, ev, despite the ! - values, iblock, isplit = LAPACK.stebz!(range, 'B', vl, vu, il, iu, -1.0, dv, ev) + values, iblock, isplit = LAPACK.stebz!(range, 'B', T(vl), T(vu), il, iu, -1.0, dv, ev) vectors = LAPACK.stein!(dv, ev, values, iblock, isplit) return Eigen(values, vectors) end From af2ce81b9c97d2379b76dabbc31e41ee9a9044da Mon Sep 17 00:00:00 2001 From: "Steven G. Johnson" Date: Thu, 27 Nov 2025 23:49:22 -0500 Subject: [PATCH 3/4] Update src/tridiag.jl MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit Co-authored-by: Mateus Araújo --- src/tridiag.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/tridiag.jl b/src/tridiag.jl index 17076f0f..42bcbc52 100644 --- a/src/tridiag.jl +++ b/src/tridiag.jl @@ -305,7 +305,7 @@ function syevr_tri_eigvals(range::AbstractChar, dv::AbstractVector{T}, ev::Abstr end end # note that this function does not actually modify dv, ev, despite the ! - values, iblock, isplit = LAPACK.stebz!(range, 'B', vl, vu, il, iu, -1.0, dv, ev) + values, iblock, isplit = LAPACK.stebz!(range, 'B', T(vl), T(vu), il, iu, -1.0, dv, ev) return values end syevr_tri_eigen(dv::AbstractVector{T}, ev::AbstractVector{T}) where {T<:BlasReal} = From f5460b56033366e1c32375127e619621610d9bc4 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Mateus=20Ara=C3=BAjo?= Date: Fri, 28 Nov 2025 08:53:34 +0100 Subject: [PATCH 4/4] fix doctest --- src/symmetriceigen.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/symmetriceigen.jl b/src/symmetriceigen.jl index a8a54037..8f850c43 100644 --- a/src/symmetriceigen.jl +++ b/src/symmetriceigen.jl @@ -189,7 +189,7 @@ julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.]) julia> eigvals(A, 2:2) 1-element Vector{Float64}: - 0.9999999999999996 + 1.0 julia> eigvals(A) 3-element Vector{Float64}: @@ -228,7 +228,7 @@ julia> A = SymTridiagonal([1.; 2.; 1.], [2.; 3.]) julia> eigvals(A, -1, 2) 1-element Vector{Float64}: - 1.0000000000000009 + 1.0000000000000002 julia> eigvals(A) 3-element Vector{Float64}: