Regular, shape-polymorphic, parallel arrays in Haskell

Edit
Delete
Tags
Autopost

Regular, shape-polymorphic, parallel arrays in Haskell

The paper Regular, shape-polymorphic, parallel arrays in Haskell describes a new Haskell library, named Repa, for multi-dimensional, regular arrays based on our work on Data Parallel Haskell and a delayed array representation avoiding unnecessary data structures. As the following table shows, the sequential performance already gets close to the performance of C programs. More importantly, the Haskell code transparently makes use of multicore hardware — something that is much more difficult in C.

The Haskell code is purely functional, using collective array operations. Here is matrix-matrix multiplication:

mmMult :: (Num e, Elt e)

=> Array DIM2 e -> Array DIM2 e -> Array DIM2 e

mmMult arr brr = sum (zipWith (*) arrRepl brrRepl)

where

trr = force (transpose2D brr)

arrRepl = replicate (Z :.All :.colsB :.All) arr

brrRepl = replicate (Z :.All :.All :.rowsA) trr

(Z:. colsA:. rowsA) = extent arr

(Z:. colsB:. rowsB) = extent brr

In addition to good performance, the library facilitates reuse through the use of shape polymorphism.

Posted

20 responses
Like
Comment

about 2 years ago paurullan responded:

Very interesting, thanks for the link!

about 2 years ago Yakov Zaytsev responded:

What C implementation of routines was taken? What about to compare with LAPACK/BLAS performance to get real world numbers.

about 2 years ago Felipe Lessa responded:

On Section 4.3, 2nd paragraph, it is said "any rank n > 1" but an example of a rank-1 array is given shortly after.

about 2 years ago Leon Zhang responded:

Interesting. I am curious how many cores of your Repa one ?

about 2 years ago Manuel Chakravarty responded:

Yakov, we wrote our own matrix multiplication (taking care to iterate in a cache-friendly manner). The Laplace code, we got off the net.

A comparison to LAPACK/BLAS, while interesting, would be besides the point of the paper. Our aim was to compare implementations that require similar effort: that of a competent programmer writing a basic implementation of the algorithm. In contrast, LAPACK/BLAS has been highly tuned by numerics experts.

about 2 years ago Manuel Chakravarty responded:

Leon, we used an 8-core machine. See the "Benchmarks" section in the paper for machine details. The limited scalability of Laplace appears to be due to saturating the memory bus.

about 2 years ago Manuel Chakravarty responded:

Felipe, that is a typo. It should read "any rank n >= 1".

about 2 years ago Jonathan responded:

This sounds great! Can't wait to see further development of this, and integration of Repa into DPH (or whatever the correct terminology is).

about 2 years ago Sigmoid responded:

That;s pretty hot performance in the singlethreaded sense too. But a comparison with BLAS/LAPACK would also imply a repeatable position that allows us to compare performance.

about 2 years ago Manuel Chakravarty responded:

Sigmoid, I agree that repeatability is important. We are currently bringing the library into shape to be released as an easy to build package on Hackage. This release will include the C implementations of the benchmarks, too.

about 2 years ago Tyson responded:

Nice work.

I'm only part way through your paper so far, but am wondering if there would be something you could do to improve locality. One common example we use in an introduction to HPC is doing matrix matrix multiplication when your memory layout (row major versus column major) matches your loop layout and when it doesn't. It makes a huge performance difference.

One possible way to avoid the row major versus column major problem is to store a large matrix as a series of smaller sub-matrixes, where the sub-matrixes are sized to easily fit in the L1 cache. As many matrix operations decompose nicely into sub-matrix operations, this would also seem to be a good fit for functional style programming.

http://en.wikipedia.org/wiki/Strassen_algorithm

Related, perhaps something like Z-order would be helpful

http://en.wikipedia.org/wiki/Z-order_(curve)

about 2 years ago Tyson responded:

I finished your paper now, and it doesn't seem like my earlier comment is likely that relevant. Overall a very nice piece of work, and I look forward to it being a regular part of Haskell.

Three minor comments

- section 3.2 discusses DArray's without them ever being introduced anywhere,

- something like Redim or Reshape might be a better name for your Slice class as it is used for both slice and replicate, and

- why not have a foldl, foldr, and fold, where the former two are for when order of operations matter, and the latter when they don't.

I was also thinking something like the Slice class trick would give a nice way of specifying which dimensions to apply shape polymorphic functions to. That is, instead of having do an index permutation to being the dimensions I wish to operate to the end, it would be nice to be able to just do

