diff --git a/NEWS.md b/NEWS.md index 788374aa38244..48807023238c1 100644 --- a/NEWS.md +++ b/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]) diff --git a/base/linalg/symmetric.jl b/base/linalg/symmetric.jl index f9f16d8ef7307..8c3c25728cd7d 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) @@ -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) diff --git a/base/sparse/cholmod.jl b/base/sparse/cholmod.jl index 4939f1af57890..e68417cad7a36 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 e8a79439724bd..5404a8ebb7bc6 100644 --- a/test/linalg/symmetric.jl +++ b/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 diff --git a/test/sparse/cholmod.jl b/test/sparse/cholmod.jl index 256732c71c308..f3aa5fe341049 100644 --- a/test/sparse/cholmod.jl +++ b/test/sparse/cholmod.jl @@ -692,6 +692,14 @@ end @test F\b ≈ ones(m + n) end +@testset "Test that imaginary parts in Hermitian{T,SparseMatrixCSC{T}} are ignored" begin + 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 + @testset "\\ '\\ and .'\\" begin # Test that \ and '\ and .'\ work for Symmetric and Hermitian. This is just # a dispatch exercise so it doesn't matter that the complex matrix has