Skip to content

Commit

Permalink
Avoid stack overflow in triangular eigvecs (#55497)
Browse files Browse the repository at this point in the history
This fixes a stack overflow in 
```julia
julia> using LinearAlgebra, StaticArrays

julia> U = UpperTriangular(SMatrix{2,2}(1:4))
2×2 UpperTriangular{Int64, SMatrix{2, 2, Int64, 4}} with indices SOneTo(2)×SOneTo(2):
 1  3
 ⋅  4

julia> eigvecs(U)
Warning: detected a stack overflow; program state may be corrupted, so further execution might be unreliable.
ERROR: StackOverflowError:
Stacktrace:
 [1] eigvecs(A::UpperTriangular{Float32, SMatrix{2, 2, Float32, 4}}) (repeats 79984 times)
   @ LinearAlgebra ~/.julia/juliaup/julia-nightly/share/julia/stdlib/v1.12/LinearAlgebra/src/triangular.jl:2749
```
After this,
```julia
julia> eigvecs(U)
2×2 Matrix{Float32}:
 1.0  1.0
 0.0  1.0
```
  • Loading branch information
jishnub authored and KristofferC committed Sep 12, 2024
1 parent e67fe36 commit 7b46791
Show file tree
Hide file tree
Showing 2 changed files with 24 additions and 0 deletions.
8 changes: 8 additions & 0 deletions stdlib/LinearAlgebra/src/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -2782,6 +2782,14 @@ end

# Generic eigensystems
eigvals(A::AbstractTriangular) = diag(A)
# fallback for unknown types
function eigvecs(A::AbstractTriangular{<:BlasFloat})
if istriu(A)
eigvecs(UpperTriangular(Matrix(A)))
else # istril(A)
eigvecs(LowerTriangular(Matrix(A)))
end
end
function eigvecs(A::AbstractTriangular{T}) where T
TT = promote_type(T, Float32)
if TT <: BlasFloat
Expand Down
16 changes: 16 additions & 0 deletions stdlib/LinearAlgebra/test/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1198,6 +1198,22 @@ end
end
end

@testset "eigvecs for AbstractTriangular" begin
S = SizedArrays.SizedArray{(3,3)}(reshape(1:9,3,3))
for T in (UpperTriangular, UnitUpperTriangular,
LowerTriangular, UnitLowerTriangular)
U = T(S)
V = eigvecs(U)
λ = eigvals(U)
@test U * V V * Diagonal(λ)

MU = MyTriangular(U)
V = eigvecs(U)
λ = eigvals(U)
@test MU * V V * Diagonal(λ)
end
end

@testset "(l/r)mul! and (l/r)div! for generic triangular" begin
@testset for T in (UpperTriangular, LowerTriangular, UnitUpperTriangular, UnitLowerTriangular)
M = MyTriangular(T(rand(4,4)))
Expand Down

0 comments on commit 7b46791

Please sign in to comment.