diff --git a/stdlib/LinearAlgebra/src/cholesky.jl b/stdlib/LinearAlgebra/src/cholesky.jl index 60f171838c567..10aabb35928b8 100644 --- a/stdlib/LinearAlgebra/src/cholesky.jl +++ b/stdlib/LinearAlgebra/src/cholesky.jl @@ -181,12 +181,15 @@ end function _chol!(A::AbstractMatrix, ::Type{UpperTriangular}) require_one_based_indexing(A) n = checksquare(A) + realdiag = eltype(A) <: Complex @inbounds begin for k = 1:n + Akk = realdiag ? real(A[k,k]) : A[k,k] for i = 1:k - 1 - A[k,k] -= A[i,k]'A[i,k] + Akk -= realdiag ? abs2(A[i,k]) : A[i,k]'A[i,k] end - Akk, info = _chol!(A[k,k], UpperTriangular) + A[k,k] = Akk + Akk, info = _chol!(Akk, UpperTriangular) if info != 0 return UpperTriangular(A), info end @@ -205,12 +208,15 @@ end function _chol!(A::AbstractMatrix, ::Type{LowerTriangular}) require_one_based_indexing(A) n = checksquare(A) + realdiag = eltype(A) <: Complex @inbounds begin for k = 1:n + Akk = realdiag ? real(A[k,k]) : A[k,k] for i = 1:k - 1 - A[k,k] -= A[k,i]*A[k,i]' + Akk -= realdiag ? abs2(A[k,i]) : A[k,i]*A[k,i]' end - Akk, info = _chol!(A[k,k], LowerTriangular) + A[k,k] = Akk + Akk, info = _chol!(Akk, LowerTriangular) if info != 0 return LowerTriangular(A), info end diff --git a/stdlib/LinearAlgebra/test/cholesky.jl b/stdlib/LinearAlgebra/test/cholesky.jl index 888ff065b5839..86a78a9e954fe 100644 --- a/stdlib/LinearAlgebra/test/cholesky.jl +++ b/stdlib/LinearAlgebra/test/cholesky.jl @@ -440,4 +440,11 @@ end end end +@testset "issue #37356, diagonal elements of hermitian generic matrix" begin + B = Hermitian(hcat([one(BigFloat) + im])) + @test Matrix(cholesky(B)) ≈ B + C = Hermitian(hcat([one(BigFloat) + im]), :L) + @test Matrix(cholesky(C)) ≈ C +end + end # module TestCholesky