Reland "Reroute Symmetric/Hermitian + Diagonal through triangular"

This backports the following commits:
commit 9690961c426ce2640d7db6c89952e69f87873a93
Author: Jishnu Bhattacharya <[email protected]>
Date:   Mon Apr 29 21:43:31 2024 +0530

    Add upper/lowertriangular functions and use in applytri (#53573)

    We may use the fact that a `Diagonal` is already triangular to avoid
    adding a wrapper.

    Fixes the specific example in
    #53564, although not the
    broader issue. This is because it changes the operation from a
    `UpperTriangular + UpperTriangular` to a `UpperTriangular + Diagonal`,
    which uses broadcasting. The latter operation may also allow one to
    define more efficient methods.

commit 77821cdddb968eeabf31ccb6b214ccf59a604c68
Author: Jishnu Bhattacharya <[email protected]>
Date:   Wed Aug 28 00:53:31 2024 +0530

    Remove Diagonal-triangular specialization

commit 621fb2e739a04207df63857700aca3562b41b5eb
Author: Jishnu Bhattacharya <[email protected]>
Date:   Wed Aug 28 00:50:49 2024 +0530

    Restrict broadcasting to strided-diag Diagonal

commit 58eb2045ddb5dbbfdb759c06239ca54751e73d71
Author: Jishnu Bhattacharya <[email protected]>
Date:   Wed Aug 28 00:44:47 2024 +0530

    Add tests for partly filled parent

commit 5aa6080a580bfbc9453e94a06f3e379e4517b316
Author: Jishnu Bhattacharya <[email protected]>
Date:   Tue Aug 27 20:42:07 2024 +0530

    Reroute Symmetric/Hermitian + Diagonal through triangular

stdlib/LinearAlgebra/src/diagonal.jl

@@ -250,21 +250,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)
@@ -991,3 +976,6 @@ end
 function Base.muladd(A::Diagonal, B::Diagonal, z::Diagonal)
     Diagonal(A.diag .* B.diag .+ z.diag)
 end
+
+uppertriangular(D::Diagonal) = D
+lowertriangular(D::Diagonal) = D

stdlib/LinearAlgebra/src/special.jl

@@ -264,6 +264,25 @@ 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
+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)

stdlib/LinearAlgebra/src/symmetric.jl

@@ -277,21 +277,21 @@ diag(A::Hermitian) = hermitian.(diag(parent(A)), sym_uplo(A.uplo))
 
 function applytri(f, A::HermOrSym)
     if A.uplo == 'U'
-        f(UpperTriangular(A.data))
+        f(uppertriangular(A.data))
     else
-        f(LowerTriangular(A.data))
+        f(lowertriangular(A.data))
     end
 end
 
 function applytri(f, A::HermOrSym, B::HermOrSym)
     if A.uplo == B.uplo == 'U'
-        f(UpperTriangular(A.data), UpperTriangular(B.data))
+        f(uppertriangular(A.data), uppertriangular(B.data))
     elseif A.uplo == B.uplo == 'L'
-        f(LowerTriangular(A.data), LowerTriangular(B.data))
+        f(lowertriangular(A.data), lowertriangular(B.data))
     elseif A.uplo == 'U'
-        f(UpperTriangular(A.data), UpperTriangular(_conjugation(B)(B.data)))
+        f(uppertriangular(A.data), uppertriangular(_conjugation(B)(B.data)))
     else # A.uplo == 'L'
-        f(UpperTriangular(_conjugation(A)(A.data)), UpperTriangular(B.data))
+        f(uppertriangular(_conjugation(A)(A.data)), uppertriangular(B.data))
     end
 end
 parentof_applytri(f, args...) = applytri(parent ∘ f, args...)

stdlib/LinearAlgebra/src/triangular.jl

