Avoid checking the diagonal for non-real elements when constructing H…

…ermitians

NEWS.md

@@ -211,6 +211,9 @@ This section lists changes that do not have deprecation warnings.
     the type of `n`). Use the corresponding mutating functions `randperm!` and `randcycle!`
     to control the array type ([#22723]).
 
+  * Hermitian now ignores any imaginary components in the diagonal instead of checking
+    the diagonal. ([#17367])
+
   * Worker-worker connections are setup lazily for an `:all_to_all` topology. Use keyword
     arg `lazy=false` to force all connections to be setup during a `addprocs` call. ([#22814])
 

base/linalg/symmetric.jl

@@ -74,13 +74,15 @@ julia> Hlower = Hermitian(A, :L)
 ```
 
 Note that `Hupper` will not be equal to `Hlower` unless `A` is itself Hermitian (e.g. if `A == A'`).
+
+All non-real parts of the diagonal will be ignored.
+
+```julia
+Hermitian(fill(complex(1,1), 1, 1)) == fill(1, 1, 1)
+```
 """
 function Hermitian(A::AbstractMatrix, uplo::Symbol=:U)
     n = checksquare(A)
-    for i=1:n
-        isreal(A[i, i]) || throw(ArgumentError(
-            "Cannot construct Hermitian from matrix with nonreal diagonals"))
-    end
     Hermitian{eltype(A),typeof(A)}(A, char_uplo(uplo))
 end
 
@@ -113,13 +115,21 @@ size(A::HermOrSym, d) = size(A.data, d)
 size(A::HermOrSym) = size(A.data)
 @inline function getindex(A::Symmetric, i::Integer, j::Integer)
     @boundscheck checkbounds(A, i, j)
-    @inbounds r = (A.uplo == 'U') == (i < j) ? A.data[i, j] : A.data[j, i]
-    r
+    @inbounds if (A.uplo == 'U') == (i < j)
+        return A.data[i, j]
+    else
+        return A.data[j, i]
+    end
 end
 @inline function getindex(A::Hermitian, i::Integer, j::Integer)
     @boundscheck checkbounds(A, i, j)
-    @inbounds r = (A.uplo == 'U') == (i < j) ? A.data[i, j] : conj(A.data[j, i])
-    r
+    @inbounds if (A.uplo == 'U') == (i < j)
+        return A.data[i, j]
+    elseif i == j
+        return eltype(A)(real(A.data[i, j]))
+    else
+        return conj(A.data[j, i])
+    end
 end
 
 function setindex!(A::Symmetric, v, i::Integer, j::Integer)
@@ -153,8 +163,15 @@ similar(A::Union{Symmetric,Hermitian}, ::Type{T}, dims::Dims{N}) where {T,N} = s
 
 # Conversion
 convert(::Type{Matrix}, A::Symmetric) = copytri!(convert(Matrix, copy(A.data)), A.uplo)
-convert(::Type{Matrix}, A::Hermitian) = copytri!(convert(Matrix, copy(A.data)), A.uplo, true)
+function convert(::Type{Matrix}, A::Hermitian)
+    B = copytri!(convert(Matrix, copy(A.data)), A.uplo, true)
+    for i = 1:size(A, 1)
+        B[i,i] = real(B[i,i])
+    end
+    return B
+end
 convert(::Type{Array}, A::Union{Symmetric,Hermitian}) = convert(Matrix, A)
+
 parent(A::HermOrSym) = A.data
 convert(::Type{Symmetric{T,S}},A::Symmetric{T,S}) where {T,S<:AbstractMatrix} = A
 convert(::Type{Symmetric{T,S}},A::Symmetric) where {T,S<:AbstractMatrix} = Symmetric{T,S}(convert(S,A.data),A.uplo)

base/sparse/cholmod.jl

@@ -928,8 +928,34 @@ convert(::Type{Sparse}, A::SparseMatrixCSC{Complex{Float32},<:ITypes}) =
     convert(Sparse, convert(SparseMatrixCSC{Complex{Float64},SuiteSparse_long}, A))
 convert(::Type{Sparse}, A::Symmetric{Float64,SparseMatrixCSC{Float64,SuiteSparse_long}}) =
     Sparse(A.data, A.uplo == 'L' ? -1 : 1)
-convert(::Type{Sparse}, A::Hermitian{Tv,SparseMatrixCSC{Tv,SuiteSparse_long}}) where {Tv<:VTypes} =
-    Sparse(A.data, A.uplo == 'L' ? -1 : 1)
+# The convert method for Hermitian is very similar to the general convert method, but we need to
+# remove any non real elements in the diagonal because, in contrast to BLAS/LAPACK these are not
+# ignored by CHOLMOD. If even tiny imaginary parts are present CHOLMOD will fail with a non-positive
+# definite/zero pivot error.
+function convert(::Type{Sparse}, AH::Hermitian{Tv,SparseMatrixCSC{Tv,SuiteSparse_long}}) where {Tv<:VTypes}
+    A = AH.data
+
+    # Here we allocate a Symmetric/Hermitian CHOLMOD.Sparse matrix so we only need to copy
+    # a single triangle of AH
+    o = allocate_sparse(A.m, A.n, length(A.nzval), true, true, AH.uplo == 'L' ? -1 : 1, Tv)
+    s = unsafe_load(o.p)
+    for i = 1:length(A.colptr)
+        unsafe_store!(s.p, A.colptr[i] - 1, i)
+    end
+    for i = 1:length(A.rowval)
+        unsafe_store!(s.i, A.rowval[i] - 1, i)
+    end
+    for j = 1:A.n
+        for ip = A.colptr[j]:A.colptr[j + 1] - 1
+            v = A.nzval[ip]
+            unsafe_store!(s.x, ifelse(A.rowval[ip] == j, real(v), v), ip)
+        end
+    end
+
+    @isok check_sparse(o)
+
+    return o
+end
 function convert(::Type{Sparse},
                  A::Union{SparseMatrixCSC{BigFloat,Ti},
                           Symmetric{BigFloat,SparseMatrixCSC{BigFloat,Ti}},

test/linalg/symmetric.jl

@@ -349,6 +349,12 @@ end
     @test eigvals(Mi7647) == eigvals(Mi7647, 0.5, 3.5) == [1.0:3.0;]
 end
 
+@testset "Hermitian wrapper ignores imaginary parts on diagonal" begin
+    A = [1.0+im 2.0; 2.0 0.0]
+    @test !ishermitian(A)
+    @test Hermitian(A)[1,1] == 1
+end
+
 @testset "Issue #7933" begin
     A7933 = [1 2; 3 4]
     B7933 = copy(A7933)
@@ -363,10 +369,10 @@ end
     @test_throws ArgumentError f(A, 1:4)
 end
 
-@testset "Issue #10671" begin
+@testset "Ignore imaginary part of Hermitian diagonal" begin
     A = [1.0+im 2.0; 2.0 0.0]
     @test !ishermitian(A)
-    @test_throws ArgumentError Hermitian(A)
+    @test diag(Hermitian(A)) == real(diag(A))
 end
 
 @testset "Issue #17780" begin

test/sparse/cholmod.jl

@@ -654,6 +654,14 @@ let A = sprandn(5, 5, 0.4) |> t -> t't + I,
     @test cholfact(A)\B ≈ A\B
 end
 
+# Test that imaginary parts in Hermitian{T,SparseMatrixCSC{T}} are ignored
+let A = sparse([1,2,3,4,1], [1,2,3,4,2], [complex(2.0,1),2,2,2,1])
+    Fs = cholfact(Hermitian(A))
+    Fd = cholfact(Hermitian(Array(A)))
+    @test sparse(Fs) ≈ Hermitian(A)
+    @test Fs\ones(4) ≈ Fd\ones(4)
+end
+
 # Make sure that ldltfact performs an LDLt (Issue #19032)
 let m = 400, n = 500,
     A = sprandn(m, n, .2),