use real diagonal in cholesky (#37382)

JuliaLang · Sep 12, 2020 · eeffe57 · eeffe57
1 parent 7ecd1f9
commit eeffe57
Show file tree

Hide file tree

Showing 2 changed files with 17 additions and 4 deletions.
diff --git 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
@@ -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