From 4706d1a7cf468c9082a04987b04477a4551b389a Mon Sep 17 00:00:00 2001 From: Sakse Dalum Date: Wed, 3 Oct 2018 20:32:38 +0800 Subject: [PATCH 1/5] Preserve types when adding/subtracting Herm/Sym/UniformScaling --- stdlib/LinearAlgebra/src/symmetric.jl | 15 +++++++++++++ stdlib/LinearAlgebra/src/uniformscaling.jl | 24 +++++++++++++++++++++ stdlib/LinearAlgebra/test/symmetric.jl | 14 +++++++++--- stdlib/LinearAlgebra/test/uniformscaling.jl | 10 +++++++++ 4 files changed, 60 insertions(+), 3 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index b5f6204361230..42ad93f272b2f 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -329,6 +329,13 @@ transpose(A::Hermitian{<:Real}) = A adjoint(A::Symmetric) = Adjoint(A) transpose(A::Hermitian) = Transpose(A) +real(A::Symmetric{<:Real}) = copy(A) +real(A::Hermitian{<:Real}) = copy(A) +real(A::Symmetric) = Symmetric(real(A.data), sym_uplo(A.uplo)) +real(A::Hermitian) = Hermitian(real(A.data), sym_uplo(A.uplo)) +imag(A::Symmetric) = Symmetric(imag(A.data), sym_uplo(A.uplo)) +imag(A::Hermitian) = Hermitian(imag(A.data), sym_uplo(A.uplo)) + Base.copy(A::Adjoint{<:Any,<:Hermitian}) = copy(A.parent) Base.copy(A::Transpose{<:Any,<:Symmetric}) = copy(A.parent) Base.copy(A::Adjoint{<:Any,<:Symmetric}) = @@ -393,6 +400,14 @@ end (-)(A::Symmetric{Tv,S}) where {Tv,S} = Symmetric{Tv,S}(-A.data, A.uplo) (-)(A::Hermitian{Tv,S}) where {Tv,S} = Hermitian{Tv,S}(-A.data, A.uplo) +## Addition/subtraction +for f in (:+, :-) + @eval $f(A::Symmetric, B::Symmetric) = Symmetric($f(A.data, B), sym_uplo(A.uplo)) + @eval $f(A::Hermitian, B::Hermitian) = Hermitian($f(A.data, B), sym_uplo(A.uplo)) + @eval $f(A::Hermitian, B::Symmetric{<:Real}) = Hermitian($f(A.data, B), sym_uplo(A.uplo)) + @eval $f(A::Symmetric{<:Real}, B::Hermitian) = Hermitian($f(A.data, B), sym_uplo(A.uplo)) +end + ## Matvec mul!(y::StridedVector{T}, A::Symmetric{T,<:StridedMatrix}, x::StridedVector{T}) where {T<:BlasFloat} = BLAS.symv!(A.uplo, one(T), A.data, x, zero(T), y) diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index 54db4d8a6b893..be23844bb5a43 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -110,6 +110,30 @@ for (t1, t2) in ((:UnitUpperTriangular, :UpperTriangular), end end +# Adding a complex UniformScaling to the diagonal of a Hermitian +# matrix breaks the hermiticity, if the UniformScaling is non-real. +# However, to preserve type stability, we do not special-case a +# UniformScaling{<:Complex} that happens to be real. +function (+)(A::Hermitian, J::UniformScaling{<:Complex}) + checksquare(A) + A_ = copytri!(copy(parent(A)), A.uplo) + B = convert(AbstractMatrix{Base._return_type(+, Tuple{eltype(A), typeof(J)})}, A_) + @inbounds for i in axes(A, 1) + B[i,i] += J + end + return B +end + +function (-)(J::UniformScaling{<:Complex}, A::Hermitian) + checksquare(A) + A_ = copytri!(copy(parent(A)), A.uplo) + B = convert(AbstractMatrix{Base._return_type(+, Tuple{eltype(A), typeof(J)})}, -A_) + @inbounds for i in axes(A, 1) + B[i,i] += J + end + return B +end + function (+)(A::AbstractMatrix, J::UniformScaling) checksquare(A) B = copy_oftype(A, Base._return_type(+, Tuple{eltype(A), typeof(J)})) diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 9e9814895d7e6..0a4aec9d4a212 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -90,6 +90,17 @@ end @test (-Hermitian(aherm))::typeof(Hermitian(aherm)) == -aherm end + @testset "Addition and subtraction for Symmetric/Hermitian matrices" begin + for f in (+, -) + @test (f(Symmetric(asym), Symmetric(aposs)))::typeof(Symmetric(asym)) == f(asym, aposs) + @test (f(Hermitian(aherm), Hermitian(apos)))::typeof(Hermitian(aherm)) == f(aherm, apos) + @test (f(Symmetric(real(asym)), Hermitian(aherm)))::typeof(Hermitian(aherm)) == f(real(asym), aherm) + @test (f(Hermitian(aherm), Symmetric(real(asym))))::typeof(Hermitian(aherm)) == f(aherm, real(asym)) + @test (f(Symmetric(asym), Hermitian(aherm))) == f(asym, aherm) + @test (f(Hermitian(aherm), Symmetric(asym))) == f(aherm, asym) + end + end + @testset "getindex and unsafe_getindex" begin @test aherm[1,1] == Hermitian(aherm)[1,1] @test asym[1,1] == Symmetric(asym)[1,1] @@ -415,9 +426,6 @@ end @test T([true false; false true]) .+ true == T([2 1; 1 2]) end - - @test_throws ArgumentError Hermitian(X) + 2im*I - @test_throws ArgumentError Hermitian(X) - 2im*I end @testset "Issue #21981" begin diff --git a/stdlib/LinearAlgebra/test/uniformscaling.jl b/stdlib/LinearAlgebra/test/uniformscaling.jl index 2e2b33a08250e..ac52bb271b073 100644 --- a/stdlib/LinearAlgebra/test/uniformscaling.jl +++ b/stdlib/LinearAlgebra/test/uniformscaling.jl @@ -173,6 +173,16 @@ let @test @inferred(J - T) == J - Array(T) @test @inferred(T\I) == inv(T) + if isa(A, Array) + T = Hermitian(randn(3,3)) + else + T = Hermitian(view(randn(3,3), 1:3, 1:3)) + end + @test @inferred(T + J) == Array(T) + J + @test @inferred(J + T) == J + Array(T) + @test @inferred(T - J) == Array(T) - J + @test @inferred(J - T) == J - Array(T) + @test @inferred(I\A) == A @test @inferred(A\I) == inv(A) @test @inferred(λ\I) === UniformScaling(1/λ) From f7241e40d24333820dfd77b37a734ad9d2e27ef8 Mon Sep 17 00:00:00 2001 From: Sakse Dalum Date: Wed, 3 Oct 2018 20:59:45 +0800 Subject: [PATCH 2/5] Make `real(::SymOrHerm{<:Real})` consistent with `real(::Array)`. --- stdlib/LinearAlgebra/src/symmetric.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 42ad93f272b2f..8d1162a5d67a7 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -329,8 +329,7 @@ transpose(A::Hermitian{<:Real}) = A adjoint(A::Symmetric) = Adjoint(A) transpose(A::Hermitian) = Transpose(A) -real(A::Symmetric{<:Real}) = copy(A) -real(A::Hermitian{<:Real}) = copy(A) +real(A::HermOrSym{<:Real}) = A real(A::Symmetric) = Symmetric(real(A.data), sym_uplo(A.uplo)) real(A::Hermitian) = Hermitian(real(A.data), sym_uplo(A.uplo)) imag(A::Symmetric) = Symmetric(imag(A.data), sym_uplo(A.uplo)) From 4b9b17cafa78815bb6349e8e317f9661531f4c10 Mon Sep 17 00:00:00 2001 From: Sakse Dalum Date: Wed, 3 Oct 2018 23:07:51 +0800 Subject: [PATCH 3/5] Fix embarrassing ambiguity --- stdlib/LinearAlgebra/src/symmetric.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index 8d1162a5d67a7..a49963f93260a 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -329,7 +329,8 @@ transpose(A::Hermitian{<:Real}) = A adjoint(A::Symmetric) = Adjoint(A) transpose(A::Hermitian) = Transpose(A) -real(A::HermOrSym{<:Real}) = A +real(A::Symmetric{<:Real}) = A +real(A::Hermitian{<:Real}) = A real(A::Symmetric) = Symmetric(real(A.data), sym_uplo(A.uplo)) real(A::Hermitian) = Hermitian(real(A.data), sym_uplo(A.uplo)) imag(A::Symmetric) = Symmetric(imag(A.data), sym_uplo(A.uplo)) From 7150bf27d3a581702d7a59f62c7fa189c235711d Mon Sep 17 00:00:00 2001 From: Sakse Dalum Date: Tue, 23 Oct 2018 16:59:54 +0800 Subject: [PATCH 4/5] More tests, remove imag(::Hermitian), simplify code --- stdlib/LinearAlgebra/src/symmetric.jl | 1 - stdlib/LinearAlgebra/src/uniformscaling.jl | 19 ++++++++++--------- stdlib/LinearAlgebra/test/symmetric.jl | 15 +++++++++++++++ 3 files changed, 25 insertions(+), 10 deletions(-) diff --git a/stdlib/LinearAlgebra/src/symmetric.jl b/stdlib/LinearAlgebra/src/symmetric.jl index a49963f93260a..60793d6b71a20 100644 --- a/stdlib/LinearAlgebra/src/symmetric.jl +++ b/stdlib/LinearAlgebra/src/symmetric.jl @@ -334,7 +334,6 @@ real(A::Hermitian{<:Real}) = A real(A::Symmetric) = Symmetric(real(A.data), sym_uplo(A.uplo)) real(A::Hermitian) = Hermitian(real(A.data), sym_uplo(A.uplo)) imag(A::Symmetric) = Symmetric(imag(A.data), sym_uplo(A.uplo)) -imag(A::Hermitian) = Hermitian(imag(A.data), sym_uplo(A.uplo)) Base.copy(A::Adjoint{<:Any,<:Hermitian}) = copy(A.parent) Base.copy(A::Transpose{<:Any,<:Symmetric}) = copy(A.parent) diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index be23844bb5a43..071de577b88e8 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -114,22 +114,23 @@ end # matrix breaks the hermiticity, if the UniformScaling is non-real. # However, to preserve type stability, we do not special-case a # UniformScaling{<:Complex} that happens to be real. -function (+)(A::Hermitian, J::UniformScaling{<:Complex}) - checksquare(A) +function (+)(A::Hermitian{T,S}, J::UniformScaling{<:Complex}) where {T,S} A_ = copytri!(copy(parent(A)), A.uplo) B = convert(AbstractMatrix{Base._return_type(+, Tuple{eltype(A), typeof(J)})}, A_) - @inbounds for i in axes(A, 1) - B[i,i] += J + @inbounds for i in diagind(B) + B[i] += J.λ end return B end -function (-)(J::UniformScaling{<:Complex}, A::Hermitian) - checksquare(A) +function (-)(J::UniformScaling{<:Complex}, A::Hermitian{T,S}) where {T,S} A_ = copytri!(copy(parent(A)), A.uplo) - B = convert(AbstractMatrix{Base._return_type(+, Tuple{eltype(A), typeof(J)})}, -A_) - @inbounds for i in axes(A, 1) - B[i,i] += J + B = convert(AbstractMatrix{Base._return_type(+, Tuple{eltype(A), typeof(J)})}, A_) + @inbounds for i in eachindex(B) + B[i] = -B[i] + end + @inbounds for i in diagind(B) + B[i] += J.λ end return B end diff --git a/stdlib/LinearAlgebra/test/symmetric.jl b/stdlib/LinearAlgebra/test/symmetric.jl index 0a4aec9d4a212..6604df1e6e165 100644 --- a/stdlib/LinearAlgebra/test/symmetric.jl +++ b/stdlib/LinearAlgebra/test/symmetric.jl @@ -164,6 +164,21 @@ end @test transpose(H) == Hermitian(copy(transpose(aherm))) end end + + @testset "real, imag" begin + S = Symmetric(asym) + H = Hermitian(aherm) + @test issymmetric(real(S)) + @test ishermitian(real(H)) + if eltya <: Real + @test real(S) === S == asym + @test real(H) === H == aherm + elseif eltya <: Complex + @test issymmetric(imag(S)) + @test !ishermitian(imag(H)) + end + end + end @testset "linalg unary ops" begin From f14426972b5fb68facbb86182c5f8be69678a2fb Mon Sep 17 00:00:00 2001 From: Sakse Dalum Date: Thu, 8 Nov 2018 14:32:10 +0800 Subject: [PATCH 5/5] =?UTF-8?q?Remove=20`.=CE=BB`?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit --- stdlib/LinearAlgebra/src/uniformscaling.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/uniformscaling.jl b/stdlib/LinearAlgebra/src/uniformscaling.jl index 071de577b88e8..ed05ff9998190 100644 --- a/stdlib/LinearAlgebra/src/uniformscaling.jl +++ b/stdlib/LinearAlgebra/src/uniformscaling.jl @@ -118,7 +118,7 @@ function (+)(A::Hermitian{T,S}, J::UniformScaling{<:Complex}) where {T,S} A_ = copytri!(copy(parent(A)), A.uplo) B = convert(AbstractMatrix{Base._return_type(+, Tuple{eltype(A), typeof(J)})}, A_) @inbounds for i in diagind(B) - B[i] += J.λ + B[i] += J end return B end @@ -130,7 +130,7 @@ function (-)(J::UniformScaling{<:Complex}, A::Hermitian{T,S}) where {T,S} B[i] = -B[i] end @inbounds for i in diagind(B) - B[i] += J.λ + B[i] += J end return B end