use symmetrized lufact in inv of Symmetric/Hermitian matrices (JuliaL…

…ang#22686)
jeffwong · Jul 24, 2017 · a3498cc · a3498cc
1 parent dfff018
commit a3498cc
Show file tree

Hide file tree

Showing 2 changed files with 45 additions and 10 deletions.
diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl
@@ -343,9 +343,24 @@ det(A::Symmetric) = det(bkfact(A))
 # ambiguity with RowVector
 \(A::HermOrSym{<:Any,<:StridedMatrix}, B::RowVector) = invoke(\, Tuple{AbstractMatrix, RowVector}, A, B)
 
-
-inv(A::Hermitian{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = Hermitian{T,S}(inv(bkfact(A)), A.uplo)
-inv(A::Symmetric{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = Symmetric{T,S}(inv(bkfact(A)), A.uplo)
+function _inv(A::HermOrSym)
+    n = checksquare(A)
+    B = inv!(lufact(A))
+    conjugate = isa(A, Hermitian)
+    # symmetrize
+    if A.uplo == 'U' # add to upper triangle
+        @inbounds for i = 1:n, j = i:n
+            B[i,j] = conjugate ? (B[i,j] + conj(B[j,i])) / 2 : (B[i,j] + B[j,i]) / 2
+        end
+    else # A.uplo == 'L', add to lower triangle
+        @inbounds for i = 1:n, j = i:n
+            B[j,i] = conjugate ? (B[j,i] + conj(B[i,j])) / 2 : (B[j,i] + B[i,j]) / 2
+        end
+    end
+    B
+end
+inv(A::Hermitian{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = Hermitian{T,S}(_inv(A), A.uplo)
+inv(A::Symmetric{T,S}) where {T<:BlasFloat,S<:StridedMatrix} = Symmetric{T,S}(_inv(A), A.uplo)
 
 eigfact!(A::RealHermSymComplexHerm{<:BlasReal,<:StridedMatrix}) = Eigen(LAPACK.syevr!('V', 'A', A.uplo, A.data, 0.0, 0.0, 0, 0, -1.0)...)
 

diff --git a/test/linalg/symmetric.jl b/test/linalg/symmetric.jl
@@ -203,12 +203,14 @@ end
             @test Symmetric(asym)\x ≈ asym\x
         end
 
-        #inversion
-        @test inv(Hermitian(asym)) ≈ inv(asym)
-        if eltya <: Real && eltya != Int
-            @test inv(Symmetric(asym)) ≈ inv(asym)
-            @test inv(Hermitian(a)) ≈ inv(full(Hermitian(a)))
-            @test inv(Symmetric(a)) ≈ inv(full(Symmetric(a)))
+        # inversion
+        for uplo in (:U, :L)
+            @test inv(Hermitian(asym, uplo)) ≈ inv(full(Hermitian(asym, uplo)))
+            @test inv(Symmetric(asym, uplo)) ≈ inv(full(Symmetric(asym, uplo)))
+            if eltya <: Real
+                @test inv(Hermitian(a, uplo)) ≈ inv(full(Hermitian(a, uplo)))
+            end
+            @test inv(Symmetric(a, uplo)) ≈ inv(full(Symmetric(a, uplo)))
         end
 
         # conversion
@@ -344,11 +346,29 @@ end
     @test issymmetric(Hermitian(B, :L))
 end
 
-@testset "$HS solver with $RHS RHS - $T" for HS  in (Hermitian, Symmetric),
+@testset "$HS solver with $RHS RHS - $T" for HS in (Hermitian, Symmetric),
         RHS in (Hermitian, Symmetric, Diagonal, UpperTriangular, LowerTriangular),
         T   in (Float64, Complex128)
     D = rand(T, 10, 10); D = D'D
     A = HS(D)
     B = RHS(D)
     @test A\B ≈ Matrix(A)\Matrix(B)
 end
+
+@testset "inversion of Hilbert matrix" begin
+    for T in (Float64, Complex128)
+        H = T[1/(i + j - 1) for i in 1:8, j in 1:8]
+        @test norm(inv(Symmetric(H))*(H*ones(8))-1) ≈ 0 atol = 1e-5
+        @test norm(inv(Hermitian(H))*(H*ones(8))-1) ≈ 0 atol = 1e-5
+    end
+end
+
+@testset "inverse edge case with complex Hermitian" begin
+    # Hermitian matrix, where inv(lufact(A)) generates non-real diagonal elements
+    for T in (Complex64, Complex128)
+        A = T[0.650488+0.0im 0.826686+0.667447im; 0.826686-0.667447im 1.81707+0.0im]
+        H = Hermitian(A)
+        @test inv(H) ≈ inv(A)
+        @test ishermitian(full(inv(H)))
+    end
+end