Skip to content

Commit

Permalink
Don't access parent of triangular matrix in powm (#52583)
Browse files Browse the repository at this point in the history
Since the values stored in the parent corresponding to the structural
zeros of a tridiagonal matrix aren't well-defined, using it in `ldiv!`
is a footgun that may lead to heisenbugs (one seen in
https://buildkite.com/julialang/julia-master/builds/31285#018c7cc7-6c77-41ac-a01b-1c7d14cb1b15).
This PR changes it to using the tridiagonal matrix directly in `ldiv!`,
which should lead to predictable results, and be bug-free. The failing
tests for #52571 pass locally with this change.

(cherry picked from commit ef549ae)
  • Loading branch information
jishnub authored and KristofferC committed Jan 5, 2024
1 parent 3120989 commit a2e4054
Show file tree
Hide file tree
Showing 2 changed files with 9 additions and 4 deletions.
8 changes: 4 additions & 4 deletions stdlib/LinearAlgebra/src/triangular.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1551,7 +1551,7 @@ end
# Higham and Lin, "An improved Schur-Padé algorithm for fractional powers of
# a matrix and their Fréchet derivatives", SIAM. J. Matrix Anal. & Appl.,
# 34(3), (2013) 1341–1360.
function powm!(A0::UpperTriangular{<:BlasFloat}, p::Real)
function powm!(A0::UpperTriangular, p::Real)
if abs(p) >= 1
throw(ArgumentError("p must be a real number in (-1,1), got $p"))
end
Expand Down Expand Up @@ -1583,22 +1583,22 @@ function powm!(A0::UpperTriangular{<:BlasFloat}, p::Real)
end
copyto!(Stmp, S)
mul!(S, A, c)
ldiv!(Stmp, S.data)
ldiv!(Stmp, S)

c = (p - j) / (j4 - 2)
for i = 1:n
@inbounds S[i, i] = S[i, i] + 1
end
copyto!(Stmp, S)
mul!(S, A, c)
ldiv!(Stmp, S.data)
ldiv!(Stmp, S)
end
for i = 1:n
S[i, i] = S[i, i] + 1
end
copyto!(Stmp, S)
mul!(S, A, -p)
ldiv!(Stmp, S.data)
ldiv!(Stmp, S)
for i = 1:n
@inbounds S[i, i] = S[i, i] + 1
end
Expand Down
5 changes: 5 additions & 0 deletions stdlib/LinearAlgebra/test/dense.jl
Original file line number Diff line number Diff line change
Expand Up @@ -1067,6 +1067,11 @@ end
end
end

@testset "BigFloat triangular real power" begin
A = Float64[3 1; 0 3]
@test A^(3/4) big.(A)^(3/4)
end

@testset "diagonal integer matrix to real power" begin
A = Matrix(Diagonal([1, 2, 3]))
@test A^2.3 float(A)^2.3
Expand Down

0 comments on commit a2e4054

Please sign in to comment.