@@ -153,6 +153,12 @@ const UpperOrUnitUpperTriangular{T,S} = Union{UpperTriangular{T,S}, UnitUpperTri
 const LowerOrUnitLowerTriangular{T,S} = Union{LowerTriangular{T,S}, UnitLowerTriangular{T,S}}
 const UpperOrLowerTriangular{T,S} = Union{UpperOrUnitUpperTriangular{T,S}, LowerOrUnitLowerTriangular{T,S}}
 
+uppertriangular(M) = UpperTriangular(M)
+lowertriangular(M) = LowerTriangular(M)
+
+uppertriangular(U::UpperOrUnitUpperTriangular) = U
+lowertriangular(U::LowerOrUnitLowerTriangular) = U
+
 imag(A::UpperTriangular) = UpperTriangular(imag(A.data))
 imag(A::LowerTriangular) = LowerTriangular(imag(A.data))
 imag(A::UpperTriangular{<:Any,<:StridedMaybeAdjOrTransMat}) = imag.(A)

stdlib/LinearAlgebra/test/diagonal.jl

@@ -1277,6 +1277,12 @@ end
     @test c == Diagonal([2,2,2,2])
 end
 
+@testset "uppertriangular/lowertriangular" begin
+    D = Diagonal([1,2])
+    @test LinearAlgebra.uppertriangular(D) === D
+    @test LinearAlgebra.lowertriangular(D) === D
+end
+
 @testset "mul/div with an adjoint vector" begin
     A = [1.0;;]
     x = [1.0]

stdlib/LinearAlgebra/test/special.jl

@@ -536,4 +536,17 @@ end
     @test v * S isa Matrix
 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

stdlib/LinearAlgebra/test/symmetric.jl

@@ -507,6 +507,31 @@ end
             @test Su - Sl == -(Sl - Su) == MSu - MSl
         end
     end
+    @testset "non-strided" begin
+        @testset "diagonal" begin
+            for ST1 in (Symmetric, Hermitian), uplo1 in (:L, :U)
+                m = ST1(Matrix{BigFloat}(undef,2,2), uplo1)
+                m.data[1,1] = 1
+                m.data[2,2] = 3
+                m.data[1+(uplo1==:L), 1+(uplo1==:U)] = 2
+                A = Array(m)
+                for ST2 in (Symmetric, Hermitian), uplo2 in (:L, :U)
+                    id = ST2(I(2), uplo2)
+                    @test m + id == id + m == A + id
+                end
+            end
+        end
+        @testset "unit triangular" begin
+            for ST1 in (Symmetric, Hermitian), uplo1 in (:L, :U)
+                H1 = ST1(UnitUpperTriangular(big.(rand(Int8,4,4))), uplo1)
+                M1 = Matrix(H1)
+                for ST2 in (Symmetric, Hermitian), uplo2 in (:L, :U)
+                    H2 = ST2(UnitUpperTriangular(big.(rand(Int8,4,4))), uplo2)
+                    @test H1 + H2 == M1 + Matrix(H2)
+                end
+            end
+        end
+    end
 end
 
 # bug identified in PR #52318: dot products of quaternionic Hermitian matrices,

stdlib/LinearAlgebra/test/triangular.jl

@@ -996,6 +996,14 @@ end
     end
 end
 
+@testset "uppertriangular/lowertriangular" begin
+    M = rand(2,2)
+    @test LinearAlgebra.uppertriangular(M) === UpperTriangular(M)
+    @test LinearAlgebra.lowertriangular(M) === LowerTriangular(M)
+    @test LinearAlgebra.uppertriangular(UnitUpperTriangular(M)) === UnitUpperTriangular(M)
+    @test LinearAlgebra.lowertriangular(UnitLowerTriangular(M)) === UnitLowerTriangular(M)
+end
+
 @testset "arithmetic with partly uninitialized matrices" begin
     @testset "$(typeof(A))" for A in (Matrix{BigFloat}(undef,2,2), Matrix{Complex{BigFloat}}(undef,2,2)')
         A[2,1] = eltype(A) <: Complex ? 4 + 3im : 4