How to set all elements in a lower triangular matrix?

General Usage
,
1

I have a lower triangular matrix which gets reused. In one of the steps, I have all of the elements (below the diagonal) set to the same value. How do I make Julia understand that I’m not trying to set values above the diagonal?

MWE:

julia> using LinearAlgebra

julia> JB =  LowerTriangular(zeros(Float64, 3, 3))
3×3 LowerTriangular{Float64, Matrix{Float64}}:
 0.0   ⋅    ⋅
 0.0  0.0   ⋅
 0.0  0.0  0.0

julia> JB .= -5
ERROR: ArgumentError: cannot set indices in the upper triangular part of an LowerTriangular matrix to a nonzero value (-5)
[...]

What I’m looking for is something like

julia> JB .= -5
3×3 LowerTriangular{Float64, Matrix{Float64}}:
 -5.0    ⋅     ⋅
 -5.0  -5.0    ⋅
 -5.0  -5.0  -5.0

where the matrix structure is still preserved, and only the non-zero elements are set.

2 
function setall!(A::LowerTriangular, value)
  for col=1:size(A,2), row=col:size(A,1) A[row,col] = value end
  A
end

setall!(JB, -5.0)

or, if you must use .=

ltindices = [c for c in CartesianIndices(JB) if c[1]>=c[2]]
JB[ltindices] .= -5.0
1 Like
3

Another way can be:

function setall2!(A::LowerTriangular, value)
    A.data .= value
    A
end

This works faster than setall! (see previous post), but at the cost of going into the internal details of LowerTriangular. Note that LowerTriangular doesn’t store the matrix in a compressed form, but instead stores (and ignores) the other half of the matrix.
Doing the index calculations is slower than setting the invisible part of the matrix.

2 Likes
4

Instead of using internals you can just call parent:

julia> T = LowerTriangular(rand(3, 3))
3×3 LowerTriangular{Float64, Matrix{Float64}}:
 0.171761   ⋅         ⋅ 
 0.381324  0.648652   ⋅ 
 0.672454  0.111463  0.15696

julia> parent(T) .= 1
3×3 Matrix{Float64}:
 1.0  1.0  1.0
 1.0  1.0  1.0
 1.0  1.0  1.0

julia> T
3×3 LowerTriangular{Float64, Matrix{Float64}}:
 1.0   ⋅    ⋅ 
 1.0  1.0   ⋅ 
 1.0  1.0  1.0
5 Likes
5

Note that fill!(parent(T), 1) is likely to be even faster

6

I think they end up being the same. The @code_native output looks identical.

7

Yes, this broadcast operation is reduced to a fill!. The latter might be faster to compile, though, but the performance will be identical.