JuliaLang/julia — pull request #5294 by lindahua (2014-01-04 04:29:08+00:00)
Title: WIP: Fast implementation of reduction along dims

This is a new implementation of reduction along dims, mainly to improve the performance.

Currently, I have implemented a generic code generation function & macro, and a specialized method for sum over contiguous arrays. Key features of the implementation:
- Cache friendly access. The input elements are accessed in the same order as it is stored in the memory, no matter what region the user specify.
- Based on three vector-kernel (instead of just a binary operation). This allows a variety of optimization to be incorporated (e.g. when SIMD is available, we can easily take advantage of it by rewriting the vector kernels using SIMD instructions).
- Functions in this implementation recursively calls over slices of the array of lower-rank. In this way, it no longer needs the trick of caching functions in a dictionary (which the original implementation relies on).
- The implementation is carefully tested.

This implementation is based on the gist: https://gist.github.com/lindahua/8251556
(with some modification to take care of corner cases, such as empty input).

The script in the gist above runs a benchmark to compare the performance this new implementation and the original one. On my macbook pro, it yields 2x - 16x performance gain (depending on the region):

```julia
# the result is obtained by running it 10 times over an array of size (100,100,100,10)
# et0: elapsed seconds using original method
# et1: elapsed seconds using new method
#
region = 1 : et0 = 0.1245, et1 = 0.0605, gain = 2.0587x
region = 2 : et0 = 0.3095, et1 = 0.0885, gain = 3.4970x
region = 3 : et0 = 0.4836, et1 = 0.0708, gain = 6.8286x
region = 4 : et0 = 0.5405, et1 = 0.1278, gain = 4.2284x
region = (1,2) : et0 = 0.1098, et1 = 0.0464, gain = 2.3665x
region = (1,3) : et0 = 0.1959, et1 = 0.0485, gain = 4.0368x
region = (1,4) : et0 = 0.1445, et1 = 0.0502, gain = 2.8805x
region = (2,3) : et0 = 0.5159, et1 = 0.0785, gain = 6.5687x
region = (2,4) : et0 = 0.4115, et1 = 0.0784, gain = 5.2470x
region = (3,4) : et0 = 1.0379, et1 = 0.0645, gain = 16.0804x
region = (1,2,3) : et0 = 0.1063, et1 = 0.0461, gain = 2.3080x
region = (1,2,4) : et0 = 0.1077, et1 = 0.0460, gain = 2.3426x
region = (1,3,4) : et0 = 0.1587, et1 = 0.0476, gain = 3.3352x
region = (2,3,4) : et0 = 1.1473, et1 = 0.0767, gain = 14.9629x
region = (1,2,3,4) : et0 = 0.1110, et1 = 0.0459, gain = 2.4203x
```

This is still work in progress. I will continue adding other methods, such as maximum, minimum, all, any, prod, etc.

I may also write a short document to explain the rationale behind the implementation.