From d1eb60925950daf6106eb6e8459991b21a252e35 Mon Sep 17 00:00:00 2001
From: Andreas Noack <andreasnoackjensen@gmail.com>
Date: Sun, 10 Jul 2016 22:13:00 -0400
Subject: [PATCH] Avoid checking the diagonal for non-real elements when
 constructing Hermitians

---
 base/linalg/symmetric.jl  | 84 +++++++++++++--------------------------
 base/sparse/cholmod.jl    | 24 ++++++++++-
 test/linalg/symmetric.jl  | 12 +++---
 test/sparsedir/cholmod.jl | 15 +++++--
 4 files changed, 69 insertions(+), 66 deletions(-)

diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl
index 0a52b6041adbf9..1eba7620d4bd2f 100644
--- a/base/linalg/symmetric.jl
+++ b/base/linalg/symmetric.jl
@@ -43,13 +43,15 @@ eigfact(Hupper)
 
 `eigfact` will use a method specialized for matrices known to be Hermitian.
 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
 
@@ -60,13 +62,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
 
 similar{T}(A::Symmetric, ::Type{T}) = Symmetric(similar(A.data, T))
@@ -82,9 +92,16 @@ end
 
 # 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!(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)
 full(A::Union{Symmetric,Hermitian}) = convert(Array, A)
+
 parent(A::HermOrSym) = A.data
 convert{T,S<:AbstractMatrix}(::Type{Symmetric{T,S}},A::Symmetric{T,S}) = A
 convert{T,S<:AbstractMatrix}(::Type{Symmetric{T,S}},A::Symmetric) = Symmetric{T,S}(convert(S,A.data),A.uplo)
@@ -115,53 +132,8 @@ ctranspose(A::Hermitian) = A
 trace(A::Hermitian) = real(trace(A.data))
 
 #tril/triu
-function tril(A::Hermitian, k::Integer=0)
-    if A.uplo == 'U' && k <= 0
-        return tril!(A.data',k)
-    elseif A.uplo == 'U' && k > 0
-        return tril!(A.data',-1) + tril!(triu(A.data),k)
-    elseif A.uplo == 'L' && k <= 0
-        return tril(A.data,k)
-    else
-        return tril(A.data,-1) + tril!(triu!(A.data'),k)
-    end
-end
-
-function tril(A::Symmetric, k::Integer=0)
-    if A.uplo == 'U' && k <= 0
-        return tril!(A.data.',k)
-    elseif A.uplo == 'U' && k > 0
-        return tril!(A.data.',-1) + tril!(triu(A.data),k)
-    elseif A.uplo == 'L' && k <= 0
-        return tril(A.data,k)
-    else
-        return tril(A.data,-1) + tril!(triu!(A.data.'),k)
-    end
-end
-
-function triu(A::Hermitian, k::Integer=0)
-    if A.uplo == 'U' && k >= 0
-        return triu(A.data,k)
-    elseif A.uplo == 'U' && k < 0
-        return triu(A.data,1) + triu!(tril!(A.data'),k)
-    elseif A.uplo == 'L' && k >= 0
-        return triu!(A.data',k)
-    else
-        return triu!(A.data',1) + triu!(tril(A.data),k)
-    end
-end
-
-function triu(A::Symmetric, k::Integer=0)
-    if A.uplo == 'U' && k >= 0
-        return triu(A.data,k)
-    elseif A.uplo == 'U' && k < 0
-        return triu(A.data,1) + triu!(tril!(A.data.'),k)
-    elseif A.uplo == 'L' && k >= 0
-        return triu!(A.data.',k)
-    else
-        return triu!(A.data.',1) + triu!(tril(A.data),k)
-    end
-end
+tril(A::HermOrSym, k::Integer = 0) = tril!(full(A), k)
+triu(A::HermOrSym, k::Integer = 0) = triu!(full(A), k)
 
 -{Tv,S<:AbstractMatrix}(A::Symmetric{Tv,S}) = Symmetric{Tv,S}(-A.data, A.uplo)
 
diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl
index 5dd273ceb7e959..6f9759dd7033bf 100644
--- a/base/sparse/cholmod.jl
+++ b/base/sparse/cholmod.jl
@@ -886,7 +886,29 @@ function convert{Tv<:VTypes,Ti<:ITypes}(::Type{Sparse}, A::SparseMatrixCSC{Tv,Ti
 end
 convert{Ti<:ITypes}(::Type{Sparse}, A::SparseMatrixCSC{Complex{Float32},Ti}) = 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{Tv<:VTypes}(::Type{Sparse}, A::Hermitian{Tv,SparseMatrixCSC{Tv,SuiteSparse_long}}) = Sparse(A.data, A.uplo == 'L' ? -1 : 1)
+function convert{Tv<:VTypes}(::Type{Sparse}, AH::Hermitian{Tv,SparseMatrixCSC{Tv,SuiteSparse_long}})
+    A = AH.data
+
+    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{Ti<:ITypes}(::Type{Sparse},
     A::Union{SparseMatrixCSC{BigFloat,Ti},
              Symmetric{BigFloat,SparseMatrixCSC{BigFloat,Ti}},
diff --git a/test/linalg/symmetric.jl b/test/linalg/symmetric.jl
index 3a33105b5372ff..8f95a6a1916404 100644
--- a/test/linalg/symmetric.jl
+++ b/test/linalg/symmetric.jl
@@ -199,6 +199,12 @@ let n=10
     end
 end
 
+# Hermitian wrapper ignores imaginary parts
+let A = [1.0+im 2.0; 2.0 0.0]
+    @test !ishermitian(A)
+    @test Hermitian(A)[1,1] == 1
+end
+
 #Issue #7647: test xsyevr, xheevr, xstevr drivers
 for Mi7647 in (Symmetric(diagm(1.0:3.0)),
                Hermitian(diagm(1.0:3.0)),
@@ -225,12 +231,6 @@ for f in (eigfact, eigvals, eig)
     end
 end
 
-#Issue 10671
-let A = [1.0+im 2.0; 2.0 0.0]
-    @test !ishermitian(A)
-    @test_throws ArgumentError Hermitian(A)
-end
-
 # Unary minus for Symmetric matrices
 let A = Symmetric(randn(5,5))
     B = -A
diff --git a/test/sparsedir/cholmod.jl b/test/sparsedir/cholmod.jl
index f7b4461911a70f..7cc548f378b362 100644
--- a/test/sparsedir/cholmod.jl
+++ b/test/sparsedir/cholmod.jl
@@ -649,6 +649,15 @@ let x = rand(5)
 end
 
 # Real factorization and complex rhs
-A = sprandn(5,5,0.4) |> t -> t't + I
-B = complex(randn(5,2), randn(5,2))
-@test cholfact(A)\B ≈ A\B
+let A = sprandn(5,5,0.4) |> t -> t't + I, B = complex(randn(5,2), randn(5,2))
+    @test cholfact(A)\B ≈ A\B
+end
+
+# Test that imaginary parts in Hermitian{T,SparseMatrixCSC{T}} is ignored
+let A = sparse([1,2,3,4,1], [1,2,3,4,2], [2.0,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
+