From d0f333ad9e607587f9ef15754294756954ad5741 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 27 Aug 2024 20:42:07 +0530 Subject: [PATCH 1/4] Reroute Symmetric/Hermitian + Diagonal through triangular --- stdlib/LinearAlgebra/src/diagonal.jl | 15 --------------- stdlib/LinearAlgebra/src/special.jl | 21 +++++++++++++++++++++ 2 files changed, 21 insertions(+), 15 deletions(-) diff --git a/stdlib/LinearAlgebra/src/diagonal.jl b/stdlib/LinearAlgebra/src/diagonal.jl index 92f399bb774ff..2f2d6097b1fb6 100644 --- a/stdlib/LinearAlgebra/src/diagonal.jl +++ b/stdlib/LinearAlgebra/src/diagonal.jl @@ -272,21 +272,6 @@ end (+)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag + Db.diag) (-)(Da::Diagonal, Db::Diagonal) = Diagonal(Da.diag - Db.diag) -for f in (:+, :-) - @eval function $f(D::Diagonal{<:Number}, S::Symmetric) - return Symmetric($f(D, S.data), sym_uplo(S.uplo)) - end - @eval function $f(S::Symmetric, D::Diagonal{<:Number}) - return Symmetric($f(S.data, D), sym_uplo(S.uplo)) - end - @eval function $f(D::Diagonal{<:Real}, H::Hermitian) - return Hermitian($f(D, H.data), sym_uplo(H.uplo)) - end - @eval function $f(H::Hermitian, D::Diagonal{<:Real}) - return Hermitian($f(H.data, D), sym_uplo(H.uplo)) - end -end - (*)(x::Number, D::Diagonal) = Diagonal(x * D.diag) (*)(D::Diagonal, x::Number) = Diagonal(D.diag * x) (/)(D::Diagonal, x::Number) = Diagonal(D.diag / x) diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl index 28a1b4b1b2eab..e352f016ef9b0 100644 --- a/stdlib/LinearAlgebra/src/special.jl +++ b/stdlib/LinearAlgebra/src/special.jl @@ -266,6 +266,27 @@ function (-)(A::UniformScaling, B::Diagonal) Diagonal(Ref(A) .- B.diag) end +for f in (:+, :-) + @eval function $f(D::Diagonal{<:Number}, S::Symmetric) + uplo = sym_uplo(S.uplo) + return Symmetric(parentof_applytri($f, Symmetric(D, uplo), S), uplo) + end + @eval function $f(S::Symmetric, D::Diagonal{<:Number}) + uplo = sym_uplo(S.uplo) + return Symmetric(parentof_applytri($f, S, Symmetric(D, uplo)), uplo) + end + @eval function $f(D::Diagonal{<:Real}, H::Hermitian) + uplo = sym_uplo(H.uplo) + return Hermitian(parentof_applytri($f, Hermitian(D, uplo), H), uplo) + end + @eval function $f(H::Hermitian, D::Diagonal{<:Real}) + uplo = sym_uplo(H.uplo) + return Hermitian(parentof_applytri($f, H, Hermitian(D, uplo)), uplo) + end + @eval $f(D::Diagonal, U::UpperOrLowerTriangular{<:Any, <:StridedMatrix}) = broadcast($f, D, U) + @eval $f(U::UpperOrLowerTriangular{<:Any, <:StridedMatrix}, D::Diagonal) = broadcast($f, U, D) +end + ## Diagonal construction from UniformScaling Diagonal{T}(s::UniformScaling, m::Integer) where {T} = Diagonal{T}(fill(T(s.λ), m)) Diagonal(s::UniformScaling, m::Integer) = Diagonal{eltype(s)}(s, m) From c032a9f11282f38e5590e6fadd069ec17475db70 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 28 Aug 2024 00:44:47 +0530 Subject: [PATCH 2/4] Add tests for partly filled parent --- stdlib/LinearAlgebra/test/special.jl | 13 +++++++++++++ 1 file changed, 13 insertions(+) diff --git a/stdlib/LinearAlgebra/test/special.jl b/stdlib/LinearAlgebra/test/special.jl index 9bb84ba0e9d03..f52a5fcf3ef12 100644 --- a/stdlib/LinearAlgebra/test/special.jl +++ b/stdlib/LinearAlgebra/test/special.jl @@ -786,4 +786,17 @@ end end end +@testset "Partly filled Hermitian and Diagonal algebra" begin + D = Diagonal([1,2]) + for S in (Symmetric, Hermitian), uplo in (:U, :L) + M = Matrix{BigInt}(undef, 2, 2) + M[1,1] = M[2,2] = M[1+(uplo == :L), 1 + (uplo == :U)] = 3 + H = S(M, uplo) + HM = Matrix(H) + @test H + D == D + H == HM + D + @test H - D == HM - D + @test D - H == D - HM + end +end + end # module TestSpecial From c70b517a8e940b9d262b52c3c6b8178fbfd3bbdf Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 28 Aug 2024 00:50:49 +0530 Subject: [PATCH 3/4] Restrict broadcasting to strided-diag Diagonal --- stdlib/LinearAlgebra/src/special.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl index e352f016ef9b0..57ac586c6b189 100644 --- a/stdlib/LinearAlgebra/src/special.jl +++ b/stdlib/LinearAlgebra/src/special.jl @@ -283,8 +283,8 @@ for f in (:+, :-) uplo = sym_uplo(H.uplo) return Hermitian(parentof_applytri($f, H, Hermitian(D, uplo)), uplo) end - @eval $f(D::Diagonal, U::UpperOrLowerTriangular{<:Any, <:StridedMatrix}) = broadcast($f, D, U) - @eval $f(U::UpperOrLowerTriangular{<:Any, <:StridedMatrix}, D::Diagonal) = broadcast($f, U, D) + @eval $f(D::Diagonal{<:Any, <:StridedVector}, U::UpperOrLowerTriangular{<:Any, <:StridedMatrix}) = broadcast($f, D, U) + @eval $f(U::UpperOrLowerTriangular{<:Any, <:StridedMatrix}, D::Diagonal{<:Any, <:StridedVector}) = broadcast($f, U, D) end ## Diagonal construction from UniformScaling From 4bd3de1181d7fe251714198c6d0ec33dd0382b4f Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Wed, 28 Aug 2024 00:53:31 +0530 Subject: [PATCH 4/4] Remove Diagonal-triangular specialization --- stdlib/LinearAlgebra/src/special.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/stdlib/LinearAlgebra/src/special.jl b/stdlib/LinearAlgebra/src/special.jl index 57ac586c6b189..3cf19db1542a6 100644 --- a/stdlib/LinearAlgebra/src/special.jl +++ b/stdlib/LinearAlgebra/src/special.jl @@ -283,8 +283,6 @@ for f in (:+, :-) uplo = sym_uplo(H.uplo) return Hermitian(parentof_applytri($f, H, Hermitian(D, uplo)), uplo) end - @eval $f(D::Diagonal{<:Any, <:StridedVector}, U::UpperOrLowerTriangular{<:Any, <:StridedMatrix}) = broadcast($f, D, U) - @eval $f(U::UpperOrLowerTriangular{<:Any, <:StridedMatrix}, D::Diagonal{<:Any, <:StridedVector}) = broadcast($f, U, D) end ## Diagonal construction from UniformScaling