laplace equation benchmark performance

[ViralBShah]

Track performance of the laplace equation benchmark. Discussion here:

https://groups.google.com/d/topic/julia-dev/VDlr3_6Szos/discussion

[ViralBShah]

Benchmark has been added to the perf2 suite.

https://github.com/JuliaLang/julia/tree/master/test/perf2

[JeffBezanson]

I don't think the vectorized version will ever be as fast as the other one.
These performance issues are getting out of hand. They need to be analyzed and narrowed down.

[ViralBShah]

Yes, I was just thinking the same. The vectorized version does far too much work, but there is certainly room for improvement.

-viral

On 17-Aug-2012, at 3:57 PM, Jeff Bezanson wrote:

I don't think the vectorized version will ever be as fast as the other one.
These performance issues are getting out of hand. They need to be analyzed and narrowed down.

—
Reply to this email directly or view it on GitHub.

[ViralBShah]

Original multi-language performance comparison:

    Octave (0):  14.9932 seconds
     GNU R (0):  56.4161 seconds
     Julia (0):  39.1679 seconds
              gcc     intel
  1      :  13.4774     *     NumPy 
  2      :      *     1.6633  Cilk Plus 
  3,   4 :   2.4963   1.9179  Vectorized Fortran (pure) 

[ViralBShah]

Devectorized julia for the same number of iterations (2^15) now runs in 3.8 seconds for me.

Does numpy implement views and ends up evaluating lazily in this case, and hence end up being much faster? Since octave ends up taking the same amount of time, I presume that we just have a lot of room for improvement in julia.

[toivoh]

Yes, numpy indexing with scalars and slices (ranges) always produces views.

[ViralBShah]

Does octave do the same? I didn't think so - since it is matching numpy's performance. I'll write the benchmark in a few other languages for comparison.

[toivoh]

I would imagine that it does copy on write, since it strives to emulate matlab. But since I'm not an octave user, I really don't know.

[StefanKarpinski]

I really think it's arguable that our slices should be views.

[toivoh]

+1
Except for copy on write (which I'm definitely not suggesting),
I'm not aware of any other sane way to get performance with slices.

[JeffBezanson]

I could get on board with that. We either have to make all ops on SubArray fast, or build views in to the Array type. This is actually a tough decision. Of course we want it to be possible to write something like SubArray and have it be fast.
An issue that immediately comes up is that linear indexing is no longer fast with views, though it is insanely fast for plain dense arrays. This makes me think we need a smarter index abstraction so there is a fast and generic way to iterate over more kinds of arrays. I'm willing to go to extreme lengths to get this. Syntax, changes to the compiler IR, everything is on the table. Rewriting array.jl is no big deal. We can do that in maybe a day.

[timholy]

An issue that immediately comes up is that linear indexing is no longer fast with views

In my currently-stalled rewrite of the Image library, I have a substitute brewing:

type Counter{N}
    max::NTuple{N, Int}
end
Counter(sz) = Counter{length(sz)}(to_tuple(sz))

function start{N}(c::Counter{N})
    state = ones(Int,N)
    state[1] = 0 # because of start/done/next sequence, start out one behind
    return state
end
function done(c::Counter, state)
    # we do the increment as part of "done" to make exit-testing more efficient
    state[1] += 1
    i = 1
    max = c.max
    while state[i] > max[i] && i < length(state)
        state[i] = 1
        i += 1
        state[i] += 1
    end
    state[end] > c.max[end]
end
function next(c::Counter, state)
    return state, state
end

This can be pretty fast, particularly if you handle the inner loop like this:

szv = [sz...]  # vector from tuple
szv[1] = 1
for c = Counter(szv)
    cprod = prod(c)
    for i = 1:sz[1]
        p[index] = i*cprod
        index += 1
    end
end

Then for most applications the performance is essentially indistinguishable from linear indexing. (However, there are traps: this is the usage that prompted me to report performance issue 73d8ca4). Assuming it could be generalized to work on general ranges, and to return a linear index, would that help?

I think in my initial version of this I tried returning both the linear index and the coordinate representation, and also ran into performance problems due to tuple creation, but memory of exactly what was going on is now hazy.

[timholy]

Should have specified: on the inner loop I was demoing how you'd set each element to the product of coordinates, i.e., p[i,j,k] = i*j*k. And you'd need to set index = 1 at the beginning.

[JeffBezanson]

This seems to rely on making the first dimension a special case for performance, or at least assuming it is large enough for the inner loop to dominate. I think we need something more drastic. We may also need to depend on better compiler optimizations like allocating tuples on the stack or in registers.