`rmul!(A::UnitUpperTriangular, D::Diagonal)` cannot change `A` to be the result

[mcabbott]

Bug noted here exists at least since v1.6.7 (edit: with latest code added by #1228):

julia> using LinearAlgebra

julia> A = UnitUpperTriangular([1 0; 0 1])
2×2 UnitUpperTriangular{Int64, Matrix{Int64}}:
 1  0
 ⋅  1

julia> rmul!(A, Diagonal([2,3]))
2×2 UpperTriangular{Int64, Matrix{Int64}}:
 2  0
 ⋅  3

julia> A
2×2 UnitUpperTriangular{Int64, Matrix{Int64}}:
 1  0
 ⋅  1

julia> @which rmul!(A, Diagonal([2,3]))
rmul!(A::UnitUpperTriangular{<:Any, <:StridedMatrix{T} where T}, D::Diagonal)
     @ LinearAlgebra ~/.julia/juliaup/julia-1.13.0-rc1+0.aarch64.apple.darwin14/Julia-1.13.app/Contents/Resources/julia/share/julia/stdlib/v1.13/LinearAlgebra/src/diagonal.jl:388

Clearer version -- A.data is in fact mutated, but the diagonal of A does not (and cannot) take the mathematically correct value:

julia> A = UnitUpperTriangular([10 10; 10 10])
2×2 UnitUpperTriangular{Int64, Matrix{Int64}}:
 1  10
 ⋅   1

julia> rmul!(A, Diagonal([2,3]))
2×2 UpperTriangular{Int64, Matrix{Int64}}:
 2  30
 ⋅   3

julia> A
2×2 UnitUpperTriangular{Int64, Matrix{Int64}}:
 1  30
 ⋅   1

julia> A.data
2×2 Matrix{Int64}:
  2  30
 10   3

[odow]

I'm not even sure this is a bug. Doing:

julia> A = rmul!(A, Diagonal([2,3]))
2×2 UpperTriangular{Int64, Matrix{Int64}}:
 2  0
 ⋅  3

appropriately converts from UnitUpperTriangular to UpperTriangular. We can't update in-place without changing the type.

There's an argument for throwing an error instead, as-is seems okay.

[mcabbott]

I guess it comes down to whether rmul!'s contract is that it must mutate A to take the value (like mul!), or whether it is simply allowed to mutate it as scratch space (like eigvals!.) Those two have unambiguous documentation:

help?> mul!

  mul!(Y, A, B) -> Y

  Calculates the matrix-matrix or matrix-vector product A B and stores the result in Y,

...

help?> eigvals!

  eigvals!(A; permute::Bool=true, scale::Bool=true, sortby) -> values

  Same as eigvals, but saves space by overwriting the input A, instead of creating a copy. The
  permute, scale, and sortby keywords are the same as for eigen.

  │ Note
  │
  │  The input matrix A will not contain its eigenvalues after eigvals! is called on it - A
  │  is used as a workspace.

Note that other methods of rmul! do give errors if they cannot mutate the array, rather than just returning an UpperTriangular:

julia> A
2×2 UnitUpperTriangular{Int64, Matrix{Int64}}:
 1  30
 ⋅   1

julia> rmul!(A, 1.0)
2×2 UnitUpperTriangular{Int64, Matrix{Int64}}:
 1  30
 ⋅   1

julia> rmul!(A, 33)
ERROR: ArgumentError: cannot set index (1, 1) on the diagonal of a UnitUpperTriangular matrix to a non-unit value (33)
Stacktrace:
 [1] throw_nonuniterror(Tn::Symbol, x::Any, i::Int64, j::Int64)
   @ LinearAlgebra ~/.julia/juliaup/julia-1.13.0-rc1+0.aarch64.apple.darwin14/Julia-1.13.app/Contents/Resources/julia/share/julia/stdlib/v1.13/LinearAlgebra/src/triangular.jl:292
 [2] setindex!
   @ ~/.julia/juliaup/julia-1.13.0-rc1+0.aarch64.apple.darwin14/Julia-1.13.app/Contents/Resources/julia/share/julia/stdlib/v1.13/LinearAlgebra/src/triangular.jl:320 [inlined]

[adienes]

should error IMO

[dkarrasch]

should error IMO

I agree. Most users will call rmul!(A, D) and continue to work with A (and rightfully expect this to work correctly).