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
    https://github.com/JuliaLang/julia/issues/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

Author: jishnub

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,14 @@ 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
+
+Base.dataids(A::UpperOrLowerTriangular) = Base.dataids(A.data)
+
 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,255 @@ end
     @test v * S isa Matrix
 end
 
+@testset "copyto! between matrix types" begin
+    dl, d, du = zeros(Int,4), [1:5;], zeros(Int,4)
+    d_ones = ones(Int,size(du))
+
+    @testset "from Diagonal" begin
+        D = Diagonal(d)
+        @testset "to Bidiagonal" begin
+            BU = Bidiagonal(similar(d, BigInt), similar(du, BigInt), :U)
+            BL = Bidiagonal(similar(d, BigInt), similar(dl, BigInt), :L)
+            for B in (BL, BU)
+                copyto!(B, D)
+                @test B == D
+            end
+
+            @testset "mismatched size" begin
+                for B in (BU, BL)
+                    B .= 0
+                    copyto!(B, Diagonal(Int[1]))
+                    @test B[1,1] == 1
+                    B[1,1] = 0
+                    @test iszero(B)
+                end
+            end
+        end
+        @testset "to Tridiagonal" begin
+            T = Tridiagonal(similar(dl, BigInt), similar(d, BigInt), similar(du, BigInt))
+            copyto!(T, D)
+            @test T == D
+
+            @testset "mismatched size" begin
+                T .= 0
+                copyto!(T, Diagonal([1]))
+                @test T[1,1] == 1
+                T[1,1] = 0
+                @test iszero(T)
+            end
+        end
+        @testset "to SymTridiagonal" begin
+            for du2 in (similar(du, BigInt), similar(d, BigInt))
+                S = SymTridiagonal(similar(d), du2)
+                copyto!(S, D)
+                @test S == D
+            end
+
+            @testset "mismatched size" begin
+                S = SymTridiagonal(zero(d), zero(du))
+                copyto!(S, Diagonal([1]))
+                @test S[1,1] == 1
+                S[1,1] = 0
+                @test iszero(S)
+            end
+        end
+    end
+
+    @testset "from Bidiagonal" begin
+        BU = Bidiagonal(d, du, :U)
+        BUones = Bidiagonal(d, oneunit.(du), :U)
+        BL = Bidiagonal(d, dl, :L)
+        BLones = Bidiagonal(d, oneunit.(dl), :L)
+        @testset "to Diagonal" begin
+            D = Diagonal(zero(d))
+            for B in (BL, BU)
+                @test copyto!(D, B) == B
+                D .= 0
+            end
+            for B in (BLones, BUones)
+                errmsg = "cannot copy a Bidiagonal with a non-zero off-diagonal band to a Diagonal"
+                @test_throws errmsg copyto!(D, B)
+                @test iszero(D)
+            end
+
+            @testset "mismatched size" begin
+                for uplo in (:L, :U)
+                    D .= 0
+                    copyto!(D, Bidiagonal(Int[1], Int[], uplo))
+                    @test D[1,1] == 1
+                    D[1,1] = 0
+                    @test iszero(D)
+                end
+            end
+        end
+        @testset "to Tridiagonal" begin
+            T = Tridiagonal(similar(dl, BigInt), similar(d, BigInt), similar(du, BigInt))
+            for B in (BL, BU, BLones, BUones)
+                copyto!(T, B)
+                @test T == B
+            end
+
+            @testset "mismatched size" begin
+                T = Tridiagonal(oneunit.(dl), zero(d), oneunit.(du))
+                for uplo in (:L, :U)
+                    T .= 0
+                    copyto!(T, Bidiagonal([1], Int[], uplo))
+                    @test T[1,1] == 1
+                    T[1,1] = 0
+                    @test iszero(T)
+                end
+            end
+        end
+        @testset "to SymTridiagonal" begin
+            for du2 in (similar(du, BigInt), similar(d, BigInt))
+                S = SymTridiagonal(similar(d, BigInt), du2)
+                for B in (BL, BU)
+                    copyto!(S, B)
+                    @test S == B
+                end
+                errmsg = "cannot copy a non-symmetric Bidiagonal matrix to a SymTridiagonal"
+                @test_throws errmsg copyto!(S, BUones)
+                @test_throws errmsg copyto!(S, BLones)
+            end
+
+            @testset "mismatched size" begin
+                S = SymTridiagonal(zero(d), zero(du))
+                for uplo in (:L, :U)
+                    copyto!(S, Bidiagonal([1], Int[], uplo))
+                    @test S[1,1] == 1
+                    S[1,1] = 0
+                    @test iszero(S)
+                end
+            end
+        end
+    end
+
+    @testset "from Tridiagonal" begin
+        T = Tridiagonal(dl, d, du)
+        TU = Tridiagonal(dl, d, d_ones)
+        TL = Tridiagonal(d_ones, d, dl)
+        @testset "to Diagonal" begin
+            D = Diagonal(zero(d))
+            @test copyto!(D, T) == Diagonal(d)
+            errmsg = "cannot copy a Tridiagonal with a non-zero off-diagonal band to a Diagonal"
+            D .= 0
+            @test_throws errmsg copyto!(D, TU)
+            @test iszero(D)
+            errmsg = "cannot copy a Tridiagonal with a non-zero off-diagonal band to a Diagonal"
+            @test_throws errmsg copyto!(D, TL)
+            @test iszero(D)
+
+            @testset "mismatched size" begin
+                D .= 0
+                copyto!(D, Tridiagonal(Int[], Int[1], Int[]))
+                @test D[1,1] == 1
+                D[1,1] = 0
+                @test iszero(D)
+            end
+        end
+        @testset "to Bidiagonal" begin
+            BU = Bidiagonal(zero(d), zero(du), :U)
+            BL = Bidiagonal(zero(d), zero(du), :L)
+            @test copyto!(BU, T) == Bidiagonal(d, du, :U)
+            @test copyto!(BL, T) == Bidiagonal(d, du, :L)
+
+            BU .= 0
+            BL .= 0
+            errmsg = "cannot copy a Tridiagonal with a non-zero superdiagonal to a Bidiagonal with uplo=:L"
+            @test_throws errmsg copyto!(BL, TU)
+            @test iszero(BL)
+            @test copyto!(BU, TU) == Bidiagonal(d, d_ones, :U)
+
+            BU .= 0
+            BL .= 0
+            @test copyto!(BL, TL) == Bidiagonal(d, d_ones, :L)
+            errmsg = "cannot copy a Tridiagonal with a non-zero subdiagonal to a Bidiagonal with uplo=:U"
+            @test_throws errmsg copyto!(BU, TL)
+            @test iszero(BU)
+
+            @testset "mismatched size" begin
+                for B in (BU, BL)
+                    B .= 0
+                    copyto!(B, Tridiagonal(Int[], Int[1], Int[]))
+                    @test B[1,1] == 1
+                    B[1,1] = 0
+                    @test iszero(B)
+                end
+            end
+        end
+    end
+
+    @testset "from SymTridiagonal" begin
+        S2 = SymTridiagonal(d, ones(Int,size(d)))
+        for S in (SymTridiagonal(d, du), SymTridiagonal(d, zero(d)))
+            @testset "to Diagonal" begin
+                D = Diagonal(zero(d))
+                @test copyto!(D, S) == Diagonal(d)
+                D .= 0
+                errmsg = "cannot copy a SymTridiagonal with a non-zero off-diagonal band to a Diagonal"
+                @test_throws errmsg copyto!(D, S2)
+                @test iszero(D)
+
+                @testset "mismatched size" begin
+                    D .= 0
+                    copyto!(D, SymTridiagonal(Int[1], Int[]))
+                    @test D[1,1] == 1
+                    D[1,1] = 0
+                    @test iszero(D)
+                end
+            end
+            @testset "to Bidiagonal" begin
+                BU = Bidiagonal(zero(d), zero(du), :U)
+                BL = Bidiagonal(zero(d), zero(du), :L)
+                @test copyto!(BU, S) == Bidiagonal(d, du, :U)
+                @test copyto!(BL, S) == Bidiagonal(d, du, :L)
+
+                BU .= 0
+                BL .= 0
+                errmsg = "cannot copy a SymTridiagonal with a non-zero off-diagonal band to a Bidiagonal"
+                @test_throws errmsg copyto!(BU, S2)
+                @test iszero(BU)
+                @test_throws errmsg copyto!(BL, S2)
+                @test iszero(BL)
+
+                @testset "mismatched size" begin
+                    for B in (BU, BL)
+                        B .= 0
+                        copyto!(B, SymTridiagonal(Int[1], Int[]))
+                        @test B[1,1] == 1
+                        B[1,1] = 0
+                        @test iszero(B)
+                    end
+                end
+            end
+        end
+    end
+end
+
+@testset "BandIndex indexing" begin
+    for D in (Diagonal(1:3), Bidiagonal(1:3, 2:3, :U), Bidiagonal(1:3, 2:3, :L),
+                Tridiagonal(2:3, 1:3, 1:2), SymTridiagonal(1:3, 2:3))
+        M = Matrix(D)
+        for band in -size(D,1)+1:size(D,1)-1
+            for idx in 1:size(D,1)-abs(band)
+                @test D[BandIndex(band, idx)] == M[BandIndex(band, idx)]
+            end
+        end
+        @test_throws BoundsError D[BandIndex(size(D,1),1)]
+    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

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,