From 3fa31c53c69a982e06fbacf1ca6c412f13dd8829 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Thu, 24 Apr 2025 00:21:49 +0530
Subject: [PATCH 01/13] Bounds-checking in triangular indexing branches

---
 src/triangular.jl  | 65 ++++++++++++++++++++++++++++++++++++----------
 test/triangular.jl | 47 +++++++++++++++++++++++++++++++++
 2 files changed, 98 insertions(+), 14 deletions(-)

diff --git a/src/triangular.jl b/src/triangular.jl
index 26ad4204..6e9120b1 100644
--- a/src/triangular.jl
+++ b/src/triangular.jl
@@ -238,10 +238,22 @@ Base.isassigned(A::UpperOrLowerTriangular, i::Int, j::Int) =
 Base.isstored(A::UpperOrLowerTriangular, i::Int, j::Int) =
     _shouldforwardindex(A, i, j) ? Base.isstored(A.data, i, j) : false
 
-@propagate_inbounds getindex(A::Union{UnitLowerTriangular{T}, UnitUpperTriangular{T}}, i::Int, j::Int) where {T} =
-    _shouldforwardindex(A, i, j) ? A.data[i,j] : ifelse(i == j, oneunit(T), zero(T))
-@propagate_inbounds getindex(A::Union{LowerTriangular, UpperTriangular}, i::Int, j::Int) =
-    _shouldforwardindex(A, i, j) ? A.data[i,j] : diagzero(A,i,j)
+@propagate_inbounds function getindex(A::Union{UnitLowerTriangular{T}, UnitUpperTriangular{T}}, i::Int, j::Int) where {T}
+    if _shouldforwardindex(A, i, j)
+        A.data[i,j]
+    else
+        @boundscheck checkbounds(A, i, j)
+        ifelse(i == j, oneunit(T), zero(T))
+    end
+end
+@propagate_inbounds function getindex(A::Union{LowerTriangular, UpperTriangular}, i::Int, j::Int)
+    if _shouldforwardindex(A, i, j)
+        A.data[i,j]
+    else
+        @boundscheck checkbounds(A, i, j)
+        @inbounds diagzero(A,i,j)
+    end
+end
 
 _shouldforwardindex(U::UpperTriangular, b::BandIndex) = b.band >= 0
 _shouldforwardindex(U::LowerTriangular, b::BandIndex) = b.band <= 0
@@ -250,10 +262,20 @@ _shouldforwardindex(U::UnitLowerTriangular, b::BandIndex) = b.band < 0
 
 # these specialized getindex methods enable constant-propagation of the band
 Base.@constprop :aggressive @propagate_inbounds function getindex(A::Union{UnitLowerTriangular{T}, UnitUpperTriangular{T}}, b::BandIndex) where {T}
-    _shouldforwardindex(A, b) ? A.data[b] : ifelse(b.band == 0, oneunit(T), zero(T))
+    if _shouldforwardindex(A, b)
+        A.data[b]
+    else
+        @boundscheck checkbounds(A, b)
+        ifelse(b.band == 0, oneunit(T), zero(T))
+    end
 end
 Base.@constprop :aggressive @propagate_inbounds function getindex(A::Union{LowerTriangular, UpperTriangular}, b::BandIndex)
-    _shouldforwardindex(A, b) ? A.data[b] : diagzero(A.data, b)
+    if _shouldforwardindex(A, b)
+        A.data[b]
+    else
+        @boundscheck checkbounds(A, b)
+        @inbounds diagzero(A, b)
+    end
 end
 
 _zero_triangular_half_str(::Type{<:UpperOrUnitUpperTriangular}) = "lower"
@@ -265,14 +287,20 @@ _zero_triangular_half_str(::Type{<:LowerOrUnitLowerTriangular}) = "upper"
     throw(ArgumentError(
         lazy"cannot set index in the $Ts triangular part ($i, $j) of an $Tn matrix to a nonzero value ($x)"))
 end
-@noinline function throw_nononeerror(T, @nospecialize(x), i, j)
+@noinline function throw_nonuniterror(T, @nospecialize(x), i, j)
+    check_compatible_type(T, x)
     Tn = nameof(T)
     throw(ArgumentError(
         lazy"cannot set index on the diagonal ($i, $j) of an $Tn matrix to a non-unit value ($x)"))
 end
+function check_compatible_type(T, @nospecialize(x))
+    ET = eltype(T)
+    convert(ET, x) # check that the types are compatible with setindex!
+end
 
 @propagate_inbounds function setindex!(A::UpperTriangular, x, i::Integer, j::Integer)
     if i > j
+        @boundscheck checkbounds(A, i, j)
         iszero(x) || throw_nonzeroerror(typeof(A), x, i, j)
     else
         A.data[i,j] = x
@@ -282,9 +310,11 @@ end
 
 @propagate_inbounds function setindex!(A::UnitUpperTriangular, x, i::Integer, j::Integer)
     if i > j
