# Iterating over elements of upper triangular matrix, but cartesian indices are needed

Basically what I (think) I need is this function:

Of course I can just implement it, two lines of code. But that seems something that probably is needed quite frequently in matrix computations, so I am surprised I didn’t find it as a Base function.

First I thought that most people would prefer constructing an iterator like (for a 3x3 matrix, for example):

``````julia> for c in Iterators.filter(c -> c[1] < c[2], CartesianIndices((1:3,1:3)))
@show c
end
c = CartesianIndex(1, 2)
c = CartesianIndex(1, 3)
c = CartesianIndex(2, 3)
``````

That would be fine for my current purposes, except that it does not play well with `@threads`:

``````julia> Threads.@threads for c in Iterators.filter(c -> c[1] < c[2], CartesianIndices((1:3,1:3)))
@show c
end
Stacktrace:
[1] wait
[3] top-level scope

nested task error: MethodError: no method matching length(::Base.Iterators.Filter{var"#19#20", CartesianIndices{2, Tuple{UnitRange{Int64}, UnitRange{Int64}}}})

``````

What I am missing here? Should I just use a custom implementation of the indexing, is there a Base function already for that, or some alternative way to iterate that is accepted by `@threads`?

I don’t think there’s a built-in iterator for this, although it’s been discussed (including for sparse arrays) e.g. here but probably elsewhere too. Maybe someone knows a better link.

While it won’t directly help with threading, something like `(CartesianIndex(x,y) for x in 1:3 for y in x+1:3)` (i.e Iterators.flatten) may be more efficient than filtering. You can use `Iterators.partition` and `Threads.@spawn` to multi-thread this.

3 Likes

Thanks. The custom function resulted to be the most practical alternative for now. For the records, it is this one:

``````@inline function iuppert(k::Integer,n::Integer)
i = n - 1 - floor(Int,sqrt(-8*k + 4*n*(n-1) + 1)/2 - 0.5)
j = k + i + ( (n-i+1)*(n-i) - n*(n-1) )÷2
return i, j
end
``````
``````julia> iuppert(1,3)
(1,2)

julia> iuppert(3,3)
(2,2)

julia> iuppert(3,4)
(1,4)

``````

I also needed an iterator over indices of an upper triangular matrix, so after finding this thread, I wrote this small package: GitHub - perrutquist/TriangularIndices.jl

It uses a slightly different formula (compared to the above) which generates the indices for an upper triangle in column-first order.