Fix solve for complex Hermitian with non-vanishing imaginary part on diagonal (#54276)

(cherry picked from commit e80a880)

Author: dkarrasch

stdlib/LinearAlgebra/src/lu.jl

@@ -85,6 +85,11 @@ function lu!(A::StridedMatrix{T}, ::RowMaximum; check::Bool = true) where {T<:Bl
 end
 function lu!(A::HermOrSym{T}, pivot::Union{RowMaximum,NoPivot,RowNonZero} = lupivottype(T); check::Bool = true) 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 = check)
 end
 # for backward compatibility

stdlib/LinearAlgebra/test/symmetric.jl

@@ -261,8 +261,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