\documentclass{llncs}
\usepackage{code}
\usepackage{url}
\usepackage{pstricks,pst-plot}
\usepackage{epsfig}

\newcommand{\simon}[1]{[{\em {\bf SLPJ}: #1}]}
\newcommand{\depth}[1]{{\cal D}[\![ #1 ]\!]}
\newcommand{\maxx}[1]{\raisebox{-1mm}{$\stackrel{\textstyle max}{\scriptstyle #1}$}\,}


%pstricks stuff
\newgray{bggray}{.8}
\newgray{fggray}{.6}
\newgray{fglightgray}{.8}
\newgray{dgray}{.5}

\begin{document}
\makeatactive

\title{Notes on NDP}

\author{Simon Peyton Jones and Satnam Singh}
\maketitle

\section{Nested data parallelism}

\begin{quote}
This chapter was written in collaboration with Manuel Chakravarty, Gabriele Keller, and Roman Leshchinskiy
(University of New South Wales, Sydney).
\end{quote}
The two major ways of exploiting parallelism that we have seen so far
each have their disadvantages:
\begin{itemize}
\item 
The {\tt par}/{\tt seq} style is semantically transparent,
but it is hard to ensure that the granularity is consistently large
enough to be worth spawning new threads.
\item
Explicitly-forked threads, communicating using {\tt MVar}s or STM 
give the programmer precise control over granularity, but at the cost
of a new layer of semantic complexity: there are now many threads,
each mutating shared memory.  Reasoning about all the interleavings of
these threads is hard, especially if there are a lot of them.
\end{itemize}
Furthermore, neither is easy to implement on a distributed-memory machine,
\simon{I'm assuming you've discussed shared vs distributed memory}
because any pointer can point to any value, so spatial locality is poor.  
It is possible to support
this anarchic memory model on a distributed-memory archiecture, as Glasgow
Parallel Haskell has shown \cite{gum-europar00}, but it is very hard to get 
reliable, predicatable, and scalable performance.  In short, we have no good
\emph{performance model}, which is a Bad Thing if your main purpose in writing
a parallel program is to improve performance.

In this chapter we will explore another parallel programming paradigm: \emph{data
parallelism}.  The basic idea of data parallelism is simple:
\begin{quote}
\emph{Do the same thing, in parallel, to every element of a large collection of values.}
\end{quote}
Not every program can be expressed in this way, but data parallelism is very attractive 
for those that can, because:
\begin{itemize}
\item Everything remains purely functional, like {\tt par}/{\tt seq}, 
so there is no new semantic complexity.
\item Granularity is very good: to a first approximation, we get just
one thread (with its attendent overheads) for each physical processor, rather
than one thread for each data item (of which there are zillions).
\item Locality is very good: the data can be physically partitioned across
the processors without random cross-heap pointers.
\end{itemize}
As a result, we get an excellent performance model.  

\subsection{Flat data parallelism}

Data parallelism sounds good doesn't it?  Indeed, data-parallel
programming is widely and successfully used in mainstream languages
such as High-Performance Fortran.  However, there's a catch: the
application has to fit the data-parallel programming paradigm, and
only a fairly narrow class of applications do so.  But this
narrow-ness is largely because mainstream data-parallel technology
only supports so-called \emph{flat} data parallelism.  Flat data
parallelism works like this
\begin{quote}
Apply the same \emph{sequential}
function {\tt f}, in parallel, to every element of a large collection of values {\tt a}.
Not only is {\tt f} sequential, but it has a similar run-time for each element of 
the collection.  
\end{quote}
Here is how we might write such a loop in Data Parallel Haskell:
\begin{verbatim}
  sumSq :: [: Float :] -> Float
  sumSq a = sumP [: x*x | x <- a :]
\end{verbatim}
The data type {\tt [: Float :]} is  pronounced ``parallel vector of {\tt Float}''.  We use
a bracket notation reminiscent of lists, because parallel vectors are similar 
to lists in that consist of an sequence of elements.  Many functions available for lists
are also available for parallel vectors.  For example
\begin{verbatim}
  mapP     :: (a -> b) -> [:a:] -> [:b:]
  zipWithP :: (a -> b -> c) -> [:a:] -> [:b:] -> [:c:]
  sumP     :: Num a => [:a:] -> a

  (+:+)    :: [:a:] -> [:a:] -> [:a:]
  filterP  :: (a -> Bool) -> [:a:] -> [:a:]
  anyP     :: (a -> Bool) -> [:a:] -> Bool
  concatP  :: [:[:a:]:] -> [:a:]
  nullP    :: [:a:] -> Bool
  lengthP  :: [:a:] -> Int
  (!:)     :: [:a:] -> Int -> a   -- Zero-based indexing
\end{verbatim}
These functions, and many more, are exported by {\tt Data.Array.Parallel}.
Just as we have list comprehensions, we also have parallel-array comprehensions, of which one is
used in the above example.  But, just as with list comprehensions, array comprehensions are
syntactic sugar, and we could just as well have written
\begin{verbatim}
  sumSq :: [: Float :] -> Float
  sumSq a = sumP (mapP (\x -> x*x) a)
\end{verbatim}
Notice that there is no @forkIO@, and no @par@. The parallelism comes implicitly from
use of the primitives operating on parallel vectors, such as {\tt mapP}, {\tt sumP}, and so on.

Flat data parallelism is not restricted to consuming a single array.  For example,
here is how we might take the product of two vectors, by multiplying corresponding 
elements and adding up the results:
\begin{verbatim}
  vecMul :: [:Float:] -> [:Float:] -> Float
  vecMul a b = sumP [: x*y | x <- a | y <- b :]
\end{verbatim}
The array comprehension uses a second vertical bar ``@|@'' to indicate that
we interate over @b@ in lockstep with @a@.  (This same facility is available for ordinary
list comprehensions too.)  As before the comprehension is just syntactic sugar, and
we could have equivalently written this:
\begin{verbatim}
  vecMul :: [:Float:] -> [:Float:] -> Float
  vecMul a b = sumP (zipWithP (*) a b)
\end{verbatim}

\subsection{Pros and cons of flat data parallelism}

If you can express your program using flat data parallelism, we can implement it 
really well on a N-processor machine:
\begin{itemize}
\item Divide {\tt a} into N chunks, one for each processor.
\item Compile a sequential loop that applies {\tt f} successively to each element of
a  chunk
\item Run this loop on each procesor
\item Combine the results.
\end{itemize}
Notice that the granularity is good (there is one large-grain thread per processor);
locality is good (the elements of {\tt a} are accessed successively);
load-balancing is good (each processor does $1/N$ of the work).
Furthermore the algorithm works well even if {\tt f} itself does very little work to
each element, a situation that is a killer if we spawn a new thread for each invocation
of {\tt f}.

In exchange for this great implementation, the programming model is horrible: \emph{all the
paralellism must come from a single parallel loop}. 
This restriction makes the programming model is very non-compositional.  If you
have an existing function {\tt g} written using the data-parallel @mapP@, 
you can't call {\tt g} from another data-parallel map (e.g. @mapP g a@),
because the argument to @mapP@ must be a \emph{sequential} function.  

Furthermore, just as the control structure must be flat, so must the
data structure.  We cannot allow {\tt a} to contain rich nested
structure (e.g. the elements of {\tt a} cannot themselves be vectors),
or else similar-run-time promise of {\tt f} could not be guaranteed,
and data locality would be lost.

\subsection{Nested data parallelism}

In the early 90's, Guy Blelloch described \emph{nested} data-parallel programming.
The idea is similar:
\begin{quote}
Apply the same 
function {\tt f}, in parallel, to every element of a large collection of values {\tt a}.
However, {\tt f} may \emph{itself} be a (nested) data-parallel function, and does not need to 
have a similar run-time for each element of {\tt a}.
\end{quote}
For example, here is how we might multiply a matrix by a vector:
\begin{verbatim}
  type Vector = [:Float:]
  type Matrix = [:Vector:]

  matMul :: Matrix -> Vector -> Vector
  matMul m v = [: vecMul r v | r <- m :]
\end{verbatim}
That is, for each row of the matrix, multiply it by the vector @v@ using @vecMul@.
Here we are calling a data-parallel function @vecMul@ from inside a data-parallel operation
(the comprehension in @matMul@).

In very regular examples like this, consisting of visible, nested
loops, modern FORTRAN compilers can collapse a loop nest into one
loop, and partition the loop across the processors.  It is not
entirely trivial to do this, but it is well within the reach of
compiler technology.  But the flattening process only works for the
simplest of cases.  A typical complication is the matrices may be
\emph{sparse}.

A sparse vector (or matrix) is one in which almost all the elements are zero.  We may
represent a sparse vector by a (dense) vector of pairs:
\begin{code}
  type SparseVector = [: (Int, Float) :]
\end{code}
In this representation, only non-zero elements of the vector are represented, by a pair
of their index and value.  A sparse matrix can now be represented by a (dense) vector of
rows, each of which is a sparse vector:
\begin{code}
  type SparseMatrix = [: SparseVector :]
\end{code}
Now we may write @vecMul@ and @matMul@ for sparse arguments thus\footnote{
Incidentally, although these functions are very short, they are important in
some applications.  For example, multiplying a sparse matrix by a dense vector 
(i.e. @sparseMatMul@) is the inner
loop of the NAS Conjugate Gradient benchmark, consuming 95\% of runtime \cite{prins:irregular}.}:
\begin{code}
  sparseVecMul :: SparseVector -> Vector -> Float
  sparseVecMul sv v = sumP [: x * v!:i | (i,x) <- sv :]
 
  sparseMatMul :: SparseMatrix -> Vector -> Vector
  sparseMatMul sm v = [: sparseVecMul r v | r <- sm :]
\end{code}
We use the indexing operator @(!:)@ to index the dense vector @v@.  In
this code, the control structure is the same as before (a nested loop,
with both levels being data-parallel), but now the data structure is
much less regular, and it is \emph{much} less obvious how to flatten the program into 
a single data-parallel loop, in such a way that the work is evenly distributed over
N processors, regardless of the distribution of non-zero data in the matrix.

Blelloch's remarkable contribution was to show that it is possible to
take \emph{any} program written using nested data parallelism (easy to
write but hard to implement efficiently), and transform it
systematically into a program that uses flat data parallelism (hard to
write but easy to implement efficiently).  He did this for a
special-purpose functional language, NESL, designed specifically to
demonstrate nested data parallelism.  

As a practical programming language, however, NESL is very limited: it
is a first-order language, it has only a fixed handful of data types,
it is implemented using an interpreter, and so on.  Fortunately, in a
series of papers, Manuel Chakravarty, Gabriele Keller and Roman
Leshchinskiy have generalised Blelloch's transformation to a modern,
higher order functional programming language with user-defined
algebraic data types -- in other words, Haskell.
Data Parallel Haskell is a research prototype implementation of 
all these ideas, in the Glasgow Haskell Compiler, GHC.

The matrix-multiply examples may have suggested to you that
Data Parallel Haskell is intended primarily for scientific
applications, and that the nesting depth of parallel computations 
is statically fixed.  However the programming paradigm is much
more flexible than that.  In the rest of this chapter we will give a 
series of examples of programming in Data Parallel Haskell, designed to help
you gain familiarity with the programming style.

Most (in due course, all) of these examples can be found at in the Darcs repository
\url{http://darcs.haskell.org/packages/ndp}, in the sub-directory \url{examples/}.
You can also find a dozen or so other examples of data-parallel algorithms written in NESL at 
\url{http://www.cs.cmu.edu/~scandal/nesl/algorithms.html}.

\subsection{Word search}

Here is a tiny verion of a web search engine.  A @Document@ is a
vector of words, each of which is a string.  The task is to find
all the occurences of a word in a large collection of documents,
returning the matched documents and the matching word positions in 
those documents.  So here is the type signature for @search@:
\begin{code}
  type Document = [: String :]
  type DocColl  = [: Document :]
  search :: DocColl -> String -> [: (Document, [:Int:]) :]
\end{code}
We start by solving an easier problem, that of finding
all the occurrences of a word in a single document:
\begin{code}
  wordOccs :: Document -> String -> [:Int:]
  wordOccs d s = [: i | (i,s2) <- zipP [:1..lengthP d:] d
                     , s == s2 :]
\end{code}
Here we use a \emph{filter} in the array comprehension, that selects just
those pairs @(i,s2)@ for which @s==s2@.  Because this is an array comprehension,
the implied filtering is performed in data parallel.  The @(i,s2)@ pairs are chosen
from a vector of pairs, itself constructed by zipping the document with 
the vector of its indices.  The latter vector @[: 1..lengthP d :]@ is 
again analogous to the list notation @[1..n]@, which generate the list of
values between @1@ and @n@.  
As you can see, in both of these cases (filtering and enumeration)
Data Parallel Haskell tries hard to make 
parallel arrays and vectors as notationally similar as possible.

With this function in hand, it is easy to build @search@:
\begin{code}
  search :: [: Document :] -> String -> [: (Document, [:Int:]) :]
  search ds s = [: (d,is) | d <- ds
                          , let is = wordOccs d s
                          , not (nullP is) :]
\end{code}


\subsection{Prime numbers}

Let us consider the problem of computing the prime numbers up to
a fixed number @n@, using the sieve of Erathosthenes.  You may know
the cunning solution using lazy evaluation, thus:
\begin{code}
  primes :: [Int]
  primes = 2 : [x | x <- [3..]
                  , not (any (`divides` x) (smallers x))]
         where
           smallers x = takeWhile (\p -> p*p <= x) primes

  divides :: Int -> Int -> Bool
  divides a b = b `mod` a == 0
\end{code}
(In fact, this code is \emph{not} the sieve of Eratosthenes,
as Melissa O'Neill's elegant article shows \cite{oneill:sieve}, 
but it will serve our
purpose here.)
Notice that when considering a candidate prime @x@, we check that is
is not divisible by any prime smaller than the square root of @x@.
This test involves using @primes@, the very list the definition produces.

How can we do this in parallel?  In principle we want to test a whole 
batch of numbers in parallel for prime factors, so we must specify how
big the batch is:
\begin{code}
  primesUpTo :: Int -> [: Int :]
  primesUpTo 1 = [: :]
  primesUpTo 2 = [: 2 :]
  primesUpTo n = smallers +:+ 
                 [: x | x <- [: ns+1..n :]
                      , not (anyP (`divides` x) smallers) :]
    where
      ns       = intSqrt n
      smallers = primesUpTo ns
\end{code}
As in the case of @wordOccs@,
we use a boolean condition in a comprehension to filter the candidate primes. 
This time, however, computing the condition itself is a nested data-parallel computation
(as it was in @search@).  
used here to filter candidate primes @x@.  

To compute @smallers@ we make a recursive call to @primesUpTo@.  This
makes @primesUpTo@ unlike all the previous examples: the depth of
data-parallel nesting is determined \emph{dynamically}, rather than
being statically fixed to depth two.
It should be clear that the struture of the parallelism is now much more
complicated than before, and well out of the reach of mainstream flat data-parallel
systems.  But it has abundant data parallelism, and will execute with scalable
performance on a parallel processor. 

\subsection{Quicksort}

In all the examples so far the ``branching factor'' has been large.
That is, each data-parallel operations has worked on a large
collection.  What happens if the collection is much smaller?  For
example, a divide-and-conquer algorithm usually divides a problem into
a handful (perhaps only two) sub-problems, solves them, and combines
the results.  If we visualise the tree of tasks for a divide-and-conquer algorithm,
it will have a small branching factor at each node, and may be highly un-balanced.

Is this amenable to nested data parallelism?  Yes, it is.  Quicksort
is a classic divide-and-conquer algorithm, and one that we have already studied.
Here it is, expressed in Data Parallel Haskell:
\begin{code}
  qsort :: [: Double :] -> [: Double :]
  qsort xs | lengthP xs <=  1 = xs
           | otherwise        = rs!:0 +:+ eq +:+ rs!:1
           where 
             p = xs !: (lengthP xs `div` 2)
	     lt = [:x | x <- xs, x < p :]
             eq = [:x | x <- xs, x == p:]
	     gr = [:x | x <- xs, x > p :] 
	     rs = mapP qsort [: lt, gr :]
\end{code}
The crucial step here is the use of @mapP@ on a \emph{two-element} array @[: lt, gr :]@.
This says ``in data-parallel, apply @qsort@ to @lt@ and @gr@''.  The fact that there are
only two elements in the vector does not matter.
If you visualise the binary tree of sorting tasks that quicksort generates, then each
horizontal layer of the tree is done in data-parallel, even though each layer consists of many
unrelated sorting tasks.

\subsection{Barnes Hut}
All our previous examples worked on simple flat or nested
collections. Let's now have a look at an algorithm based on a 
more complex structure, in which the elements of a parallel array come from
a \emph{recursive} and \emph{user-defined} algebraic data type.

In the following, we present an implementation\footnote{Our description here is based
heavily on that in \cite{chakravarty-et-al:nepal}.} of a simple version of the
Barnes-Hut $n$-body algorithm\cite{barnes-hut:nbody}, which is a
representative of an important class of parallel algorithms covering
applications like simulation and radiocity computations.  These algorithms
consist of two main steps: first, the data is clustered in a hierarchical tree
structure; then, the data is traversed according to the hierarchical structure
computed in the first step.  In general, we have the situation that the
computations that have to be applied to data on the same level of the tree can
be executed in parallel.
Let us first have a look at  the Barnes-Hut algorithm and the
data structures that are required, before we discuss the actual implementation in parallel Haskell.




%Hierarchical $n$-body codes, such as the Barnes-Hut
%algorithm~\cite{barnes-hut:nbody}, exhibit a similar parallel structure as
%\<qsort\> and there is considerable interest in their high-performance
%implementation.  We showed that rose trees of the form
%
%\begin{haskell}
%  data Tree = Node Particle \plist{Tree}
%\end{haskell}
%%
%lend themselves to a particularly elegant nested data parallel implementation
%of the Barnes-Hut algorithm, which our program transformations can compile into
%efficient code~\cite{keller-chak:tflat}.

%\begin{verbatim}
%hier das Material von Flattening Tress - 
%davon kann man sicher was umarbeiten
%[insbesondere statt der vanilla Nesl Version,
%die mit Vektoren arbeitet, gleich die Version
%mit rekursiven Typen beschreiben]
%\end{verbatim}

An $n$-body algorithm determines the interaction between a set of particles by
computing the forces which act between each pair of particles. A precise
solution therefore requires the computations of $n^2$ forces, which is not
feasible for large numbers of particles. The Barnes-Hut algorithm minimizes
the number of force calculations by grouping particles hierarchically into
\emph{cells} according to their spatial position.  The hierarchy is
represented by a tree.  This allows approximating the accelerations induced by
a group of particles on distant particles by using the centroid of that
group's cell.  The algorithm has two phases: (1) The tree is constructed from
a particle set, and (2) the acceleration for each particle is computed in a
down-sweep over the tree. Each particle is represented by a value of type
@MassPoint@, a pair of position in the two dimensional space and mass:
\begin{code}
  type Vec       = (Double, Double) 
  type Area      = (Vec, Vec) 
  type Mass      = Double 
  type MassPoint = (Vec, Mass) 
\end{code}
We represent the tree as a node which contains the centroid and a parallel
array of subtrees:

\begin{code}
  data Tree = Node MassPoint [:Tree:]
\end{code}
Notice that a @Tree@ contains a parallel array of @Tree@.

Each iteration of @bhTree@ takes the current particle set and the area in
which the particles are located as parameters. It first splits the area into
four subareas @subAs@ of equal size. It then subdivides the particles into
four subsets according to the subarea they are located in. Then, @bhTree@ is
called recursively for each subset and subarea.  The resulting four trees are
the subtrees of the tree representing the particles of the area, and the
centroid of their roots is the centroid of the complete area. Once an area
contains only one particle, the recursion terminates.
Figure~\ref{fig:ex-area} shows such a decomposition of an area for a given set
of particles, and Figure~\ref{fig:ex-tree} displays the resulting tree
structure. 
{%\small
\begin{code}
  bhTree :: [:MassPnt:] -> Area -> Tree
  bhTree p  area = Node p [::]
  bhTree ps area =
    let
       subAs = splitArea area
       pgs   = splitParticles ps subAs
       subts = [: bhTree pg a| pg <- pgs | a <- subAs :]
       cd    = centroid [:mp | Node mp _ <- subts :] 
     in Node cd subts
\end{code}
}                         
% %  if \#ps == 1 then [Cell (ps[0], [])]            \com{cell has only one particle}
% %  else let \informal{split particles \prgm{ps} into spatial groups \prgm{pgs}}
% %           subtrees  = \{bhTree (pg): pg in pgs\};
% %           children  = \{subtree[0] : subtree in subtrees\};             \(\,(\star)\)
% %           cd        = centroid (\{mp : Cell (mp, is) in children\});
% %           sizes     = \{\#subtree : subtree in subtrees\}                \(\,(\star)\)
% %  in [Cell (cd, sizes)] ++ flatten (subtrees)  \com{whole tree as flat vector} \((\star)\)


\begin{figure}[tb]
\begin{center}
\psset{yunit=3ex,
       xunit=4.2ex,
       linewidth=.6pt}     
\begin{pspicture}(10,7)
%\psaxes(20,8)
\rput(0,-1){\psframe(8,8)}
\rput(4,-1){\psline(0,8)}
\rput(2,-1){\psline(0,8)}
\rput(6,-1){\psline(0,4)}
%
\rput(0,1){\psline(8,0)}  
\rput(0,3){\psline(8,0)}
\rput(0,5){\psline(4,0)}
%
\rput(0,0){\psline(2,0)}
\rput(1,-1){\psline(0,2)}
%
\rput(0,0){\psline(2,0)}
\rput(1,-1){\psline(0,2)}
%
\rput(2,2){\psline(2,0)}
\rput(3,1){\psline(0,2)}
%
\rput(3.5,2.5){$p_6$}
\rput(2.5,2.5){$p_7$}
\rput(5,2){$p_4$}
\rput(7,0){$p_5$}
\rput(3,4){$p_2$}
\rput(3,6){$p_3$}
\rput(1.5,0.5){$p_8$}
\rput(1.5,-0.5){$p_9$}
\rput(6,5){$p_1$}
%
\end{pspicture}
\end{center}
\caption{Hierarchical division of an area into subareas}
\label{fig:ex-area}
\end{figure}

% %
The tree computed by @bhTree@ is then used to compute the forces that act on
each particle by a function @accels@. It first splits the set of particles
into two subsets: @fMps@, which contains the particles far away (according
to a given criteria), and @cMps@, which contains those close to the centroid
stored in the root of the tree. For all particles in @fMps@, the
acceleration is approximated by computing the interaction between the particle
and the centroid. Then, @accels@ is called recursively for with @cMps@ and
each of the subtrees. The computation terminates once there are no particles
left in the set.
%
\begin{code}
  accels:: Tree -> [:MassPoint:] -> [:Vec:]
  accels _                [::] = [::]
  accels (Node cd subts)  mps  = 
    let
       (fMps, cMps) = splitMps mps
       fAcs         = [:accel  cd mp | mp <- fMps:]   
       cAcs         = [:accels t cMps| t <- subts:]   
    in combine farAcs closeAcs

  accel :: MassPoint -> MassPoint -> Vec
  -- Given two particles, the function accel computes the
  -- acceleration that one particle exerts on the other
\end{code}


\begin{figure}[tb]
\begin{center}
\psset{yunit=3ex,
       xunit=4.2ex,
       linewidth=.6pt}     
\begin{pspicture}(10,6)

\rput(-10,-3){%
  \rput(9.5,2.5){$p_6$}
  \rput(10.6,2.50){$p_7$}
  \rput(11.4,2.5){$p_8$}
  \rput(12.5,2.5){$p_9$}
                                %
  \rput(10.4,4){{$c_4$}}
  \rput(11.6,4){{$c_5$}}
  \rput(12.5,4){$p_2$}
  \rput(13.5,4){$p_3$}
  \rput(14.5,4){$p_4$}
  \rput(15.5,4){$p_5$}
                                %
  \rput(11,5.5){{$c_1$}}
  \rput(13,5.5){{$c_2$}}
  \rput(15,5.5){{$c_3$}}
  \rput(16.5,5.5){$p_1$}
                                %
  \rput(14,7){$c_0$}
                                %
  \rput(11.4,5.8){\psline(2.2,1.1)}
  \rput(16.4,5.8){\psline(-2.2,1.1)}
  \rput(13.7,6.6){\psline(-0.7,-0.7)}
  \rput(14.3,6.6){\psline(0.7,-0.7)}
                                %
  \rput(10.8,5.2){\psline(-0.4,-0.8)}
  \rput(12.8,5.2){\psline(-0.4,-0.8)}
  \rput(14.8,5.2){\psline(-0.4,-0.8)}
                                %
  \rput(11.2,5.2){\psline(0.4,-0.8)}
  \rput(13.2,5.2){\psline(0.4,-0.8)}
  \rput(15.2,5.2){\psline(0.4,-0.8)}
                                %
  \rput(10.5,3.7){\psline(0.2,-0.7)}
  \rput(11.5,3.7){\psline(-0.2,-0.7)}
                                %
  \rput(10.2,3.8){\psline(-.8,-0.8)}
  \rput(11.7,3.8){\psline(0.8,-0.8)}
}
%

\end{pspicture}
\end{center}
\caption{Example of a Barnes-Hut tree.}
\label{fig:ex-tree}
\end{figure}


% %\vspace*{-1.5ex}
% %\begin{programl}\small
% %  function accels (tree, mps) : ([Cell], [MassPnt]) -> [Vec] =
% %    if \#mps == 0 then []
% %    else let Cell (cd, chNos) = tree[0];             \com{get root}        \(\;(\star)\)
% %             \informal{split \prgm{mps} into \prgm{closeMps} and \prgm{farMps} (direct force calculation)}
% %             farAcs    = \{accel (cd, mp) : mp in farMps\};
% %             subTrees  = partition (drop (tree, 1), chNos);            \((\star)\)
% %             closeAcss = \{accels (tree, closeMps) : tree in subTrees\};
% %    in \informal{combine \prgm{farAcs} and \prgm{closeAcss}}
% %\end{programl} 
% %\vspace*{-1.5ex}
% %
% %It computes the acceleration for the mass points in \prgm{farMps} directly
% %(using the function \prgm{accel}) and recurses into the tree for those in
% %\prgm{closeMps}.  The function \prgm{drop} omits the first element of a vector
% %and \prgm{partition} forms a nested vector according to the lengths passed in
% %the second argument (it is the same as \(\Pvec\) in the next section).  The
% %type \prgm{Vec} represents `vectors' in the sense of physics.

% %The lines marked by \((\star)\) in the functions are (partially) artifacts of
% %maintaining the tree as a vector.  Figure~\ref{fig:ex-tree} depicts the
% %grouping of an exemplary particle distribution and the corresponding tree.
% %



% %
% % \begin{figure}[tb]
% % \begin{center}
% % \begin{minipage}{0.495\textwidth}
% % \epsfxsize=\textwidth\epsfbox{data.ps}
% % \end{minipage}
% % \end{center}
% % \caption{Runtime  of the Barnes-Hut NBody algorithm on the Cray
% %   T3E}
% % \label{bh-bench-runtime}
% % \end{figure}


% % \begin{figure}[tb]
% % \begin{center}
% % \begin{minipage}{0.495\textwidth}
% % \epsfxsize=\textwidth\epsfbox{speedup3.ps}
% % \end{minipage}
% % \end{center}
% % \caption{Speedup of the Barnes-Hut NBody algorithm on the Cray
% %   T3E}
% % \label{bh-bench}
% % \end{figure}



The tree is both built and traversed level by level, i.e., all nodes in one
level of the tree are processed in a single parallel step, one level after the
other. This information is important for the compiler to achieve good data
locality and load balance, because it implies that each processor should have
approximately the same number of masspoints of each level. We can see the tree
as having a sequential dimension to it, its depth, and a parallel dimension,
the breadth, neither of which can be predicted statically. The programmer
conveys this information to the compiler by the choice the data structure: By
putting all subtrees into a parallel array in the type definition, the
compiler assumes that all subtrees are going to be processed in parallel. The
depth of the tree is modelled by the recursion in the type, which is
inherently sequential.  
% The type transformation in the compilation phase,
% then, transforms the tree into a list of arrays connected by global pointers,
% where each of the arrays is distributed over the processors involved in the
% computation. 


\begin{figure}[tb]
\begin{center}
\psset{yunit=3ex,
      xunit=4.2ex,
      linewidth=.6pt}     
    \begin{pspicture}(10,9)

\rput(2,8){Processor 1}
\rput(8,8){Processor 2}

\rput(4.5,6.5){$c_0$}
\rput(5.5,6.5){$c_0$}

% level 1
\rput(3.2,5.2){\psline(0,0)(0,0.2)(3.6,0.2)(3.6,0)}
\multirput(3,4)(1,0){4}{\psframe[linecolor=bggray](1,1)}

\rput(3.5,4.5){$c_1$}
\rput(4.5,4.5){$c_2$}
\rput(5.5,4.5){$c_3$}
\rput(6.5,4.5){$p_1$}


% level 2
\rput(2.2,3.2){\psline(0,0)(0,0.2)(1.6,0.2)(1.6,0)}
\rput(4.2,3.2){\psline(0,0)(0,0.2)(1.6,0.2)(1.6,0)}
\rput(6.2,3.2){\psline(0,0)(0,0.2)(1.6,0.2)(1.6,0)}

\rput(2.5,2.5){$c_4$}
\rput(3.5,2.5){$c_5$}
\rput(4.5,2.5){$p_2$}
\rput(5.5,2.5){$p_3$}
\rput(6.5,2.5){$p_4$}
\rput(7.5,2.5){$p_5$}


\multirput(2,2)(1,0){6}{\psframe[linecolor=bggray](1,1)}

% level 3
\rput(3.2,1.2){\psline(0,0)(0,0.2)(1.6,0.2)(1.6,0)}
\rput(5.2,1.2){\psline(0,0)(0,0.2)(1.6,0.2)(1.6,0)}

\multirput(3,0)(1,0){4}{\psframe[linecolor=bggray](1,1)}
\rput(3.5,0.5){$p_5$}
\rput(4.5,0.5){$p_6$}
\rput(5.5,0.5){$p_7$}
\rput(6.5,0.5){$p_8$}


\rput(5,-0.5 ){\psline[linecolor=fggray,linestyle=dotted,linewidth=1pt](0,8)}




\rput(4,6.5){\psbezier{->}(0,0)(-1.8,-0.5)(-1.2,-1.8)}
\rput(2.6,4.5){\psbezier{->}(0,0)(-1.8,-0.5)(-1.2,-1.8)}
\rput(1.6,2.4){\psbezier{->}(0,0)(-0.6,-1.3)(1.2,-1.8)}
\end{pspicture}
\end{center}
\caption{Distribution of the values of the flattened tree}
\label{fig:ex-dist}
\end{figure}

% \subsection{Wang's Algorithm}
% \subsubsection{Solving Tridiagonal Systems of Linear Equations}

% In addition to the obvious uses of sum types, the extension of
% flattening to the full range of Haskell types
% allows a declarative type-based control of data distribution.
% Consider the operational implications for an array of arrays
% @[:[:Int:]:]@ versus an array of (sequential) lists
% @[:[Int]:]@.  On a distributed memory machine, values of the
% former will be evenly distributed over the available processing
% elements; in particular, if the subarrays vary substantially in size,
% they may be split up across processor boundaries to facilitate
% parallel operations over all elements of the nested array
% simultaneously.  In contrast, arrays of lists are optimised for
% sequential operations over the sublists; although, the sequential
% processing of all the sublists is expected to proceed in parallel.  One
% application where the distinction of parallel and sequential
% data-structures is useful is the parallel solution of tridiagonal
% systems of linear equations as proposed by Wang~\cite{Wang:1981:PMT}.

% Tridiagonal systems of linear equations are a special form of sparse
% linear systems occuring in numerous scientific applications.
% % such as special Poisson equations.
% Such system can be solved sequentially in linear time by first
% eliminating the elements of the lower diagonal by a top-down traversal,
% and then eliminating the upper diagonal by traversing the matrix from bottom to
% top.  Unfortunately, in each step a pivot row is needed that is computed
% just in the step before, so the algorithm is completely sequential.

% In the parallel solution proposed by Wang,
% the matrix is subdivided into blocks of consecutive rows, which
% are then processed simultaneously.
% The algorithm runs in 
% three phases. 
% %parallel and one sequential phase. 
% First, all rows of a block are traversed top-down and then bottom-up to eliminate the
% lower and upper diagonal, respectively. However, since the first row
% in each but the first block still contains the lower diagonal
% element, a vertical chain of fill-in elements appears in this column. 
% As the matrix is symmetric, a chain of fill-ins also occurs on the right in all but the last
% block in the bottom-up traversal. 
% The non-zero elements of the matrix after the first phase are shown in Fig.~\ref{fig:wang}.
% % % \begin{figure}[htb]
% % % \begin{center}
% % %     \leavevmode
% % %     \epsfig{file=wang.eps,scale=0.4} 
% % % \caption{Situation with 3 blocks after first parallel phase in Wang's algorithm}
% % %   \label{fig:wang}
% % % \end{center}
% % % \end{figure}

% To diagonalise the matrix, the left and right chains of fill-ins must be
% eliminated. 

% The first block's last row contains non-zeros suitable 
% for elimination of all left fill-ins in the second block. 
% Once the left chain element of the second block's last row has been
% eliminated, this updated row can be used as a pivot for the 
% elimination of the left fill-in chain in the third block etc.
% Thus, a pipelining phase is necessary over all blocks to
% propagate suitable pivot rows for the elimination of the left chains of
% fill-ins. Analogously, pivots can be propagated upwards starting with
% the last block to eliminate the right chains of extra non-zeros. 

% In each block, once the pivot row from the preceding block is available, the
% fill-in elements may be eliminated in any order. There are no sequential
% inner-block dependencies. However, as described above, there is a sequential
% dependency among the blocks.  Elimination of the left chain can start only
% after the pivot row from the previous block is available, but this is the case
% only after the left fill-in of the previous block's last row has been
% eliminated already.  Thus, it is important that during pipelining, only the
% first and last rows of each block are touched, because eliminating all
% fill-ins first before propagating pivots to the next block would mean a
% completely sequential traversal of the matrix.

% After the pipelining phase, there are pivot rows for each block that can
% be used to eliminate both the left and the right chains of fill-ins. 
% Like in the first phase, all blocks can be processed in parallel.
% Again, one top-down and one bottom-up traversal are necessary to
% obtain the desired diagonal structure. 

% \subsubsection{Encoding Wang's Algorithm in parallel Haskell}

% In parallel Haskell, we model an equation with a tuple-type @TriRow@ containing the
% three diagonal elements, the two potential chain elements, and the right-hand
% side.
% \begin{code}
% type TriRow    = (Float,Float,Float,Float,Float,Float)
%                --left, lower, main, upper, right, rhs
% \end{code}
% A row block is a list of rows, i.e., of type
% @[TriRow]@. The whole matrix is a parallel array of row blocks,
% abbreviated by the type @Matrix@. 
% \begin{code}
% type Matrix = [:[TriRow]:] -- a parallel array of lists of rows
% \end{code}

% The following encodes the top-level function of Wang's algorithm.

% \begin{code}
%   solve :: Matrix -> [:[Float]:]
%   solve m = 
%     let
%       res  = [:elimLowerUpper x | x <- m:] -- Phase 1
%       frv  = [:f | (_,f,_) <- res:] 
%       lrv  = [:l | (_,_,l) <- res:] 
%       rowv = [:r | (r,_,_) <- res:]
%       (fpl,lpl) = pipeline (fromPArray frv)
%                            (fromPArray lrv) -- Phase 2
%       (fpv,lpv) = (toPArray fpl, toPArray lpl)
%       dm        = [:elimLeftRight r fp lp | r  <- rowv | 
%                                             fp <- fpv  | lp <- lpv :] -- Phase 3
%      in  mapP (map (\(TriRow _ _ maine _ _ rhs) -> rhs/maine)
% \end{code}
% %\vspace*{-1.5ex}
% % %
% % The functions @elimLowerUpper@ and @elimLeftRight@ are ordinary, recursive
% % list-traversals, eliminating elements on each row both in the descending and
% % ascending phase of re\-cur\-sion---we omit the details of their definition
% % here, as they do not use parallelism. However, these traversals are executed
% % in parallel for all blocks. The function @elimLowerUpper@ is of type
% % @[TriRow] \to ([TriRow],TriRow, TriRow)@. It returns the updated row block
% % plus the two rows needed for the pipelining phase.  As the pipelining is
% % sequential, lists are used and so the arrays with the first and last pivot
% % rows are converted by the primitive @pArrayToList@.  The function
% % @pipeline@ is again an ordinary list traversal, realizing the desired pivot
% % generation and propagation.  The lists of new pivot rows are transformed into
% % parallel arrays using @listTopArray@, so that the third phase can work in
% % parallel on all blocks to eliminate the fill-in values.


\subsection{A performance model}

\simon{This seciton feels a bit wooly, but I think we owe the students something.}

One of the main advantages of the data parallel programming model is that it comes with
a \emph{performance model} that lets us make reasonable predictions about the behaviour of
the program on a parallel machine, including its \emph{scalability} -- that is, how performance
changes as we add procesors.  So what is this performance model?

First, we must make explicit something we have glossed over thus far:
data-parallel arrays are strict.  More precisely, if any element of a
parallel diverges, then all elements diverge\footnote{%
What if the elements are pairs?  See Leshchinskiy's thesis
for the details \cite{leshchinskiy:thesis}.}.  This makes sense,
because if we demand any element of a parallel array then we must
compute them all in data parallel; and if that computation diverges we
are justified in not returning any of them.  The same constraint means
that we can represent parallel arrays very efficiently.  For example,
an array of floats, @[:Float:]@, is represnted by a contiguous array
of unboxed floating-point numbers.  There are no pointers, and
iterating over the array has excellent spatial locality.

In reasoning about performance, Blelloch \cite{blelloch:cacm} characterises the \emph{work} and \emph{depth} of
the program:
\begin{itemize}
\item The \emph{work}, $W$, of the program is the time it would take to execute on a single processor.
\item The \emph{depth}, $D$, of the program is the time it woudl take to execute on an 
infinite number processors, under the assumption that the additional processors leap into action
when (but only when) a @mapP@, or other data-parallel primitive, is executed.
\end{itemize}
If you think of the unrolled data-flow diagram for the program, the work is the number of
nodes in the data-flow diagram, while the depth is the longest path from input to output.

Of course, we do not have an infinite number of processors.  Suppose instead tha we have $P$ processors.  
Then if everything worked perfectly, the work work be precisely evenly balanced across the processors
and the execution time $T$ would be $W/P$.  That will not happen if the depth $D$ is very large.  So in
fact, we have
$$ W/P \leq T \leq W/P + L*D $$
where $L$ is a constant that grows with the latency of communication in the machine.
Even this is a wild approximation, because it takes no account of bandwidth
limitations.  For example, between each of the recursive calls in the Quicksort example there must
be some data movement to bring together the elements less than, equal to, and greater than the pivot.
Nevertheless, if the network bandwidth of the parallel machine is high (and on serious multiprocessors 
it usually is) the model gives a resonable approximation.

\begin{figure}
$$
\begin{array}{rcll}
\depth{k} & = & 0 & \mbox{where $k$ is a constant}\\
\depth{x} & = & 0 & \mbox{where $x$ is a variable}\\
\depth{e_1 + e_2} & = & 1 + \depth{e_1} + \depth{e_2} \\
\\
\depth{@if@ ~e_1~@then@~e_2~@else@~e_3} & = & \depth{e_1} + \depth{e_2} & \mbox{if}~e_1 = @True@ \\
  & = & \depth{e_1} + \depth{e_3} & \mbox{if}~e_1 = @False@ \\
\depth{@let@ ~x@=@ e~@in@~b} & = & \depth{b[e/x]} \\
\\
\depth{e_1 ~@+:+@~e_2} & = & 1 + \depth{e_1} + \depth{e_2} \\
\depth{@concatP@~e} & = & 1 + \depth{e} \\
\depth{@mapP@~f~e} & = & 1 + \depth{e} + \maxx{x \in e}\depth{f~x} \\
\depth{@filterP@~f~e} & = & 1 + \depth{e} + \depth{f} \\
\\
\depth{@sumP@~e}   & = & 1 + \depth{e}   + log(length(e))
\end{array}
$$
\caption{Depth model for closed expressions} \label{fig:depth}
\end{figure}

How can we compute work and depth?  
It is much easier to reason about the work of a program in a strict setting than in a lazy one, because
all sub-expressions are evaluated.  This is why the performance model of the data-parallel part of
DPH is more tractable than for Haskell itself.  

The computation of depth is where we take account of data parallelism.  
Figure~\ref{fig:depth} shows the equations for calculating the depth of a closed
expression $e$, where $\depth{e}$ means ``the depth of $e$''.
These equations embody the following ideas:
\begin{itemize}
\item By default execution is sequential.  Hence, the depth of an addition is the sum of the
depths of its arguments.
\item The parallel primitive @mapP@, and its relatives such as @filterP@,
can take advantage of parallelism, so the depth is the worst depth encountered for any element.
\item The parallel reduction primitive @sumP@, and its relatives, 
take time logarithmic in the length of the array.
\end{itemize}
The rule for @mapP@ dirctly embodies the idea that nested data parallelism is flattened.  For example,
suppose $e :: @[:[:Float:]:]@$.  Then, applying the rules we see that 
$$\begin{array}{rcl}
\depth{@mapP@~f~(@concatP@~e} & = & 1 + \depth{@concatP@~e} + \maxx{x \in {\scriptsize @concatP@\, e}}\depth{f~x} \\[2mm]
  & = & 1 + 1 + \depth{e} + \maxx{x \in {\scriptsize @concatP@}\, e}\depth{f~x} \\[2mm]
  & = & 2 + \depth{e} + \maxx{xs \in e}~\maxx{x \in xs}\depth{f~x} \\[3mm]
\depth{@mapP@~(@mapP@~f)~e} & = & 1 + \depth{e} + \maxx{xs \in e}\depth{@mapP@~f~xs} \\[2mm]
 & = & 1 + \depth{e} + 1 + \maxx{xs \in e}~\maxx{x \in xs}\depth{f~x} \\[2mm]
  & = & 2 + \depth{e} + \maxx{xs \in e}~\maxx{x \in xs}\depth{f~x}
\end{array}$$
Notice that although the second case is a \emph{nested} data-parallel computation, it has the same
depth expression as the first: the data-parallel nesting is flattened.

These calculations are obviously very approximate, certainly so far as constant factors are concerned.
For example, in the inequality for execution time, $$W/P \leq T \leq W/P + L*D $$
we do not know the value of the latency-related constant $L$.  
However, what we are primarily looking for is the \emph{Asymptotic Scalability} (AS) property: 
\begin{quote}
A program has the Asymptotic Scalability property if
$D$ grows asymptotically more slowly than $W$, as the size of the 
problem increases.
\end{quote}
If this is so then, for a sufficiently large problem and assuming sufficient network bandwidth,
performance should scale linearly with the number of processors.

For example, the functions @sumSq@ and @search@ both have constant depth, 
so both have the AS property, and (assuming sufficient bandwidth) performance should scale linearly with 
the number of processors after some fairly low threshold.

For Quicksort, an inductive argument shows that the depth is logarithmic in the size of the array, assuming
the pivot is not badly chosen.  So $W = O(n log n)$ and $D= O(log n)$, and Quicksort has the AS property.

For computing primes, the depth is smaller: $D = O(log log n)$.  Why?  Because at every step we take
the square root of $n$, so that at depth $d$ we have $n = 2^{2^d}$.  Almost all the work is done at the
top level.  The work at each level involves
comparing all the numbers between $\sqrt{n}$ and $n$ with each prime smaller than $\sqrt{n}$.  There
are approximately $\sqrt{n}/log n$ primes smaller than $\sqrt{n}$, so the total work is roughly
$W = O(n^{3/2}/log n)$.  So again we have the AS property.

Leshchinskiy \emph{et al} \cite{leshchinskiy-et-al:costing-ndp} give further details of the cost model.

\subsection{How it works}

NESL's key insight is that it is possible to transform a program that uses \emph{nested} data-parallelism into
one that uses only \emph{flat} data parallelism.  While this little miracle happens behind the scenes,
it is instructive to have some idea how it works, just as a car driver may find some knowledge of 
internal combustion engines even if he is not a skilled mechanic.
The description here is necessarily brief, but the reader may find a slightly more detailed
overview in \cite{chakravarty-et-al:ndp-status}, and in the papers cited there.

We call the nested-to-flat tansformation the \emph{vectorisation} transform.  It has two parts:
\begin{itemize}
\item Transform the \emph{data} so that all parallel arrays contain only primitive, flat data, such
as @Int@, @Float@, @Double@.
\item Transform the \emph{code} to manipulate this flat data.
\end{itemize}
To begin with, let us focus on the first of these topics.  We may consider it as the driving force,
because nesting of data-parallel operations is often driven by nested data structures.

\subsubsection{Transforming the data}

As we have already discussed, a parallel array of @Float@ is represented by a contiguous array
of honest-to-goodness IEEE floating point numbers; and similarly for @Int@ and @Double@.
It is as if we could define the parallel-array type by cases, thus:
\begin{code}
  data instance [: Int    :] = PI Int ByteArray
  data instance [: Float  :] = PF Int ByteArray
  data instance [: Double :] = PD Int ByteArray
\end{code}
In each case the @Int@ field is the size of the array.
These @data@ declarations are unusual because they are \emph{non-parametric}: the representation
of an array depends on the type of the elements\footnote{
None of this is visible to the programmer, but the @data instance@ notation is in fact
available to the programmer in recent versions of GHC 
\cite{schrijvers-et-al:open-type-functions,chakravarty:assoc-syn}.
Why?  Because GHC has a typed intermediate language so we needed to figure out how to give
a \emph{typed} account of the vectorisation transformation, and once that is done it seems
natural to offer it to the programmer.  Furthermore, much of the low-level 
support code for nested data parallelism is itself written in Haskell, and operates directly on the
post-vectorisation array representation.
}.

Matters become even more interesting when we how to represent a parallel array of pairs.  We must
not represent it as a vector of pointers to heap-allocated pairs, scattered randomly around the 
address space.  We get much better locality if we intead represented it as a \emph{pair of arrays} thus:
\begin{code}
  data instance [: (a,b) :] = PP [:a:] [:b:]
\end{code}
What about a parallel array of parallel arrays?  Again, we must avoid a vector of pointers.
Instead, the natural representation is obtained by literally concatenating the (representation of)
the sub-vectors into one giant vector, together with a vector of indices to indicate where each of 
the sub-vectors begins.  \simon{Satnam, diagram needed.}
\begin{code}
  data instance [: [:a:] :] = PA [:Int:] [:a:]
\end{code}
By way of example, recall the data types for sparse matrices:
\begin{code}
  type SparseMatrix = [: SparseVector :]
  type SparseVector = [: (Int, Float) :]
\end{code}
Now consider this tiny matrix, consisting of two short documents:
\begin{code}
  m :: SparseMatrix
  m = [: [:(1,2.0), (7,1.9):], [:(3,3.0):] :]
\end{code}
This would be represented as follows:
\begin{code}
  PA [:0,2:] (PP [:1,   7,   3  :]
                 [:1.0, 1.9, 3.0:])
\end{code}
The array as the leaves are themeselves represented as byte arrays:
\begin{code}
  PA (PI 2 #<0x0,0x2>)
     (PP (PI 3 #<0x0,    0x7,     0x3>)
         (PF 3 #<0x9383, 0x92818, 0x91813>))
\end{code}
Here we have invented a fanciful notation for literal @ByteArray@s (not supported
by GHC, let alone Haskell) to stress the fact that in the end everything boils down
to literal bytes.  (The hexadecimal encodings of floating point numbers are also 
made up because the real ones have many digits!)

We have not discussed how to represent arrays of sum types (such as @Bool@, @Maybe@, or lists),
nor of function types --- see \cite{chakravarty-keller:more-types} and \cite{leshchinskiy:thesis}
respectively.

\subsubsection{Vectorising the code}

As you can see, data structures are transformed quite radically by the vectorisation transform,
and it follows that the code must be equally radically transformed.  Space precludes proper treatment
here; a good starting point is Keller's thesis \cite{keller:thesis}.

A data-parallel program has many array-valued sub-expressions.  For example, in @sumSq@ we see
\begin{code}
  sumSq a = sumP [: x*x | x <- a :]
\end{code}
However, if @a@ is a big array, it would be silly to compute a new, equally big array of squares,
only to immediately consume it with @sumP@. It would be much better for each processor to zip down its
chunk of @a@, adding the square of each element into a running total, and for each processor's total to be
combined.  

The elimination of intermediate arrays is called \emph{fusion} and is
crucial to improve the constant factor of Data Parallel Haskell.  It
turns out that vectorisation introduces many \emph{more} intermediate
arrays, which makes fusion even more important.  These constant
factors are extremely important in practice: if there is a slow-down of a factor of 50
relative to C, then even if you get linear speedup by adding
processors, Data Parallel Haskell is unlikely to become popular.

\subsection{Running Data Parallel Haskell}

GHC 6.6 and 6.8 come with support for Data Parallel Haskell syntax, and a \emph{purely sequential}
implementation of the operations.  So you can readily try out all of the examples in this paper,
and ones of your own devising thus:
\begin{itemize}
\item Use @ghc@ or @ghci@ version 6.6.x or 6.8.x.
\item Use flags @-fparr@ and @-XParallelListComp@. 
\item Import module @GHC.PArr@.
\end{itemize}
Some support for genuinely-parallel Data Parallel Haskell, including the all-important vectorisation 
transformation, will be in GHC 6.10 (planned release: autumn 2008).
It is not yet clear just how complete the support will be at that time.
At the time of writing, for example, type classes are not vectorised, and neither
are lists. Furthermore, in a full implementation we will want need support for 
partial vectorisation \cite{chakravarty-et-al:partial-vect}, among other things.

As a result, although all the examples in this paper should
work when run sequentially, they may not all vectorise as written, even in GHC 6.10.

A good source of working code is in the Darcs repository \url{http://darcs.haskell.org/packages/ndp},
whose sub-directory \url{examples/} contains many executable examples.

\subsection{Further reading}

\simon{Manuel, is this a reasonable set of pointers?}

Blelloch and Sabot originated the idea of compiling nested data parallelism into flat data parallelism
\cite{blelloch-sabot:ndp}, but an easier starting point is probably Blelloch subsequence CACM paper 
``Programming parallel algorithms'' \cite{blelloch:cacm}, and the NESL 
language manual \cite{blelloch:nesl-manual}.

Keller's thesis \cite{keller:thesis} formalised an intermediate language that models the central aspects
of data parallelism, and formalised the key vectorisation transformation.  She also studied array
\emph{fusion}, to eliminate unnecessary intermediate arrays.  Leshchinskiy's thesis \cite{leshchinskiy:thesis}
extended this work to cover higher order languages (\cite{leshchinskiy-et-al:higher-order-flattening} gives a 
paper-sized summary), while Chakravarty and Keller explain a further generalisation to handle
user-defined algebraic data tyeps \cite{chakravarty-keller:more-types}.

Data Parallel Haskell is an ongoing research project \cite{chakravarty-et-al:ndp-status}.
The Manticore projet at Chicago shares similar goals \cite{reppy-et-al:manticore}.

\bibliographystyle{alpha}
\bibliography{ndp-tutorial}

\end{document}