+        @boundscheck checkbounds(A, i, j)
         iszero(x) || throw_nonzeroerror(typeof(A), x, i, j)
     elseif i == j
-        x == oneunit(x) || throw_nononeerror(typeof(A), x, i, j)
+        @boundscheck checkbounds(A, i, j)
+        x == oneunit(eltype(A)) || throw_nonuniterror(typeof(A), x, i, j)
     else
         A.data[i,j] = x
     end
@@ -293,6 +323,7 @@ end
 
 @propagate_inbounds function setindex!(A::LowerTriangular, x, i::Integer, j::Integer)
     if i < j
+        @boundscheck checkbounds(A, i, j)
         iszero(x) || throw_nonzeroerror(typeof(A), x, i, j)
     else
         A.data[i,j] = x
@@ -302,9 +333,11 @@ end
 
 @propagate_inbounds function setindex!(A::UnitLowerTriangular, x, i::Integer, j::Integer)
     if i < j
+        @boundscheck checkbounds(A, i, j)
         iszero(x) || throw_nonzeroerror(typeof(A), x, i, j)
     elseif i == j
-        x == oneunit(x) || throw_nononeerror(typeof(A), x, i, j)
+        @boundscheck checkbounds(A, i, j)
+        x == oneunit(eltype(A)) || throw_nonuniterror(typeof(A), x, i, j)
     else
         A.data[i,j] = x
     end
