Haven't tested timing carefully, but something like the following should be
pretty decent:
julia> immutable LowerTriangularIterator{A<:AbstractMatrix}
data::A
end
julia> Base.start(iter::LowerTriangularIterator) = (1,1)
start (generic function with 49 methods)
julia> Base.done(iter::LowerTriangularIterator, state) = state[2] >
size(iter.data, 2)
done (generic function with 51 methods)
julia> @inline function Base.next(iter::LowerTriangularIterator, state)
item = iter.data[state[1], state[2]]
newstate = ifelse(state[1] == size(iter.data, 1), (state[2]+1,
state[2]+1), (state[1]+1, state[2]))
item, newstate
end
To make it vectorize (i.e., SIMD), you might also need definitions like these:
https://github.com/JuliaLang/julia/blob/aa24ee0d0f334e8d9ecbced64bd975a6cc5f240d/base/multidimensional.jl#L134-L150
but perhaps vectorization isn't applicable to your problem anyway.
--Tim
On Sunday, November 15, 2015 05:06:23 PM Seth wrote:
> Thank you, Tim - I'll definitely check it out. Right now I can't seem to
> find an efficient way of doing it, but there are a few tricks in
> Iterators.jl that might work.
>
> On Sunday, November 15, 2015 at 1:33:29 PM UTC-8, Tim Holy wrote:
> > This may also be a case where writing an iterator that just visits the
> > elements with i >= j might be what you want. You can study the
> > implementation
> > of CartesianIndex in Base for inspiration.
> >
> > --Tim
> >
> > On Saturday, November 14, 2015 02:18:45 PM hustf wrote:
> > > I believe this work-in-progress type may be adapted? I wanted to make
> >
> > this
> >
> > > faster by using tuples or immutable fixed vectors to store the data or
> > > whatnot, but I guess I'm currently stuck. Help wanted.
> > >
> > > Compared to looping unnecessarily over a full matrix, this speeds things
> >
> > up
> >
> > > by around 35%.
> > >
> > > immutable Maq7 <:AbstractMatrix{Int}
> > > data::Array{Int64,2} # upper triangular except diagonal, stored as a
> > > horizontal vector
> > > width::Int64
> > > function Maq7(mat::Matrix{Int64})
> > > if issym(mat)
> > > n=size(mat,1)
> > > if n>1
> > > v=mat[2:n,1]
> > > for col = 2:(n-1)
> > > append!(v,mat[(col+1):n,col])
> > > end
> > > new(v', size(mat,1))
> > > else
> > > error("Ouch, not > (1x1)")
> > > end
> > > else
> > > error("Ouch, not square symmetric")
> > > end
> > > end
> > > end
> > > Base.size(A::Maq7)= A.width, A.width
> > > Base.length(A::Maq7)= div(A.width*(A.width-1), 2)
> > > Base.getindex(A::Maq7, i::Int) = Base.getindex(A.data, i)
> > > function Base.getindex(A::Maq7, r::Int , c::Int )
> > > if c > r
> > > return getindex(A.data, div( (r-1)*(2*A.width -r) , 2 ) + c - r )
> > > elseif r > c
> > > return getindex(A.data, div( (c-1)*(2*A.width -c) , 2 ) + r - c )
> > > else
> > > return zero(eltype(A.data))
> > > end
> > > end
> > > Base.linearindexing(::Type{Maq7}) = Base.LinearFast()
> > >
> > >
> > >
> > > Example:
> > > julia> ms
> > >
> > > 4x4 Array{Int64,2}:
> > > 7 1 2 3
> > > 1 8 4 5
> > > 2 4 9 6
> > > 3 5 6 10
> > >
> > > julia> m=Maq7(ms)
> > >
> > > 4x4 Maq7:
> > > 0 1 2 3
> > > 1 0 4 5
> > > 2 4 0 6
> > > 3 5 6 0
> > >
> > > julia> collect(m)
> > >
> > > 6-element Array{Int64,1}:
> > > 1
> > > 2
> > > 3
> > > 4
> > > 5
> > > 6
> > >
> > > julia> length(m)
> > > 6