apply (Any :. 0 :. All) (fold (+)) arr

and have the fold apply along the second last dimension.

The integer could be the seed for the fold, or you could introduce a new empty data type to indicate the dimension). Even more generally, you could make it so a subset of dimensions could be specified for applying multi-dimensional polymorphic functions.

about 2 years ago Manuel Chakravarty responded:

Tyson, thank you very much for the feedback. Using a Slice-like approach to determining the dimensions a shape polymorphic function operates on is an interesting idea.

It would be possible to implement that as an additional convenience on top of the presented API. So, I don't think we will discuss it in the paper (given the strict page limit), but it seems like a nice idea for the library.

about 2 years ago Khudyakov responded:

From practical standpoint of view. Would repa interoperate with ”vector“ package? Multitude of incompatible array data types is already annoying.

about 2 years ago Manuel Chakravarty responded:

Khudyakov, the short answer is that it doesn't interoperate with "vector" at the moment, but it will do that in the future.

The longer answer is that both repa and vector are spin offs of the Data Parallel Haskell project and its "dph" package (that is included in GHC). Package vector started life as an experimental new version of the lowest level of the sequential component of package dph and it will sometime in the next months replace the old sequential array component of package dph. Once this happens, repa (which is implemented on top of package dph) will use vector and the arrays will interoperate.

almost 2 years ago David Duke responded:

Elegant design, and very encouraging results. Hoping to try this soon with some of our visualization examples. Three questions related to the draft paper and library:

- in the mmMult example, the extent of arr is Z :. colsA :. rowsA; this doesn't seem to match Sec 4.1 where a 4 x 5 matrix (rows x cols, I assumed) has extent Z :. 4 :. 5 (or indeed, how !: behaves).

- "traverse" doesn't require Elt e' (the result type for the element function), but map, zip, etc do require "Elt t" for the resulting element type. As the output array is delayed, it wasn't clear if this is strictly necessary.

- More generally, suppose I'm working with values sampled over a regular grid; the values themselves will held in a manifest array, but if the grid is regular, the geometric position of grid points could be represented neatly by a delayed array. However, if I want to transform the location data using repa functions, their signatures require an instance of Elt for the point type, even if I never intend to make the "location" array manifest. On first glance, only the "force" function directly requires an element type that is an instance of Elt. While its simple to work around this particular example, is there any scope for decoupling element constraints from array operations (perhaps via an array type class) - perhaps also allowing further representations, such as sub-arrays as mentioned previously by Tyson?

almost 2 years ago Manuel Chakravarty responded:

David,

- cols vs rows, we did change this around and must have failed to change it consistently.

- You are right that some of these contexts are not strictly necessary. However, it does in some cases depend on which of multiple implementations you choose. We felt that it is best to keep the API invariant wrt. to such changes in the implementation of the functions. This is at the cost of some flexibility, but we didn't really see any need for the flexibility and so went for the stable interface.

- The force function requires the class constraints, but there are other cases, too — for example, a prescan if you want to avoid needless recomputation of intermediate values.

We did actually discuss the trade off between having superfluous constraints and an independence of the interface from the implementation. We chose conservative option (stable interface), but this was to a large extent as we are unsure about what the best option is. So, we are happy to hear arguments for a different design.

over 1 year ago jan crow responded:

are the slides of repa presentation available somewhere?
thanks

over 1 year ago Ben Lippmeier responded:

Yes, they are here:
http://www.cse.unsw.edu.au/~benl/talks/repa-icfp2010-slides.pdf

10 months ago Suresh Kamat responded:

7.8s for 1024x1024 matrices?! If you just take the original unmodified matmul code and run it through a polyhedral compiler like pluto (it'll do the tiling for you and generate optimized C code) that can be compiled with any compiler (http://pluto-compiler.sourceforge.net/):>../../polycc matmul.c --tile --l2tile --unroll --ufactor=8 -o matmul.tiled.c
gcc -O3 matmul.tiled.c -o tiled -lm

./tiled
0.591599s

The above is on one core on an Core i7 870 2.93 GHz
This is more than 10 times faster than the number you report. I think your cache blocking is not really giving any benefit. Note that the effort involved in getting the above performance is zero since everything is automatic. You should really be comparing with this as the above represents the state-of-the-art from the imperative side. The 7.8s you get translates to 275 MFLOPs of performance! around 2% of the machine peak for double-precision floating point perf.

Just Testing

About Me

Regular, shape-polymorphic, parallel arrays in Haskell