Fix solve for complex Hermitian with non-vanishing imaginary part o…

…n diagonal (#54276)
JuliaLang · Apr 27, 2024 · e80a880 · e80a880
1 parent 83929ad
commit e80a880
Show file tree

Hide file tree

Showing 2 changed files with 10 additions and 2 deletions.
diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl
@@ -95,6 +95,11 @@ end
 function lu!(A::HermOrSym{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T);
         check::Bool = true, allowsingular::Bool = false) where {T}
     copytri!(A.data, A.uplo, isa(A, Hermitian))
+    @inbounds if isa(A, Hermitian) # realify diagonal
+        for i in axes(A, 1)
+            A.data[i,i] = A[i,i]
+        end
+    end
     lu!(A.data, pivot; check, allowsingular)
 end
 # for backward compatibility

diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl
@@ -264,8 +264,11 @@ end
                     @testset "inverse edge case with complex Hermitian" begin
                         # Hermitian matrix, where inv(lu(A)) generates non-real diagonal elements
                         for T in (ComplexF32, ComplexF64)
-                            A = T[0.650488+0.0im 0.826686+0.667447im; 0.826686-0.667447im 1.81707+0.0im]
-                            H = Hermitian(A)
+                            # data should have nonvanishing imaginary parts on the diagonal
+                            M = T[0.279982+0.988074im  0.770011+0.870555im
+                                    0.138001+0.889728im  0.177242+0.701413im]
+                            H = Hermitian(M)
+                            A = Matrix(H)
                             @test inv(H) ≈ inv(A)
                             @test ishermitian(Matrix(inv(H)))
                         end