From 1c25c497b53090eda04befef7c61abcef2d782c9 Mon Sep 17 00:00:00 2001
From: Daniel Karrasch <daniel.karrasch@posteo.de>
Date: Sat, 27 Apr 2024 13:52:44 +0200
Subject: [PATCH] Fix solve for complex `Hermitian`  with non-vanishing
 imaginary part on diagonal (#54276)

(cherry picked from commit e80a88061cfb2891f3bc051005dabdd12c09bb8f)
---
 stdlib/LinearAlgebra/src/lu.jl         | 5 +++++
 stdlib/LinearAlgebra/test/symmetric.jl | 7 +++++--
 2 files changed, 10 insertions(+), 2 deletions(-)

diff --git a/stdlib/LinearAlgebra/src/lu.jl b/stdlib/LinearAlgebra/src/lu.jl
index 72755e0eb6799..d2e82af5d6409 100644
--- 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
index d3b24ccf78b0b..4d9fce18cf5a6 100644
--- 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