Interpret an array expression

Accelerate is an embedded language that distinguishes between vanilla arrays (e.g. in Haskell memory on the CPU) and embedded arrays (e.g. in device memory on a GPU), as well as the computations on both of these. Since Accelerate is an embedded language, programs written in Accelerate are not compiled by the Haskell compiler (GHC). Rather, each Accelerate backend is a runtime compiler which generates and executes parallel SIMD code of the target language at application runtime.

The type constructor Acc represents embedded collective array operations. A term of type Acc a is an Accelerate program which, once executed, will produce a value of type a (an Array or a tuple of Arrays). Collective operations of type Acc a comprise many scalar expressions, wrapped in type constructor Exp, which will be executed in parallel. Although collective operations comprise many scalar operations executed in parallel, scalar operations cannot initiate new collective operations: this stratification between scalar operations in Exp and array operations in Acc helps statically exclude nested data parallelism, which is difficult to execute efficiently on constrained hardware such as GPUs.

A simple example

As a simple example, to compute a vector dot product we can write:

dotp :: Num a => Vector a -> Vector a -> Acc (Scalar a)
dotp xs ys =
  let
      xs' = use xs
      ys' = use ys
  in
  fold (+) 0 ( zipWith (*) xs' ys' )

The function dotp consumes two one-dimensional arrays (Vectors) of values, and produces a single (Scalar) result as output. As the return type is wrapped in the type Acc, we see that it is an embedded Accelerate computation - it will be evaluated in the object language of dynamically generated parallel code, rather than the meta language of vanilla Haskell.

As the arguments to dotp are plain Haskell arrays, to make these available to Accelerate computations they must be embedded with the use function.

An Accelerate backend is used to evaluate the embedded computation and return the result back to vanilla Haskell. Calling the run function of a backend will generate code for the target architecture, compile, and execute it. For example, the following backends are available:

• accelerate-llvm-native: for execution on multicore CPUs
• accelerate-llvm-ptx: for execution on NVIDIA CUDA-capable GPUs

See also Exp, which encapsulates embedded scalar computations.

Avoiding nested parallelism

As mentioned above, embedded scalar computations of type Exp can not initiate further collective operations.

Suppose we wanted to extend our above dotp function to matrix-vector multiplication. First, let's rewrite our dotp function to take Acc arrays as input (which is typically what we want):

dotp :: Num a => Acc (Vector a) -> Acc (Vector a) -> Acc (Scalar a)
dotp xs ys = fold (+) 0 ( zipWith (*) xs ys )

We might then be inclined to lift our dot-product program to the following (incorrect) matrix-vector product, by applying dotp to each row of the input matrix:

mvm_ndp :: Num a => Acc (Matrix a) -> Acc (Vector a) -> Acc (Vector a)
mvm_ndp mat vec =
  let Z :. rows :. cols  = unlift (shape mat)  :: Z :. Exp Int :. Exp Int
  in  generate (index1 rows)
               (\row -> the $ dotp vec (slice mat (lift (row :. All))))

Here, we use generate to create a one-dimensional vector by applying at each index a function to slice out the corresponding row of the matrix to pass to the dotp function. However, since both generate and slice are data-parallel operations, and moreover that slice depends on the argument row given to it by the generate function, this definition requires nested data-parallelism, and is thus not permitted. The clue that this definition is invalid is that in order to create a program which will be accepted by the type checker, we must use the function the to retrieve the result of the dotp operation, effectively concealing that dotp is a collective array computation in order to match the type expected by generate, which is that of scalar expressions. Additionally, since we have fooled the type-checker, this problem will only be discovered at program runtime.

In order to avoid this problem, we can make use of the fact that operations in Accelerate are rank polymorphic. The fold operation reduces along the innermost dimension of an array of arbitrary rank, reducing the rank (dimensionality) of the array by one. Thus, we can replicate the input vector to as many rows there are in the input matrix, and perform the dot-product of the vector with every row simultaneously:

mvm :: A.Num a => Acc (Matrix a) -> Acc (Vector a) -> Acc (Vector a)
mvm mat vec =
  let Z :. rows :. cols = unlift (shape mat) :: Z :. Exp Int :. Exp Int
      vec'              = A.replicate (lift (Z :. rows :. All)) vec
  in
  A.fold (+) 0 ( A.zipWith (*) mat vec' )

Note that the intermediate, replicated array vec' is never actually created in memory; it will be fused directly into the operation which consumes it. We discuss fusion next.

Fusion

Array computations of type Acc will be subject to array fusion; Accelerate will combine individual Acc computations into a single computation, which reduces the number of traversals over the input data and thus improves performance. As such, it is often useful to have some intuition on when fusion should occur.

The main idea is to first partition array operations into two categories:

1. Element-wise operations, such as map, generate, and backpermute. Each element of these operations can be computed independently of all others.
2. Collective operations such as fold, scanl, and stencil. To compute each output element of these operations requires reading multiple elements from the input array(s).

Element-wise operations fuse together whenever the consumer operation uses a single element of the input array. Element-wise operations can both fuse their inputs into themselves, as well be fused into later operations. Both these examples should fuse into a single loop:

If the consumer operation uses more than one element of the input array (typically, via generate indexing an array multiple times), then the input array will be completely evaluated first; no fusion occurs in this case, because fusing the first operation into the second implies duplicating work.

On the other hand, collective operations can fuse their input arrays into themselves, but on output always evaluate to an array; collective operations will not be fused into a later step. For example:

Here the element-wise sequence (use + generate + zipWith) will fuse into a single operation, which then fuses into the collective fold operation. At this point in the program the fold must now be evaluated. In the final step the map reads in the array produced by fold. As there is no fusion between the fold and map steps, this program consists of two "loops"; one for the use + generate + zipWith + fold step, and one for the final map step.

You can see how many operations will be executed in the fused program by Show-ing the Acc program, or by using the debugging option -ddump-dot to save the program as a graphviz DOT file.

As a special note, the operations unzip and reshape, when applied to a real array, are executed in constant time, so in this situation these operations will not be fused.

Tips

• Since Acc represents embedded computations that will only be executed when evaluated by a backend, we can programatically generate these computations using the meta language Haskell; for example, unrolling loops or embedding input values into the generated code.
• It is usually best to keep all intermediate computations in Acc, and only run the computation at the very end to produce the final result. This enables optimisations between intermediate results (e.g. array fusion) and, if the target architecture has a separate memory space, as is the case of GPUs, to prevent excessive data transfers.

Arrays consists of nested tuples of individual Arrays, currently up to 15-elements wide. Accelerate computations can thereby return multiple results.

Run a complete embedded array program using the reference interpreter.

This is runN specialised to an array program of one argument.

Prepare and execute an embedded array program.