From 4a5cbbb1637f56eaf8f6d158a0cf1d2177842379 Mon Sep 17 00:00:00 2001 From: Andreas Noack 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 | 35 ++++++++++++++++++++++++++--------- base/sparse/cholmod.jl | 30 ++++++++++++++++++++++++++++-- test/linalg/symmetric.jl | 10 ++++++++-- test/sparse/cholmod.jl | 8 ++++++++ 4 files changed, 70 insertions(+), 13 deletions(-) diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index 1a890da1a5007..383eff2757804 100644 --- a/base/linalg/symmetric.jl +++ b/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) @@ -150,9 +160,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!(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) full(A::Union{Symmetric,Hermitian}) = convert(Array, 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) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index b29f4ed701103..8ce25f4c2eb09 100644 --- a/base/sparse/cholmod.jl +++ b/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}}, diff --git a/test/linalg/symmetric.jl b/test/linalg/symmetric.jl index 296d599390d89..0042d6d692d6f 100644 --- a/test/linalg/symmetric.jl +++ b/test/linalg/symmetric.jl @@ -325,6 +325,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) @@ -339,10 +345,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 # Unary minus for Symmetric/Hermitian matrices diff --git a/test/sparse/cholmod.jl b/test/sparse/cholmod.jl index 936f6059e6a37..66070931936da 100644 --- a/test/sparse/cholmod.jl +++ b/test/sparse/cholmod.jl @@ -650,6 +650,14 @@ A = sprandn(5,5,0.4) |> t -> t't + I B = complex.(randn(5,2), randn(5,2)) @test cholfact(A)\B ≈ A\B +# 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)