Fix scaling unit triangular matrices (#1346)

After this, the following works:
```julia
julia> U = UnitUpperTriangular(fill(SMatrix{2,2}([1 2; 3 4]), 2, 2))
2×2 UnitUpperTriangular{SMatrix{2, 2, Int64, 4}, Matrix{SMatrix{2, 2, Int64, 4}}}:
 [1 0; 0 1]  [1 2; 3 4]
     ⋅       [1 0; 0 1]

julia> U * 2
2×2 UpperTriangular{SMatrix{2, 2, Int64, 4}, Matrix{SMatrix{2, 2, Int64, 4}}}:
 [2 0; 0 2]  [2 4; 6 8]
     ⋅       [2 0; 0 2]

julia> U = UnitUpperTriangular(@smatrix ones(4,4))
4×4 UnitUpperTriangular{Float64, SMatrix{4, 4, Float64, 16}} with indices SOneTo(4)×SOneTo(4):
 1.0  1.0  1.0  1.0
  ⋅   1.0  1.0  1.0
  ⋅    ⋅   1.0  1.0
  ⋅    ⋅    ⋅   1.0

julia> U * 2
4×4 UpperTriangular{Float64, MMatrix{4, 4, Float64, 16}} with indices SOneTo(4)×SOneTo(4):
 2.0  2.0  2.0  2.0
  ⋅   2.0  2.0  2.0
  ⋅    ⋅   2.0  2.0
  ⋅    ⋅    ⋅   2.0
```

Author: ViralBShah

src/triangular.jl

@@ -1323,6 +1323,15 @@ end
 # Generic routines #
 ####################
 
+function _set_diag!(B::UpperOrLowerTriangular, x)
+    # get a mutable array to modify the diagonal
+    Bm = parent(B) isa StridedArray ? B : copy!(similar(B), B)
+    for i in diagind(Bm.data, IndexStyle(Bm.data))
+        Bm.data[i] = x
+    end
+    Bm
+end
+
 for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
                    (LowerTriangular, UnitLowerTriangular))
     tstrided = t{<:Any, <:StridedMaybeAdjOrTransMat}
@@ -1336,10 +1345,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
 
         function (*)(A::$unitt, x::Number)
             B = $t(A.data)*x
-            for i in axes(A, 1)
-                B.data[i,i] = x
-            end
-            return B
+            _set_diag!(B, oneunit(eltype(A)) * x)
         end
 
         (*)(x::Number, A::$t) = $t(x*A.data)
@@ -1351,10 +1357,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
 
         function (*)(x::Number, A::$unitt)
             B = x*$t(A.data)
-            for i in axes(A, 1)
-                B.data[i,i] = x
-            end
-            return B
+            _set_diag!(B, x * oneunit(eltype(A)))
         end
 
         (/)(A::$t, x::Number) = $t(A.data/x)
@@ -1366,11 +1369,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
 
         function (/)(A::$unitt, x::Number)
             B = $t(A.data)/x
-            invx = inv(x)
-            for i in axes(A, 1)
-                B.data[i,i] = invx
-            end
-            return B
+            _set_diag!(B, oneunit(eltype(A)) / x)
         end
 
         (\)(x::Number, A::$t) = $t(x\A.data)
@@ -1382,11 +1381,7 @@ for (t, unitt) in ((UpperTriangular, UnitUpperTriangular),
 
         function (\)(x::Number, A::$unitt)
             B = x\$t(A.data)
-            invx = inv(x)
-            for i in axes(A, 1)
-                B.data[i,i] = invx
-            end
-            return B
+            _set_diag!(B, x \ oneunit(eltype(A)))
         end
     end
 end

test/triangular.jl

@@ -934,4 +934,14 @@ end
     end
 end
 
+@testset "block unit triangular scaling" begin
+    m = SizedArrays.SizedArray{(2,2)}([1 2; 3 4])
+    U = UnitUpperTriangular(fill(m, 4, 4))
+    M = Matrix{eltype(U)}(U)
+    @test U/2 == M/2
+    @test 2\U == 2\M
+    @test U*2 == M*2
+    @test 2*U == 2*M
+end
+
 end # module TestTriangular