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
ERROR: TaskFailedException
Stacktrace:
 [1] wait
   @ ./task.jl:322 [inlined]
 [2] threading_run(func::Function)
   @ Base.Threads ./threadingconstructs.jl:34
 [3] top-level scope
   @ ./threadingconstructs.jl:93

    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.