@@ -560,7 +593,7 @@ for (T, UT) in ((:UpperTriangular, :UnitUpperTriangular), (:LowerTriangular, :Un
     @eval @inline function _copy!(A::$UT, B::$T)
         for dind in diagind(A, IndexStyle(A))
             if A[dind] != B[dind]
-                throw_nononeerror(typeof(A), B[dind], Tuple(dind)...)
+                throw_nonuniterror(typeof(A), B[dind], Tuple(dind)...)
             end
         end
         _copy!($T(parent(A)), B)
@@ -741,7 +774,7 @@ function _triscale!(A::UpperOrUnitUpperTriangular, B::UnitUpperTriangular, c::Nu
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, c, A, (j,j))
+        @inbounds _modify!(_add, B[BandIndex(0,j)] * c, A, (j,j))
         for i in firstindex(B.data,1):(j - 1)
             @inbounds _modify!(_add, B.data[i,j] * c, A.data, (i,j))
         end
@@ -752,7 +785,7 @@ function _triscale!(A::UpperOrUnitUpperTriangular, c::Number, B::UnitUpperTriang
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, c, A, (j,j))
+        @inbounds _modify!(_add, c * B[BandIndex(0,j)], A, (j,j))
         for i in firstindex(B.data,1):(j - 1)
             @inbounds _modify!(_add, c * B.data[i,j], A.data, (i,j))
         end
@@ -783,7 +816,7 @@ function _triscale!(A::LowerOrUnitLowerTriangular, B::UnitLowerTriangular, c::Nu
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, c, A, (j,j))
+        @inbounds _modify!(_add, B[BandIndex(0,j)] *c, A, (j,j))
         for i in (j + 1):lastindex(B.data,1)
             @inbounds _modify!(_add, B.data[i,j] * c, A.data, (i,j))
         end
@@ -794,7 +827,7 @@ function _triscale!(A::LowerOrUnitLowerTriangular, c::Number, B::UnitLowerTriang
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, c, A, (j,j))
+        @inbounds _modify!(_add, c * B[BandIndex(0,j)], A, (j,j))
         for i in (j + 1):lastindex(B.data,1)
             @inbounds _modify!(_add, c * B.data[i,j], A.data, (i,j))
         end
@@ -804,6 +837,7 @@ end
 
 function _trirdiv!(A::UpperTriangular, B::UpperTriangular, c::Number)
     checksize1(A, B)
+    isunit = B isa UnitUpperTriangular
     for j in axes(B,2)
         for i in firstindex(B,1):j
             @inbounds A.data[i, j] = B.data[i, j] / c
@@ -813,6 +847,7 @@ function _trirdiv!(A::UpperTriangular, B::UpperTriangular, c::Number)
 end
 function _trirdiv!(A::LowerTriangular, B::LowerTriangular, c::Number)
     checksize1(A, B)
+    isunit = B isa UnitLowerTriangular
     for j in axes(B,2)
         for i in j:lastindex(B,1)
             @inbounds A.data[i, j] = B.data[i, j] / c
@@ -822,6 +857,7 @@ function _trirdiv!(A::LowerTriangular, B::LowerTriangular, c::Number)
 end
 function _trildiv!(A::UpperTriangular, c::Number, B::UpperTriangular)
     checksize1(A, B)
+    isunit = B isa UnitUpperTriangular
     for j in axes(B,2)
         for i in firstindex(B,1):j
             @inbounds A.data[i, j] = c \ B.data[i, j]
@@ -831,6 +867,7 @@ function _trildiv!(A::UpperTriangular, c::Number, B::UpperTriangular)
 end
 function _trildiv!(A::LowerTriangular, c::Number, B::LowerTriangular)
     checksize1(A, B)
+    isunit = B isa UnitLowerTriangular
     for j in axes(B,2)
         for i in j:lastindex(B,1)
             @inbounds A.data[i, j] = c \ B.data[i, j]
diff --git a/test/triangular.jl b/test/triangular.jl
index 21fc2d69..ee2a20a3 100644
--- a/test/triangular.jl
+++ b/test/triangular.jl
@@ -966,4 +966,51 @@ end
     end
 end
 
+@testset "indexing checks" begin
+    @testset "getindex" begin
+        U = UnitUpperTriangular(P)
+        @test_throws BoundsError U[0,0]
+        @test_throws BoundsError U[1,0]
+        @test_throws BoundsError U[BandIndex(0,0)]
+        @test_throws BoundsError U[BandIndex(-1,0)]
+
+        U = UpperTriangular(P)
+        @test_throws BoundsError U[1,0]
+        @test_throws BoundsError U[BandIndex(-1,0)]
+
+        L = UnitLowerTriangular(P)
+        @test_throws BoundsError L[0,0]
+        @test_throws BoundsError L[0,1]
+        @test_throws BoundsError U[BandIndex(0,0)]
+        @test_throws BoundsError U[BandIndex(1,0)]
+
+        L = LowerTriangular(P)
+        @test_throws BoundsError L[0,1]
+        @test_throws BoundsError L[BandIndex(1,0)]
+    end
+    @testset "setindex!" begin
+        P = [1 2; 3 4]
+        A = SizedArrays.SizedArray{(2,2)}(P)
+        M = fill(A, 2, 2)
+        U = UnitUpperTriangular(M)
+        @test_throws "Cannot `convert` an object of type Int64" U[1,1] = 1
+        L = UnitLowerTriangular(M)
+        @test_throws "Cannot `convert` an object of type Int64" L[1,1] = 1
+
+        U = UnitUpperTriangular(P)
+        @test_throws BoundsError U[0,0] = 1
+        @test_throws BoundsError U[1,0] = 0
+
+        U = UpperTriangular(P)
+        @test_throws BoundsError U[1,0] = 0
+
+        L = UnitLowerTriangular(P)
+        @test_throws BoundsError L[0,0] = 1
+        @test_throws BoundsError L[0,1] = 0
+
+        L = LowerTriangular(P)
+        @test_throws BoundsError L[0,1] = 0
+    end
+end
+
 end # module TestTriangular

From bb648e4144c44a57d1f710dfc1b6c1ec2b83f806 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Mon, 28 Apr 2025 15:51:15 +0530
Subject: [PATCH 02/13] Fix test

---
 test/triangular.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/triangular.jl b/test/triangular.jl
index ee2a20a3..c11a4da8 100644
--- a/test/triangular.jl
+++ b/test/triangular.jl
@@ -967,6 +967,7 @@ end
 end
 
 @testset "indexing checks" begin
+    P = [1 2; 3 4]
     @testset "getindex" begin
         U = UnitUpperTriangular(P)
         @test_throws BoundsError U[0,0]
@@ -989,7 +990,6 @@ end
         @test_throws BoundsError L[BandIndex(1,0)]
     end
     @testset "setindex!" begin
-        P = [1 2; 3 4]
         A = SizedArrays.SizedArray{(2,2)}(P)
         M = fill(A, 2, 2)
         U = UnitUpperTriangular(M)

From 1235c89271f54c2032e872c48bd041330b448fae Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Mon, 28 Apr 2025 16:11:26 +0530
Subject: [PATCH 03/13] Do not hardcode 64-bit Int

---
 test/triangular.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/test/triangular.jl b/test/triangular.jl
index c11a4da8..d87747e1 100644
--- a/test/triangular.jl
+++ b/test/triangular.jl
@@ -993,9 +993,9 @@ end
         A = SizedArrays.SizedArray{(2,2)}(P)
         M = fill(A, 2, 2)
         U = UnitUpperTriangular(M)
-        @test_throws "Cannot `convert` an object of type Int64" U[1,1] = 1
+        @test_throws "Cannot `convert` an object of type $Int" U[1,1] = 1
         L = UnitLowerTriangular(M)
-        @test_throws "Cannot `convert` an object of type Int64" L[1,1] = 1
+        @test_throws "Cannot `convert` an object of type $Int" L[1,1] = 1
 
         U = UnitUpperTriangular(P)
         @test_throws BoundsError U[0,0] = 1

From 21d9373448f7f1f7964112d4db38d33f6b4bbb2d Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Mon, 28 Apr 2025 17:39:54 +0530
Subject: [PATCH 04/13] Test for unit triangular `l/rdiv!`

---
 test/triangular.jl | 11 +++++++++++
 1 file changed, 11 insertions(+)

diff --git a/test/triangular.jl b/test/triangular.jl
index d87747e1..a9b4f2ff 100644
--- a/test/triangular.jl
+++ b/test/triangular.jl
@@ -1013,4 +1013,15 @@ end
     end
 end
 
+@testset "unit triangular l/rdiv!" begin
+    A = rand(3,3)
+    @testset for (UT,T) in ((UnitUpperTriangular, UpperTriangular),
+                            (UnitLowerTriangular, LowerTriangular))
+        UnitTri = UT(A)
+        Tri = T(LinearAlgebra.full(UnitTri))
+        @test 2 \ UnitTri ≈ 2 \ Tri
+        @test UnitTri / 2 ≈ Tri / 2
+    end
+end
+
 end # module TestTriangular

From 003ad2670c4266f2b70470585d7607dbb6b1a600 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Sun, 11 May 2025 19:13:54 +0530
Subject: [PATCH 05/13] Revert `(r/l)div!` changes

---
 src/triangular.jl | 4 ----
 1 file changed, 4 deletions(-)

diff --git a/src/triangular.jl b/src/triangular.jl
index 6e9120b1..e0dfc161 100644
--- a/src/triangular.jl
+++ b/src/triangular.jl
@@ -837,7 +837,6 @@ end
 
 function _trirdiv!(A::UpperTriangular, B::UpperTriangular, c::Number)
     checksize1(A, B)
-    isunit = B isa UnitUpperTriangular
     for j in axes(B,2)
         for i in firstindex(B,1):j
             @inbounds A.data[i, j] = B.data[i, j] / c
@@ -847,7 +846,6 @@ function _trirdiv!(A::UpperTriangular, B::UpperTriangular, c::Number)
 end
 function _trirdiv!(A::LowerTriangular, B::LowerTriangular, c::Number)
     checksize1(A, B)
-    isunit = B isa UnitLowerTriangular
     for j in axes(B,2)
         for i in j:lastindex(B,1)
             @inbounds A.data[i, j] = B.data[i, j] / c
@@ -857,7 +855,6 @@ function _trirdiv!(A::LowerTriangular, B::LowerTriangular, c::Number)
 end
 function _trildiv!(A::UpperTriangular, c::Number, B::UpperTriangular)
     checksize1(A, B)
-    isunit = B isa UnitUpperTriangular
     for j in axes(B,2)
         for i in firstindex(B,1):j
             @inbounds A.data[i, j] = c \ B.data[i, j]
@@ -867,7 +864,6 @@ function _trildiv!(A::UpperTriangular, c::Number, B::UpperTriangular)
 end
 function _trildiv!(A::LowerTriangular, c::Number, B::LowerTriangular)
     checksize1(A, B)
-    isunit = B isa UnitLowerTriangular
     for j in axes(B,2)
         for i in j:lastindex(B,1)
             @inbounds A.data[i, j] = c \ B.data[i, j]

From a37b66ead221b0c45e10aa76659d0ddccba9890d Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Sun, 11 May 2025 20:07:19 +0530
Subject: [PATCH 06/13] Revert scaling tests

---
 src/triangular.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/triangular.jl b/src/triangular.jl
index e0dfc161..ead4e5b7 100644
--- a/src/triangular.jl
+++ b/src/triangular.jl
@@ -774,7 +774,7 @@ function _triscale!(A::UpperOrUnitUpperTriangular, B::UnitUpperTriangular, c::Nu
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, B[BandIndex(0,j)] * c, A, (j,j))
+        @inbounds _modify!(_add, c, A, (j,j))
         for i in firstindex(B.data,1):(j - 1)
             @inbounds _modify!(_add, B.data[i,j] * c, A.data, (i,j))
         end
@@ -785,7 +785,7 @@ function _triscale!(A::UpperOrUnitUpperTriangular, c::Number, B::UnitUpperTriang
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, c * B[BandIndex(0,j)], A, (j,j))
+        @inbounds _modify!(_add, c, A, (j,j))
         for i in firstindex(B.data,1):(j - 1)
             @inbounds _modify!(_add, c * B.data[i,j], A.data, (i,j))
         end
@@ -816,7 +816,7 @@ function _triscale!(A::LowerOrUnitLowerTriangular, B::UnitLowerTriangular, c::Nu
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, B[BandIndex(0,j)] *c, A, (j,j))
+        @inbounds _modify!(_add, c, A, (j,j))
         for i in (j + 1):lastindex(B.data,1)
             @inbounds _modify!(_add, B.data[i,j] * c, A.data, (i,j))
         end
@@ -827,7 +827,7 @@ function _triscale!(A::LowerOrUnitLowerTriangular, c::Number, B::UnitLowerTriang
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, c * B[BandIndex(0,j)], A, (j,j))
+        @inbounds _modify!(_add, c, A, (j,j))
         for i in (j + 1):lastindex(B.data,1)
             @inbounds _modify!(_add, c * B.data[i,j], A.data, (i,j))
         end

From e76f28d7380ce0986fc5ae9a088ee29f68d45cdb Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Sun, 11 May 2025 22:43:43 +0530
Subject: [PATCH 07/13] Revert "Revert scaling tests"

This reverts commit 0b47f18a38c9c13ce47c16745967f96321c3596e.
---
 src/triangular.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/triangular.jl b/src/triangular.jl
index ead4e5b7..e0dfc161 100644
--- a/src/triangular.jl
+++ b/src/triangular.jl
@@ -774,7 +774,7 @@ function _triscale!(A::UpperOrUnitUpperTriangular, B::UnitUpperTriangular, c::Nu
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, c, A, (j,j))
+        @inbounds _modify!(_add, B[BandIndex(0,j)] * c, A, (j,j))
         for i in firstindex(B.data,1):(j - 1)
             @inbounds _modify!(_add, B.data[i,j] * c, A.data, (i,j))
         end
@@ -785,7 +785,7 @@ function _triscale!(A::UpperOrUnitUpperTriangular, c::Number, B::UnitUpperTriang
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, c, A, (j,j))
+        @inbounds _modify!(_add, c * B[BandIndex(0,j)], A, (j,j))
         for i in firstindex(B.data,1):(j - 1)
             @inbounds _modify!(_add, c * B.data[i,j], A.data, (i,j))
         end
@@ -816,7 +816,7 @@ function _triscale!(A::LowerOrUnitLowerTriangular, B::UnitLowerTriangular, c::Nu
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, c, A, (j,j))
+        @inbounds _modify!(_add, B[BandIndex(0,j)] *c, A, (j,j))
         for i in (j + 1):lastindex(B.data,1)
             @inbounds _modify!(_add, B.data[i,j] * c, A.data, (i,j))
         end
@@ -827,7 +827,7 @@ function _triscale!(A::LowerOrUnitLowerTriangular, c::Number, B::UnitLowerTriang
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, c, A, (j,j))
+        @inbounds _modify!(_add, c * B[BandIndex(0,j)], A, (j,j))
         for i in (j + 1):lastindex(B.data,1)
             @inbounds _modify!(_add, c * B.data[i,j], A.data, (i,j))
         end

From a91349332dcee9027fc7b4563bafbec4607d0fc3 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Sun, 11 May 2025 22:45:30 +0530
Subject: [PATCH 08/13] Test for scaling changes

---
 test/triangular.jl | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/test/triangular.jl b/test/triangular.jl
index a9b4f2ff..de91ddfd 100644
--- a/test/triangular.jl
+++ b/test/triangular.jl
@@ -950,6 +950,10 @@ end
     @test 2\U == 2\M
     @test U*2 == M*2
     @test 2*U == 2*M
+
+    U2 = copy(U)
+    @test rmul!(U, 1) == U2
+    @test lmul!(1, U) == U2
 end
 
 @testset "scaling partly initialized unit triangular" begin

From ff5139ead7c47c0fe57eb821f6a3d030a20ce775 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Sun, 11 May 2025 22:47:28 +0530
Subject: [PATCH 09/13] Whitespace

---
 src/triangular.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/triangular.jl b/src/triangular.jl
index e0dfc161..c483c746 100644
--- a/src/triangular.jl
+++ b/src/triangular.jl
@@ -816,7 +816,7 @@ function _triscale!(A::LowerOrUnitLowerTriangular, B::UnitLowerTriangular, c::Nu
     checksize1(A, B)
     _iszero_alpha(_add) && return _rmul_or_fill!(A, _add.beta)
     for j in axes(B.data,2)
-        @inbounds _modify!(_add, B[BandIndex(0,j)] *c, A, (j,j))
+        @inbounds _modify!(_add, B[BandIndex(0,j)] * c, A, (j,j))
         for i in (j + 1):lastindex(B.data,1)
             @inbounds _modify!(_add, B.data[i,j] * c, A.data, (i,j))
         end

From 763a99d0d6f10828783f66bae334e8d58396995e Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Thu, 19 Jun 2025 12:46:50 +0530
Subject: [PATCH 10/13] Check in setting zero elements

---
 src/triangular.jl  | 56 +++++++++++++++++++++++-----------------------
 test/triangular.jl | 13 +++++++++++
 2 files changed, 41 insertions(+), 28 deletions(-)

diff --git a/src/triangular.jl b/src/triangular.jl
index c483c746..ac86a823 100644
--- a/src/triangular.jl
+++ b/src/triangular.jl
@@ -278,30 +278,22 @@ Base.@constprop :aggressive @propagate_inbounds function getindex(A::Union{Lower
     end
 end
 
-_zero_triangular_half_str(::Type{<:UpperOrUnitUpperTriangular}) = "lower"
-_zero_triangular_half_str(::Type{<:LowerOrUnitLowerTriangular}) = "upper"
-
-@noinline function throw_nonzeroerror(T, @nospecialize(x), i, j)
-    Ts = _zero_triangular_half_str(T)
-    Tn = nameof(T)
+@noinline function throw_nonzeroerror(Tn, @nospecialize(x), i, j)
+    Ts = Tn in (:UpperTriangular, :UnitUpperTriangular) ? "lower" : "upper"
     throw(ArgumentError(
-        lazy"cannot set index in the $Ts triangular part ($i, $j) of an $Tn matrix to a nonzero value ($x)"))
+        lazy"cannot set index in the $Ts triangular part ($i, $j) of a $Tn matrix to a nonzero value ($x)"))
 end
-@noinline function throw_nonuniterror(T, @nospecialize(x), i, j)
-    check_compatible_type(T, x)
-    Tn = nameof(T)
+@noinline function throw_nonuniterror(Tn, @nospecialize(x), i, j)
     throw(ArgumentError(
-        lazy"cannot set index on the diagonal ($i, $j) of an $Tn matrix to a non-unit value ($x)"))
-end
-function check_compatible_type(T, @nospecialize(x))
-    ET = eltype(T)
-    convert(ET, x) # check that the types are compatible with setindex!
+        lazy"cannot set index ($i, $j) on the diagonal of a $Tn matrix to a non-unit value ($x)"))
 end
 
 @propagate_inbounds function setindex!(A::UpperTriangular, x, i::Integer, j::Integer)
     if i > j
         @boundscheck checkbounds(A, i, j)
-        iszero(x) || throw_nonzeroerror(typeof(A), x, i, j)
+        # the value must be convertible to the eltype for setindex! to be meaningful
+        xT = convert(eltype(A), x)
+        iszero(xT) || throw_nonzeroerror(:UpperTriangular, x, i, j)
     else
         A.data[i,j] = x
     end
@@ -309,12 +301,15 @@ end
 end
 
 @propagate_inbounds function setindex!(A::UnitUpperTriangular, x, i::Integer, j::Integer)
-    if i > j
-        @boundscheck checkbounds(A, i, j)
-        iszero(x) || throw_nonzeroerror(typeof(A), x, i, j)
-    elseif i == j
+    if i >= j
         @boundscheck checkbounds(A, i, j)
-        x == oneunit(eltype(A)) || throw_nonuniterror(typeof(A), x, i, j)
+        # the value must be convertible to the eltype for setindex! to be meaningful
+        xT = convert(eltype(A), x)
+        if i > j
+            iszero(xT) || throw_nonzeroerror(:UnitUpperTriangular, x, i, j)
+        else
+            xT == oneunit(eltype(A)) || throw_nonuniterror(:UnitUpperTriangular, x, i, j)
+        end
     else
         A.data[i,j] = x
     end
@@ -324,7 +319,9 @@ end
 @propagate_inbounds function setindex!(A::LowerTriangular, x, i::Integer, j::Integer)
     if i < j
         @boundscheck checkbounds(A, i, j)
-        iszero(x) || throw_nonzeroerror(typeof(A), x, i, j)
+        # the value must be convertible to the eltype for setindex! to be meaningful
+        xT = convert(eltype(A), x)
+        iszero(xT) || throw_nonzeroerror(:LowerTriangular, x, i, j)
     else
         A.data[i,j] = x
     end
@@ -332,12 +329,15 @@ end
 end
 
 @propagate_inbounds function setindex!(A::UnitLowerTriangular, x, i::Integer, j::Integer)
-    if i < j
+    if i <= j
         @boundscheck checkbounds(A, i, j)
-        iszero(x) || throw_nonzeroerror(typeof(A), x, i, j)
-    elseif i == j
-        @boundscheck checkbounds(A, i, j)
-        x == oneunit(eltype(A)) || throw_nonuniterror(typeof(A), x, i, j)
+        # the value must be convertible to the eltype for setindex! to be meaningful
+        xT = convert(eltype(A), x)
+        if i < j
+            iszero(xT) || throw_nonzeroerror(:UnitLowerTriangular, x, i, j)
+        else
+            xT == oneunit(eltype(A)) || throw_nonuniterror(:UnitLowerTriangular, x, i, j)
+        end
     else
         A.data[i,j] = x
     end
@@ -593,7 +593,7 @@ for (T, UT) in ((:UpperTriangular, :UnitUpperTriangular), (:LowerTriangular, :Un
     @eval @inline function _copy!(A::$UT, B::$T)
         for dind in diagind(A, IndexStyle(A))
             if A[dind] != B[dind]
-                throw_nonuniterror(typeof(A), B[dind], Tuple(dind)...)
+                throw_nonuniterror(nameof(typeof(A)), B[dind], Tuple(dind)...)
             end
         end
         _copy!($T(parent(A)), B)
diff --git a/test/triangular.jl b/test/triangular.jl
index de91ddfd..2ef3ec34 100644
--- a/test/triangular.jl
+++ b/test/triangular.jl
@@ -998,8 +998,21 @@ end
         M = fill(A, 2, 2)
         U = UnitUpperTriangular(M)
         @test_throws "Cannot `convert` an object of type $Int" U[1,1] = 1
+        non_unit_msg = "cannot set index $((1,1)) on the diagonal of a UnitUpperTriangular matrix to a non-unit value"
+        @test_throws non_unit_msg U[1,1] = A
         L = UnitLowerTriangular(M)
         @test_throws "Cannot `convert` an object of type $Int" L[1,1] = 1
+        non_unit_msg = "cannot set index $((1,1)) on the diagonal of a UnitLowerTriangular matrix to a non-unit value"
+        @test_throws non_unit_msg L[1,1] = A
+
+        for UT in (UnitUpperTriangular, UpperTriangular)
+            U = UT(M)
+            @test_throws "Cannot `convert` an object of type $Int" U[2,1] = 0
+        end
+        for LT in (UnitLowerTriangular, LowerTriangular)
+            L = LT(M)
+            @test_throws "Cannot `convert` an object of type $Int" L[1,2] = 0
+        end
 
         U = UnitUpperTriangular(P)
         @test_throws BoundsError U[0,0] = 1

From b19d33b27b013c2c8c2166677a0f6a99f40a3d24 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Thu, 19 Jun 2025 13:11:42 +0530
Subject: [PATCH 11/13] Use `_shouldforwardindex` in `setindex!`

---
 src/triangular.jl | 52 +++++++++++++++++++++++++++--------------------
 1 file changed, 30 insertions(+), 22 deletions(-)

diff --git a/src/triangular.jl b/src/triangular.jl
index ac86a823..da52ce27 100644
--- a/src/triangular.jl
+++ b/src/triangular.jl
@@ -289,57 +289,65 @@ end
 end
 
 @propagate_inbounds function setindex!(A::UpperTriangular, x, i::Integer, j::Integer)
-    if i > j
+    if _shouldforwardindex(A, i, j)
+        A.data[i,j] = x
+    else
         @boundscheck checkbounds(A, i, j)
         # the value must be convertible to the eltype for setindex! to be meaningful
-        xT = convert(eltype(A), x)
-        iszero(xT) || throw_nonzeroerror(:UpperTriangular, x, i, j)
-    else
-        A.data[i,j] = x
+        # however, the converted value is unused, and the compiler is free to remove
+        # the conversion if the call is guaranteed to succeed
+        convert(eltype(A), x)
+        iszero(x) || throw_nonzeroerror(nameof(typeof(A)), x, i, j)
     end
     return A
 end
 
 @propagate_inbounds function setindex!(A::UnitUpperTriangular, x, i::Integer, j::Integer)
-    if i >= j
+    if _shouldforwardindex(A, i, j)
+        A.data[i,j] = x
+    else
         @boundscheck checkbounds(A, i, j)
         # the value must be convertible to the eltype for setindex! to be meaningful
-        xT = convert(eltype(A), x)
+        # however, the converted value is unused, and the compiler is free to remove
+        # the conversion if the call is guaranteed to succeed
+        convert(eltype(A), x)
         if i > j
-            iszero(xT) || throw_nonzeroerror(:UnitUpperTriangular, x, i, j)
+            iszero(x) || throw_nonzeroerror(nameof(typeof(A)), x, i, j)
         else
-            xT == oneunit(eltype(A)) || throw_nonuniterror(:UnitUpperTriangular, x, i, j)
+            x == oneunit(eltype(A)) || throw_nonuniterror(nameof(typeof(A)), x, i, j)
         end
-    else
-        A.data[i,j] = x
     end
     return A
 end
 
 @propagate_inbounds function setindex!(A::LowerTriangular, x, i::Integer, j::Integer)
-    if i < j
+    if _shouldforwardindex(A, i, j)
+        A.data[i,j] = x
+    else
         @boundscheck checkbounds(A, i, j)
         # the value must be convertible to the eltype for setindex! to be meaningful
-        xT = convert(eltype(A), x)
-        iszero(xT) || throw_nonzeroerror(:LowerTriangular, x, i, j)
-    else
-        A.data[i,j] = x
+        # however, the converted value is unused, and the compiler is free to remove
+        # the conversion if the call is guaranteed to succeed
+        convert(eltype(A), x)
+        iszero(x) || throw_nonzeroerror(nameof(typeof(A)), x, i, j)
     end
     return A
 end
 
 @propagate_inbounds function setindex!(A::UnitLowerTriangular, x, i::Integer, j::Integer)
-    if i <= j
+    if _shouldforwardindex(A, i, j)
+        A.data[i,j] = x
+    else
         @boundscheck checkbounds(A, i, j)
         # the value must be convertible to the eltype for setindex! to be meaningful
-        xT = convert(eltype(A), x)
+        # however, the converted value is unused, and the compiler is free to remove
+        # the conversion if the call is guaranteed to succeed
+        convert(eltype(A), x)
         if i < j
-            iszero(xT) || throw_nonzeroerror(:UnitLowerTriangular, x, i, j)
+            iszero(x) || throw_nonzeroerror(nameof(typeof(A)), x, i, j)
         else
-            xT == oneunit(eltype(A)) || throw_nonuniterror(:UnitLowerTriangular, x, i, j)
+            x == oneunit(eltype(A)) || throw_nonuniterror(nameof(typeof(A)), x, i, j)
         end
-    else
-        A.data[i,j] = x
     end
     return A
 end

From 70a23a6d942e2f9199fcde71961f25a5e224bb69 Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Thu, 19 Jun 2025 13:26:14 +0530
Subject: [PATCH 12/13] Show index earlier in `throw_nonzeroerror`

---
 src/triangular.jl  | 13 +++++++++----
 test/triangular.jl |  4 ++--
 2 files changed, 11 insertions(+), 6 deletions(-)

diff --git a/src/triangular.jl b/src/triangular.jl
index da52ce27..04fb485e 100644
--- a/src/triangular.jl
+++ b/src/triangular.jl
@@ -278,12 +278,17 @@ Base.@constprop :aggressive @propagate_inbounds function getindex(A::Union{Lower
     end
 end
 
-@noinline function throw_nonzeroerror(Tn, @nospecialize(x), i, j)
-    Ts = Tn in (:UpperTriangular, :UnitUpperTriangular) ? "lower" : "upper"
+@noinline function throw_nonzeroerror(Tn::Symbol, @nospecialize(x), i, j)
+    zero_half = Tn in (:UpperTriangular, :UnitUpperTriangular) ? "lower" : "upper"
+    nstr = Tn === :UpperTriangular ? "n" : ""
     throw(ArgumentError(
-        lazy"cannot set index in the $Ts triangular part ($i, $j) of a $Tn matrix to a nonzero value ($x)"))
+        LazyString(
+            lazy"cannot set index ($i, $j) in the $zero_half triangular part ",
+            lazy"of a$nstr $Tn matrix to a nonzero value ($x)")
+        )
+    )
 end
-@noinline function throw_nonuniterror(Tn, @nospecialize(x), i, j)
+@noinline function throw_nonuniterror(Tn::Symbol, @nospecialize(x), i, j)
     throw(ArgumentError(
         lazy"cannot set index ($i, $j) on the diagonal of a $Tn matrix to a non-unit value ($x)"))
 end
diff --git a/test/triangular.jl b/test/triangular.jl
index 2ef3ec34..0ad4f521 100644
--- a/test/triangular.jl
+++ b/test/triangular.jl
@@ -641,11 +641,11 @@ end
         @testset "error message" begin
             A = UpperTriangular(Ap)
             B = UpperTriangular(Bp)
-            @test_throws "cannot set index in the lower triangular part" copyto!(A, B)
+            @test_throws "cannot set index (3, 1) in the lower triangular part" copyto!(A, B)
 
             A = LowerTriangular(Ap)
             B = LowerTriangular(Bp)
-            @test_throws "cannot set index in the upper triangular part" copyto!(A, B)
+            @test_throws "cannot set index (1, 2) in the upper triangular part" copyto!(A, B)
         end
     end
 

From aad49cf5074dbe09654d4df5f36e9758697b765b Mon Sep 17 00:00:00 2001
From: Jishnu Bhattacharya <jishnub.github@gmail.com>
Date: Thu, 19 Jun 2025 13:30:15 +0530
Subject: [PATCH 13/13] Check for diagonal in unittriangular indexing

---
 src/triangular.jl | 12 ++++++------
 1 file changed, 6 insertions(+), 6 deletions(-)

diff --git a/src/triangular.jl b/src/triangular.jl
index 04fb485e..64bf85e4 100644
--- a/src/triangular.jl
+++ b/src/triangular.jl
@@ -316,10 +316,10 @@ end
         # however, the converted value is unused, and the compiler is free to remove
         # the conversion if the call is guaranteed to succeed
         convert(eltype(A), x)
-        if i > j
-            iszero(x) || throw_nonzeroerror(nameof(typeof(A)), x, i, j)
-        else
+        if i == j # diagonal
             x == oneunit(eltype(A)) || throw_nonuniterror(nameof(typeof(A)), x, i, j)
+        else
+            iszero(x) || throw_nonzeroerror(nameof(typeof(A)), x, i, j)
         end
     end
     return A
@@ -348,10 +348,10 @@ end
         # however, the converted value is unused, and the compiler is free to remove
         # the conversion if the call is guaranteed to succeed
         convert(eltype(A), x)
-        if i < j
-            iszero(x) || throw_nonzeroerror(nameof(typeof(A)), x, i, j)
-        else
+        if i == j  # diagonal
             x == oneunit(eltype(A)) || throw_nonuniterror(nameof(typeof(A)), x, i, j)
+        else
+            iszero(x) || throw_nonzeroerror(nameof(typeof(A)), x, i, j)
         end
     end
     return A