if <math>A=\begin{bmatrix}a&b\\c&d\end{bmatrix}</math> and <math>B=\begin{bmatrix}e&f\\g&h\end{bmatrix}</math>

−

<math>B=\begin{bmatrix}e&f\\g&h\end{bmatrix}</math> then </ref><math>AB=\begin{bmatrix}a&b\\c&d\end{bmatrix}\begin{bmatrix}e&f\\g&h\end{bmatrix}=\begin{bmatrix}ae+bg&af+bh\\ce+dg&cf+dh\end{bmatrix}</math>

Revision as of 01:17, 10 May 2011

Repa is a Haskell library for high performance, regular, multi-dimensional parallel arrays. All numeric data is stored unboxed. Functions written with the Repa combinators are automatically parallel provided you supply "+RTS -N" on the command line when running the program.

This document provides a tutorial on array programming in Haskell using the repa package.

1.2 Index types and shapes

Before we can get started manipulating arrays, we need a grasp of repa's notion of array shape. Much like the classic 'array' library in Haskell, repa-based arrays are parameterized via a type which determines the dimension of the array, and the type of its index. However, while classic arrays take tuples to represent multiple dimensions, Repa arrays use a richer type language for array indices and shapes.

Many operations over arrays are polymorphic in the shape / dimension
component. Others require operating on the shape itself, rather than
the array. A typeclass, Shape, lets us operate uniformally
over arrays with different shape.

1.3 Shapes

To build values of `shape` type, we can use the `Z` and `:.` constructors:

> Z -- the zero-dimension
Z

For arrays of non-zero dimension, we must give a size. A common error is to leave off the type of the size,

> :t Z :. 10
Z :. 10 :: Num head => Z :. head

For arrays of non-zero dimension, we must give a size. A common error is to leave off the type of the size,

> :t Z :. 10
Z :. 10 :: Num head => Z :. head

leading to annoying type errors about unresolved instances, such as:

No instance for (Shape (Z :. head0))

To select the correct instance, you will need to annotate the size literals with their concrete type:

> :t Z :. (10 :: Int)
Z :. (10 :: Int) :: Z :. Int

is the shape of 1D arrays of length 10, indexed via Ints.

Given an array, you can always find its shape by calling extent.

1.4 Generating arrays

New repa arrays ("arrays" from here on) can be generated in many ways:

Repa arrays are instances of the Num. This means that
operations on numerical elements are lifted automagically onto arrays of
such elements. For example, (+) on two double values corresponds to
element-wise addition, (+), of the two arrays of doubles:

The swap function reorders the index space of the array.
To do this, we extract the current shape of the array, and write a function
that maps the index space from the old array to the new array. That index space function
is then passed to backpermute which actually constructs the new
array from the old one.

backpermute generates a new array from an old, when given the new shape, and a
function that translates between the index space of each array (i.e. a shape
transformer).

the type indicate that it works on the lowest two dimensions of the
array.

Other operations on index spaces include taking slices and joining
arrays into larger ones.

1.7.1 Example: matrix-matrix multiplication

A more advanced example from the Repa paper: matrix-matrix multiplication: the result of
matrix multiplication is a matrix whose elements are found by
multiplying the elements of each row from the first matrix by the
associated elements of the same column from the second matrix and
summing the result.

The idea is to expand both 2D argument arrays into 3D arrays by
replicating them across a new axis. The front face of the cuboid that
results represents the array a, which we replicate as often
as b has columns (colsB), producing
aRepl.

The top face represents t (the transposed b), which we
replicate as often as a has rows (rowsA), producing
bRepl,. The two replicated arrays have the same extent,
which corresponds to the index space of matrix multiplication