\documentclass{fsttcs}
\usepackage{amsmath}
\usepackage{version}
\usepackage{url}
\usepackage{graphicx}
\usepackage{epsfig}
\usepackage{pifont}
\usepackage{pgf,pgfarrows,pgfnodes}
\usepackage{stmaryrd}
\usepackage{code}
\usepackage{verbatim}
\usepackage{ifthen}
\usepackage{amsfonts,amssymb}
\usepackage{url}
\usepackage{alltt}

% version control
%
\includeversion{DRAFT}
%\excludeversion{DRAFT}

\def\Version{\relax}
\begin{DRAFT}%
  {
    \gdef\Version{%
      \\
      \textbf{--- DRAFT ---}\\[1ex]
      \ignorespaces}
  }
\end{DRAFT}


% Put figures in boxes
\usepackage{float}
\floatstyle{boxed}
\restylefloat{figure}

% TODO's
%
\newcommand{\authornote}[1]{{\emph{\textbf{#1}}}}
%\newcommand{\authornote}[1]{}
%\newcommand{\authornote}[1]{\emph{[#1]}}
\newcommand{\TODO}[1]{\authornote{TODO: #1}}
\newcommand{\todo}[1]{\authornote{TODO: #1}}
\newcommand{\simon}[1]{\authornote{[Simon: #1]}}
\newcommand{\manuel}[1]{\authornote{[Manuel: #1]}}
\newcommand{\gabi}[1]{\authornote{[Gabi: #1]}}

% Double semantic brackets
\newcommand{\lb}{[\![}
\newcommand{\rb}{]\!]}
\newcommand{\db}[1]{\lb #1 \rb}

%
%  pgf stuff
%  ---------
\newcommand{\mpgfrect}[4]{\pgfrect[stroke]{\pgfpoint{#1}{#2}}{\pgfpoint{#3}{#4}}}
\newcommand{\mpgfline}[4]{\pgfline{\pgfpoint{#1}{#2}}{\pgfpoint{#3}{#4}}}

\newcommand{\mpgfbox}[2]{\mpgfrect{#1}{#2}{0.2cm}{0.2cm}}    

% end pgf

\renewcommand{\textfraction}{0.2}
\renewcommand{\topfraction}{0.9}
\renewcommand{\dbltopfraction}{0.9}
\renewcommand{\floatpagefraction}{0.9}
\renewcommand{\dblfloatpagefraction}{0.9}

% Misc macros
%
\newcommand{\tup}[1]{\overline{#1}}

% Program font stuff
\newcommand{\ccode}[1]{\texttt{#1}}
\newcommand{\hskwd}[1]{\texttt{#1}}
\newcommand{\hslam}[1]{\ensuremath{\lambda}}
\newcommand{\capp}{\texttt{\$\!:}}  % ditto
\newcommand{\cappP}{\mathbin{{\$}{:}^\uparrow}}  % ditto

% Notation for type-indexed parallel arrays
\newcommand{\pacls}{PAElem}         % Parallel array class
\newcommand{\pacon}{PA}
\newcommand{\pa}[1]{\pacon\;#1}       % PArr in *vectorised* code
\newcommand{\paB}[1]{(\pa{#1})}   % Bracketed PArr in *vectorised* code

\newcommand{\patype}[1]{\mathsf{patype}\langle#1\rangle}  % Not sure about this

\newcommand{\inform}[1]{\ensuremath{\langle\textrm{#1}\rangle}}
\newcommand{\var}[1]{\ensuremath{\mathit{#1}}}
\newcommand{\sub}[1]{\ensuremath{_{\mathnormal{#1}}}}
\newcommand{\LV}[1]{#1\ensuremath{_{\mathrm{L}}}}
\newcommand{\VV}[1]{#1\ensuremath{_{\mathrm{V}}}}
\newcommand{\SE}[1]{#1\ensuremath{_{\mathrm{S}}}}
\newcommand{\Arr}[1]{#1$_{\code{\pacon}}$}

\newcommand{\meta}[1]{\ensuremath{\llbracket\text{#1}\rrbracket}}
\newcommand{\TT}[1]{\ensuremath{\mathcal{V}_t\llbracket #1 \rrbracket}}
\newcommand{\TTL}[1]{\ensuremath{\mathcal{L}_t\llbracket #1 \rrbracket}}
\newcommand{\VT}[1]{\ensuremath{\mathcal{V}\left\llbracket #1 \right\rrbracket}}
\newcommand{\LT}[2]{\ensuremath{\mathcal{L}\left\llbracket #2 \right\rrbracket\,#1}}
\newcommand{\LTprime}[2]{\ensuremath{\mathcal{L}'\left\llbracket #2 \right\rrbracket\,#1}}

\pagenumbers{383}{414}

\begin{document}


\title{Harnessing the Multicores: Nested Data Parallelism in Haskell}
\author{%
  Simon Peyton Jones\inst{1},
  Roman Leshchinskiy\inst{2},
  Gabriele Keller\inst{2},
  Manuel M. T. Chakravarty\inst{2}}
\institute{1}{%
   Microsoft Research Ltd,
   Cambridge, England,
   \email{simonpj@microsoft.com}
 }
\institute{2}{%
   Programming Languages and Systems,
   School of Computer Science and Engineering,
   University of New South Wales,
   \email{\{rl,keller,chak\}@cse.unsw.edu.au}
 }
\runningtitle{Harnessing the Multicores}
\runningauthors{Peyton Jones, Leshchinskiy, Keller, Chakravarty}

\makeatactive

\begin{abstract}
If you want to program a parallel computer, a purely functional
language like Haskell is a promising starting point.  Since the
language is pure, it is by-default safe for parallel evaluation,
whereas imperative languages are by-default unsafe.  But that doesn't
make it easy!  Indeed it has proved quite difficult to get robust,
scalable performance increases through parallel functional
programming, especially as the number of processors increases.

A particularly promising and well-studied approach to employing large
numbers of processors is data parallelism.  Blelloch's
pioneering work on NESL showed that it was possible to combine a
rather flexible programming model (nested data parallelism) with a
fast, scalable execution model (flat data parallelism).  In this paper
we describe Data Parallel Haskell, which embodies nested data
parallelism in a modern, general-purpose language, implemented in a
state-of-the-art compiler, GHC.  We focus particularly on the
vectorisation transformation, which transforms nested to flat data
parallelism.
\end{abstract}


\section{Introduction}
\label{sec:intro}

Computers are no longer getting faster; instead, we will be offered
computers containing more and more CPUs, each of which is no faster
than the previous generation.  As the number of CPUs increases, it
becomes more and more difficult for a programmer to deal with the interactions
of large numbers of threads.  Moreover, the physical limitations of bus
bandwidth will mean that memory access times will be increasingly non-uniform
(even if the address space is shared), and locality of reference will be
increasingly important.

In the world of massively-parallel computing with strong locality
requirements there is already a well-established, demonstrably
successful brand leader, namely \emph{data parallelism}.  In a
data-parallel computation one performs the \emph{same} computation on
a large collection of \emph{differing} data values.  Well-known
examples of data-parallel programming environments are High
Performance Fortran (HPF) \cite{HPF2}, the collective operations of the
Message Passing Interface (MPI) \cite{mpi2}, NVIDIA's Compute Unified
Device Architecture (CUDA) API for 
graphics processors \cite{CUDA},
and Google's map/reduce framework \cite{map-reduce}.

All these systems support only \emph{flat} data parallelism, in which the
computation that is performed on each data element must itself be (a)
sequential and (b) of a similar execution time to the computation on the other
data elements.  In practice, this severely limits the applications of
data-parallel computing, especially for sparse or irregular
problems~\cite{prins-etal:irr-fortran}.  Thus motivated, Blelloch and
Sabot developed the idea of \emph{nested} data parallelism in the early 90's,
and embodied it in their language
NESL~\cite{blelloch-sabot:compiling-collection}.

NESL was a seminal breakthrough but, fifteen years later it remains
largely un-exploited.  Our goal is to adopt the key insights of NESL,
embody them in a modern, widely-used functional programming language,
namely Haskell, and implement them in a state-of-the-art Haskell compiler
(GHC).  The resulting system, Data Parallel Haskell, will make 
nested data parallelism available to real users.

Doing so is not straightforward. NESL a first-order language, has very few
data types, was focused entirely on nested data parallelism,
and its implementation is an interpreter.  Haskell is a
higher-order language with an extremely rich type system; it already 
includes several other sorts of parallel execution; and its
implementation is a compiler.

This paper makes two main contributions:
\begin{itemize}
\item We give a tutorial, programmer's-eye view of what programming
in Data Parallel Haskell is like.  Rather than a series of tiny 
examples, we give a serious application that is very hard to 
fully parallelise in a flat data-parallel setting, namely the Barnes-Hut
algorithm for $N$-body simulation.

\item We give a detailed tutorial overview of the key vectorisation 
transformation.  There are two major innovations over NESL: one is the
non-parametric representation of arrays (Section
~\ref{sec:non-parametric}) and one is the treatment of first-class
functional values (Section ~\ref{sec:functions}).
\end{itemize}

All the technical innovations in this paper have appeared, piecemeal,
in our earlier publications.  Our hope, however, is that this paper
draws together a somewhat-complex set of technical strands into
a comprehensible whole.



\section{The programmer's view on DPH}
% introducing the language using barnes hut as example
% spend more time on explaining more complex data structures (ie tree) and the sequential &
% parallel nature of those structures
\label{sec:prgView}


\begin{figure}
\begin{code}
(!:)       :: [:a:] -> Int -> a
sliceP     :: [:a:] -> (Int,Int) -> [:a:]
replicateP :: Int -> a -> [:a:]
mapP       :: (a->b) -> [:a:] -> [:b:]
zipP       :: [:a:] -> [:b:] -> [:(a,b):]
zipWithP   :: (a->b->c) -> [:a:] -> [:b:] -> [:c:]
filterP    :: (a->Bool) -> [:a:] -> [:a:]

concatP    :: [:[:a:]:] -> [:a:]
concatMapP :: (a -> [:b:]) -> [:a:] -> [:b:]
unconcatP  :: [:[:a:]:] -> [:b:] -> [:[:b:]:]
transposeP :: [:[:a:]:] -> [:[:a:]:] 
expandP    :: [:[:a:]:] -> [:b:] -> [:b:]

combineP   :: [:Bool:] -> [:a:] -> [:a:] -> [:a:]
splitP     :: [:Bool:] -> [:a:] -> ([:a:], [:a:])
\end{code}
\caption{Type signatures for parallel array operations} \label{fig:parr-ops}
\end{figure}
We begin by describing Data Parallel Haskell (DPH) purely from the point
of view of the programmer, illustrating the description with a 
non-trivial example, the Barnes-Hut algorithm~\cite{barnes-hut:nbody}.  
GHC supports other forms of concurrency besides data parallelism,
but we focus here exclusively on the latter. Singh~\cite{singh-pj:afp08}~gives a tutorial
covering a broader scope, including semi-implicit parallelism (@par@),
explicit threads, transactional memory, as well as Data Parallel Haskell.

DPH is simply Haskell with the following extra features:
\begin{itemize}
\item A type of \emph{parallel arrays}, denoted @[:e:]@ for arrays of type @e@.
These arrays are indexed by values of type @Int@.  From a semantic point
of view an array @[:a:]@ is very similar to a list @[a]@ -- the difference is
in the execution pragmatics.
An array can contain elements of any type, including arrays and functions.
\item A large number of \emph{parallel operations} that operate collectively on
entire arrays. As far as possible,
these operations have the same names as Haskell's standard list
functions, but with the suffix @P@ added---i.e., @mapP@, @filterP@,
@unzipP@, and so forth.  Figure~\ref{fig:parr-ops} lists the operations
that we will use in this paper.
\item Syntactic sugar, called \emph{parallel array
comprehensions,} which are similar to list comprehensions but operate on
parallel arrays.
\end{itemize}
In addition to the parallel evaluation semantics, lists and parallel
arrays also differ with respect to strictness: more precisely, demand
for \emph{any} element of a parallel array results in the evaluation of \emph{all}
elements.

\subsection{$N$-Body Barnes-Hut Simulation Algorithm}

We will demonstrate the use of DPH features using the Barnes-Hut
$n$-body simulation algorithm as an example.
We discuss the algorithm in some detail because it is a particularly
striking example of the power of nested data parallelism, and of the
utility of user-defined data types in data-parallel programs.  
We will, for the sake of clarity, restrict ourselves to two dimensions
and neglect complications such as bodies that are very close to each other.

An $n$-body simulation computes the motion under gravitational forces
of $n$ bodies, or \emph{particles}.  A naive solution is to compute the force
between every pair of particles which requires $n^2$ calculations in each time
step. The Barnes-Hut algorithm reduces the work
complexity to the order of $n\log n$ interactions by grouping together
particles which are close to each other and calculating the centre of
gravity, or \emph{centroid} of the cluster. The centroids are then used
to approximate the effect the particles have on other particles which
are sufficiently far away. The stricter we are in determining what
exactly constitutes ``sufficiently far away'', the more precise the
final result is, and the algorithm can be parametrised accordingly.

The first phase of the algorithm determines
the hierarchical grouping of the particles, computes the centroids of the
clusters, and stores the result in a tree structure. To be more
precise, the area is split into four subareas of equal size, the
particles are grouped according to the subarea they are located in. We
repeat this step for each subarea, and terminate if an area contains either none or
only a single particle. 
Figure~\ref{fig:quadtree-graph} illustrates the tree construction process
for particles $p_1 \ldots p_9$. In the first
iteration, the particles are split into four groups,
depending on which quadrant they are located in. The upper right
quadrant already contains only a single particle, so it isn't divided
up any further. Both the upper left and the lower right quadrant
require only one more iteration, the lower left two iterations.


% This figure needs no box
% \floatstyle{plain}
% \restylefloat{figure}
% \begin{figure*}[tp]
%   \centering
%   \includegraphics[width=0.5\textwidth]{quadtree-graph.pdf}
%   \caption{Subdivision of area}
%   \label{fig:quadtree-graph}
% \end{figure*}
% \floatstyle{boxed}
% \restylefloat{figure}


Figure~\ref{fig:quadtree} shows the resulting hierarchical tree structure:
the root node contains the centroid of all particles and four
subtrees (since all subareas contain at least one particle).  Each of the
subtrees contains the centroid of the corresponding subarea -- which
is the particle itself in case of a singular particle.

% \begin{figure*}[tp]
%   \centering
%   \includegraphics[width=0.5\textwidth]{quadtree-tree.pdf}
%   \caption{Rosetree}
%   \label{fig:quadtree}
% \end{figure*}



The second phase of the algorithm now calculates the forces that affect each
particle $p$, by traversing the tree from the root downwards: for every
subtree, if the particle $p$ is sufficiently far away from the centroid stored
in the root of that subtree, the force on $p$ is calculated using this
centroid without looking at the rest of the tree.  Otherwise, we add up the
forces on $p$ from the subtrees of the current root -- and so on recursively.

\subsection{Encoding Barnes-Hut in DPH} \label{sec:bh-dph}

Since the only way to express parallelism
in DPH is to apply collective operations to parallel arrays, we need
to store all data that we want to process in parallel in such an
array.  For instance, the function @oneStep@, which computes one step in
the simulation, takes a parallel array of particles as arguments, and returns
an array of the same length, with the position and velocity of each
particle adjusted according to the gravitational forces:%
\linebreak[4]
\par{\small
\begin{code}
  -- Compute one step of the n-body simulation
  oneStep :: [:Particle:] -> [:Particle:]
  oneStep particles = moveParticles particles forces
    where
      tree   = buildTree initialArea particles
      forces = calcForces (lengthOf initialArea) tree particles

  buildTree     :: Area -> [:Particle:] -> Tree
  calcForces    :: Float -> Tree -> [:Particle:] -> [:Force:]
  moveParticles :: [:Particle:] -> [:Force:] -> [:Particle:]
  lengthOf      :: Area -> Float
\end{code}}
\noindent The function @oneStep@, as discussed before, is comprised of three
data parallel phases. First, @buildTree@ decomposes the particles
into sub-areas, returning the resulting @Tree@. This tree is then used by
@calcForces@ to compute the
forces on the particles, returning a new array of forces with one element for
each particle.
Finally, @moveParticles@ uses these forces to adjust the positions and
velocities the particles.

The data types involved in the computation are defined exactly as the would be
in regular Haskell.
For example, a @Particle@ is a record of its mass, its location, and
its velocity:
\par{\small
\begin{code}
  type Vector   = (Float, Float)
  type Area     = (Vector, Vector)
  type Force    = Vector
  type Velocity = Vector
  type Location = Vector

  data Particle = Particle { mass     :: Float 
                           , location :: Location
                           , velocity :: Velocity}
\end{code}}
\noindent Some functions are conveniently defined using the parallel array
counterparts of ordinary list processing functions (see
Figure~\ref{fig:parr-ops}).  For example, we can define @moveParticles@ like
this:
\par{\small
\begin{code}
  moveParticles :: [:Particle:] -> [:Force:] -> [:Particle:]
  moveParticles ps fs = zipWithP moveParticle ps fs

  moveParticle :: Particle -> Force -> Particle
  moveParticle (Particle { mass     = m
                         , location = loc
                         , velocity = vel })
               force
    = Particle { mass     = m
               , location = loc + vel   * timeStep
               , velocity = vel + accel * timeStep }
    where
      accel = force / m
\end{code}}
\noindent Now we turn our attention to the @Tree@ data type and its construction.
When building and traversing a @Tree@, we want to process its
sub-trees in parallel, and so we must use a parallel array
for the children:
\begin{code}
  data Tree = Node Mass Location [:Tree:]
              -- Rose tree for spatial decomposition
\end{code}
\noindent This time, unlike the flat array of particles (which may be
very long), the array of sub-trees has at most four elements at any level
(recall that we are working with only 2 dimensions).
To build a tree, we perform recursive descent over the area:
\par{\small 
\begin{code}
  -- Perform spatial decomposition and build the tree
  buildTree :: Area -> [:Particle:] -> Tree
  buildTree area [: p :]   = Node (mass p) (location p) [::]
  buildTree area particles = Node m l subtrees
    where
      (m,l)    = calcCentroid subtrees
      subtrees = [: buildTree a ps 
                 | a <- splitArea area 
                 , let ps = [:p | p <- particles, inArea a p:]
                 , lengthP ps > 0 :]


  inArea :: Area -> Particle -> Bool
  inArea ((lx,ly),(hx,hy)) (Particle { location = (x,y) })
    = lx <= x && x <= hx && ly <= y && y <= hy

  splitArea    :: Area -> [:Area:]
     -- splitArea returns the four sub-areas in a parallel array
  calcCentroid :: [:Tree:] -> (Mass, Location)
\end{code}}
%  -- Split an area into quadrants
%  splitArea ((llx,lly),(rux,ruy)) = [:a1, a2, a3, a4:]
%         where
%          a1 = ((llx,  lly),  (midx, midy))
%           a2 = ((llx,  midy), (midx,  ruy))
%           a3 = ((midx, midy), (rux,  ruy))
%          a4 = ((midx, lly),  (rux,  midy))
%          (midx,  midy) = ((llx + rux) / 2.0 , (lly + ruy) / 2.0) 
\noindent The first equation deals with the case of a single particle: we
simply record its mass and location.  In the recursive case, the
array comprehension for @subtrees@ iterates in parallel over @(splitArea area)@,
an array of exactly four elements.  For each such area @a@, we compute
the set of particles @ps@ that lie inside @a@ and, if that set is non-empty,
we recursively call @buildTree@.  The ``if non-empty'' test discards
sub-areas which do not contain any particles at all,
so the length of @subtrees@ can be anything between 1 and 4.
We omit the implementations of @inArea@ and @calcCentroid@, since they are
straightforward.

The nested comprehension in the where clause of @buildTree@ makes sure that
@inArea@ is called on every subarea/particle combination in a single parallel
step.  Another source of nested parallelism in @buildTree@ are the recursive
calls to the parallel function @buildTree@, which are performed simultaneously
on however many sub-areas contain particles (from one to four).  The number of
parallel steps is hence proportional to the depth of the rose tree.

Lastly, we have to write the function @calcForces@,
which, given a @Tree@ and an array of particles, calculates
the forces applied by the @Tree@ on those particles.
It can do so by dividing the particles into two groups: those that
are ``far'' from the centre of gravity of the @Tree@ (as determined by a
function @isFar@), and those that are ``near''.
Here is the code:
\par{\small
\begin{code}
  calcForces :: Float -> Tree -> [:Particle:] -> [:Force:]
  calcForces len (Node m l ts) ps
    = let 
        far_forces     = [: forceOn p m l | p <- ps, isFar len l p :]
        near_ps        = [: p | p <- ps, not (isFar len l p) :]
        near_forces_s  = [: calcForces (len / 2) t near_ps | t <- ts :]
        near_forces    = [: sumForces p_forces 
                         | p_forces <- transposeP near_forces_s :]
      in
      combineP [:isFar len l p | p <- ps:] far_forces near_forces

  forceOn    :: Particle -> Mass -> Location -> Force
  isFar      :: Float -> Location -> Particle -> Bool
  sumForces  :: [:Force:] -> Force
\end{code}}
The function @calcForces@ divides the particles into two groups: those that
are ``near'' the centroid @l@ of tree, and those that are not.  For the far
particles, we simply use @forceOn@ to compute the force on each such particle
from the tree, giving @far_forces@.  For particles near to @l@, @near_ps@, we
recursively use @calcForces@ (in parallel) to compute the force on each
particle from each sub-tree giving @near_forces_s@, a short vector with one
element for each sub-tree.  Each element is a vector with one element for each
particle, giving the force on that particle from the sub-tree. All that
remains is to transpose this nested structure, and add up the forces on each
particle.  Finally, we must re-combine the near and far forces, using
@combineP@, which interleaves two vectors as directed by a boolean mask.

\subsection{Communication and locality}

Here is an alternative, simpler way to write @calcForces@:
\begin{code}
  calcForces :: Tree -> [:Particle:] -> [:Force:]
  calcForces tree ps = mapP (calc t) ps
    where
      calc (Node m l ts) p
        | isFar l p = forceOn p m l
        | otherwise = sumForces [: calc t p | t <- ts :]
\end{code}
For each particle (the @mapP@), it recurses down the tree, stopping when the
centroid of the sub-tree is far away from the particle.  

Which version should we prefer?  Different 
ways of writing the code give rise to different patterns of data communication.
In this latter version you can see that every particle needs a copy of (at least the
top part of) the tree, so the danger here is that most of the tree ends up being
copied to most of the processors.  In the earlier version, the particles migrate
(in smaller and smaller groups) to the tree, rather than the other way around.

It undoubtedly complicates the programmer's life to have to think about these
matters, but there is no silver bullet.  Parallel programming is complicated,
and programmers \emph{must} think about concurrency and communication, as well
as correctness.  However, one of the advantages of the data-parallel style is
that it gives us a much better handle on the program's \emph{cost model} (both
computation and communication) than un-structured parallel
programming~\cite{blelloch:ppa}.

\subsection{Summary}

The algorithm we have described makes extensive use of data parallelism.
For example, @buildTree@ is called in parallel on the four sub-areas of the
area under consideration; and for each of those sub-areas we compute the relevant
subset of the particles in parallel. Similarly @calcForces@ is called in parallel
on the four sub-trees; and the computation of @far_forces@ is done in parallel over
all the particles.  In each case, the recursive calls over the sub-trees
express \emph{nested} data parallelism, because the computation that is performed on
the sub-tree is itself a data-parallel computation.  This is really quite difficult
to express using flat data parallel frameworks; indeed tree construction is
often not parallelised.

The rest of this paper uses the parallel Barnes-Hut algorithm as a
running example to explain the successive steps through which the
program is compiled to run efficiently on parallel shared-memory
machines.  

% ==========================================
\section{Compiling DPH programs}

The compiler must translate high-level nested data parallel programs, as described in
the previous section, into efficient low-level code. This translation consists
of four main steps:
\begin{itemize}
\item \emph{Desugaring} removes syntactic sugar, reducing the program to a simple
lambda language.  This intermediate language, GHC's ``Core'' language, is
still strongly typed.
\item \emph{Vectorisation} transforms \emph{nested} data parallelism 
into \emph{flat} data parallelism; it is a Core-to-Core transformation.
\item \emph{Fusion} optimises the Core program, by eliminating 
redundant synchronisation points and intermediate arrays,
thus dramatically improves locality of reference; 
\item \emph{Gang parallelism} divides the parallel operations spatially into
chunks, each chunk being executed by a thread from a \emph{gang} of threads.
Typically a gang contains a thread for each CPU.  Gang parallelism is expressed
by giving library implementations of the ``vector instructions'', rather than
by built-in compiler support.
\end{itemize}
GHC implements these steps using a large number of
Core-to-Core program transformations.  Many of these transformations have
been part of GHC's optimiser for a long time, in particular a sophisticated
inliner, worker-wrapper unboxing, and constructor
specialisation~\cite{pj:trenches,ghc-inliner,pj:unboxed,pj-et-al:playing-by-the-rules}.
In the course of the Data Parallel Haskell project, we are adding more,
array-specific transformations.  Due to GHC's generic support for program
transformations --- specifically, the inliner and rewrite
rules~\cite{ghc-inliner,pj-et-al:playing-by-the-rules} --- we can implement most of
these new transformations as library code, as opposed to extending the
compiler itself.  Indeed, apart from the vectorisation pass, the rest of the
optimisation pipeline operates in ignorance of the fact that the program being
optimised is a data parallel one.

In this paper we focus mainly on vectorisation, starting at 
Section~\ref{sec:vect-overview}, after taking a brief diversion to describe
how array comprehensions are desugared (Section~\ref{sec:desugar}).

% \begin{figure*}[tp]
%   \centering
%   \includegraphics[width=\textwidth]{TheBigPicture.pdf}
%   \caption{Structure of program transformations}
%   \label{fig:big-picture}
% \end{figure*}
% 
% \simon{I'm not sure this picture really helps.  It's very small, and 
% contains much that we do not have time to discuss.}
% Figure~\ref{fig:big-picture} illustrates the anatomy of the array-specific
% transformations.  The three major components are 
% All transformations apply
% to both data types (labelled ``-- Data --'') and operations on these data
% types (labelled ``-- Control --'').  The figure names these data types and
% operations and lists the sections in which they are introduced and explained.
% 
% \simon{This paragraph seems out of place -- too much detail}
% The basic idea is that, after flattening has eliminated all nesting, we use a
% notion of \emph{distributed types} to explicitly distinguish by type the
% sequential and the parallel components of a data parallel program.  Both
% components are subjected to their own fusion mechanism; i.e., \emph{stream
%   fusion} and \emph{communication fusion}, respectively.  Finally, we
% distribute the fused parallel components across a gang of threads, each of
% which operates on a chunk of each array.


\subsection{Desugaring array comprehensions}\label{sec:desugar}

\newcommand{\dsqual}[1]{{\cal Q}\db{#1}}
\begin{figure}
$$
\begin{array}{c}
\begin{array}{lrcl}
{\bf Expressions} & e & ::= & \ldots ~|~ @[:@ e ~@|@~ q @:]@ \\
{\bf Qualifiers} & p,q & ::= & x @<-@ e ~|~ e ~|~ p@,@q ~|~ p@|@q \\
\end{array}\\ \\
\begin{array}{rcl}
{\cal D}\db{@[:@e\,@|@\,q @:]@} & = & @mapP@~(\lambda q_v.e)~\dsqual{q}
\\ \\
\multicolumn{3}{l}{\mbox{$\dsqual{q}$ computes the parallel array of the tuples generated by $q$}} \\
\dsqual{x @<-@ e} & = & e \\
\dsqual{e} & = & @if@~e~@then@~@[:():]@~@else@~@[::]@ \\
\dsqual{p@,@q} & = & @concatMapP@~ (\lambda p_v. @mapP@~ (\lambda q_v.(p_v,q_v))~\dsqual{q}) \dsqual{p} \\
\dsqual{p @|@ q} & = & @zipP@~\dsqual{p}~\dsqual{q} \\
\\
\multicolumn{3}{l}{\mbox{$q_v$ is a tuple of the variables bound by $q$}} \\
(x @<-@ e)_v & = & x \\
(g)_v & = & () \\
(p@,@q)_v & = & (p_v,q_v) \\
(p @|@ q)_v & = & (p_v,q_v) \\
\end{array}
\end{array}$$
\caption{Desugaring rules for array comprehensions}
\label{fig:desugar}
\end{figure}

In the Barnes-Hut code we used both \emph{array comprehensions}, and
\emph{ordinary functions} over parallel arrays such as @zipP@ and
@mapP@.  However, just as in the case of list comprehensions, the
former is just a convenient syntactic sugar for the latter.  More
precisely, Figure~\ref{fig:desugar} gives rules for desugaring
array comprehensions.  They are quite standard
\cite{pj-wadler:comprehensive}, and practically identical to those for
lists, so we do not discuss them further. These rules are simple, but they
should be thought of as a specification rather than an implementation, because 
they generate somewhat inefficient code.  In GHC's actual implementation 
we use slightly more complicated rules.
% \manuel{GHC desugars array
%  comprehensions rather differently. (The combination of @[:():]@/@[::]@ as
%  flags with @concatMap@ is very costly in parallel.)  For the paper, I guess,
%  we are only interested in conveying the semantics of array comprehensions
%  and for that the translation is fine, but I wonder whether we should say
%  that the standard desugaring rules for lists are too inefficient for
%  parallel structures and hence GHC actually uses a different and more
%  complicated translation.}


% =================================================
\subsection{Informal overview of vectorisation}\label{sec:vect-overview}

The purpose of vectorisation is to take a program that uses nested data 
parallelism, and transform it into a program that uses only flat
data parallelism.  Consider this tiny example
\begin{code}
  f :: Float -> Float
  f x = x*x + 1
\end{code}
For every such function we build its \emph{lifted} version \ccode{\LV{f}} thus:
\begin{alltt}
  \LV{f} :: [:Float:] -> [:Float:]
  \LV{f} x = (x \LV{*} x) \LV{+} (replicateP n 1)
         where
           n = lengthP x
\end{alltt}
\noindent Internally, \ccode{\LV{f}} uses ``vector instructions''
like \ccode{\LV{+}} to do its work, where
\begin{alltt}
   \LV{+} :: [:Float:] -> [:Float:] -> [:Float:]
\end{alltt}
\noindent Notice that it must also replicate the constant 1 so that argument has the
type that \ccode{\LV{+}} expects.  So, roughly speaking (we give the true story later), 
to form the definition of \ccode{\LV{f}} we transform the body of @f@ in the following way:
\begin{itemize}
\item Replace a constant by a call to @replicateP@.
\item Replace a function by its lifted versions (e.g. @+@ becomes \ccode{\LV{+}}).
\item Replace a parameter (e.g. @x@) by itself.
\end{itemize}
This new definition obeys the equation \ccode{\LV{f}} = @mapP@~@f@, so it
takes an array to an array.  In effect, it is a specialised variant of 
@mapP@ -- specialised by fixing the function argument.  The idea is
that whenever we see the call (@mapP f@) we will replace it by \ccode{\LV{f}}.
But there is a problem!  Suppose we have
\begin{code}
  g :: [:Float:] -> [:Float:]
  g xs = mapP f xs
\end{code}
First we replace (@mapP f@) by \ccode{\LV{f}} to get:
\begin{alltt}
  g :: [:Float:] -> [:Float:]
  g xs = \LV{f} xs
\end{alltt}
But now we must lift @g@ too, in case there are calls to (@mapP g@).  
If we try, we get this:
\begin{alltt}
  \LV{g} :: [:[:Float:]:] -> [:[:Float:]:]
  \LV{g} xs = \LV{\LV{f}} xs
\end{alltt}
Not good: we need the \emph{doubly-lifted} version of @f@! If the depth of
nesting is not statically bounded (and it isn't in Barnes-Hut) then we are in 
trouble.  Blelloch's clever solution is to observe that we can define \ccode{\LV{\LV{f}}}
in terms of \ccode{\LV{f}}, thus:
\begin{alltt}
  \LV{\LV{f}} :: [:[:Float:]:] -> [:[:Float:]:]
  \LV{\LV{f}} xss = unconcatP xss (\LV{f} (concatP xss))
\end{alltt}
That is: first concatenate all the rows of @xss@ to make a single flat vector;
then map @f@ over that vector; then chop up the result to form a vector of
vectors again, guided by the original shape of @xss@.  (Note that the 
incoming vector might well be ``ragged'', so that not all the sub-vectors have
the same length.)
At first, this idea looks terribly inefficient, because of all the flattening
and un-flattening but, as we shall see, if we choose the right data representation,
@concatP@ and @unconcatP@ take constant time and involve no copying.  
% This last sentence is referenced from Sect 4.2 (Nested arrays)

This is the core of the vectorisation transformation. We have left many details 
vague.  What about higher order functions? What about user-defined data types?
We now start to tighten our description up.
We begin by discussing how to represent arrays (Section~\ref{sec:non-parametric})
and functions (Section~\ref{sec:functions}) in vectorised code.  These representation
choices in turn drive the vectorisation transformation (Section~\ref{sec:vectorisation}).
More details are given in previous
work~\cite{keller-chak:tflat,chak-keller:more-types,leshchinskiy-etal:ho-flattening,chak-etal:DPH}.


% =============================================
\section{Representing arrays in vectorised code}
\label{sec:non-parametric}

Standard arrays in Haskell are parametric; i.e., the array representation is
independent of the type of array elements.  This is achieved by using arrays
of pointers referring to the actual element data.  Such a \emph{boxed}
representation is very flexible, but it is also detrimental to performance.
The indirections consume additional memory, increase memory traffic, and
decrease locality of memory access.  The resulting runtime penalty can be
very significant.

The parallel arrays @[:a:]@ offered by the DPH \emph{source} language
are also parametric, as can be seen from the polymorphic type signatures
in Figure~\ref{fig:parr-ops}.  One of the tasks of the vectoriser is to change
the array representation, by systematically transforming a function that manipulates values of
type @[:(Int,Int):]@, say, to one that manipulates values of type \ccode{\pa{(Int,Int)}}.
These new \ccode{\pacon} arrays have a \emph{non-parametric} representation; that is,
\emph{the representation depends on the element type}~\cite{chak-keller:more-types}.  
For example, a value of type
\ccode{\pa{Int}} is held as a contiguous memory area containing unboxed 32-bit
integer values --- not as a block of pointers to @Int@-valued thunks, as is
the case in vanilla Haskell.  

Although \ccode{\pacon} is not visible to the user, 
such non-parametric data types are an in\-de\-pen\-ently-useful source-language feature,
already implemented in GHC, which we call an \emph{associated data type}~\cite{chakravarty:assoc-data}.
We will therefore explain \ccode{\pacon} using the notation of 
associated data types.

Since the representation of an array depends on the type of its elements,
there can be no useful polymorphic functions over \ccode{PA}.  For example, we cannot
define
\begin{alltt}
  lengthPA :: \pa{a} -> Int    -- WRONG!
\end{alltt}
because @lengthPA@, being polymorphic in @a@, knows nothing about the representation of \ccode{\pa{a}}.
This is just what type classes are for.  So we declare the type \ccode{\pa} in association
with a class \ccode{\pacls} that defines operations over the type, thus:
%
\begin{alltt}
  class \pacls a where
    data \pa{a}
    indexPA     :: \pa{a} -> Int -> a
    lengthPA    :: \pa{a} -> Int
    replicatePA :: Int -> a -> \pa{a}
    ...more operations...
\end{alltt}
% 
Given a type @a@ that is allowed to be an element of a parallel array,
there is a corresponding data type \ccode{\pa{a}}, and operations @indexPA@, @lengthPA@,
@replicatePA@, and so on.  These operations therefore have overloaded types, thus:
\begin{alltt}
  indexPA  :: \pacls a => \pa{a} -> Int -> a
  lengthPA :: \pacls a => \pa{a} -> Int
  ...etc...
\end{alltt}
All the parametric operations of Figure~\ref{fig:parr-ops} have \ccode{\pacon} variants
with the same types apart from the additional \ccode{(\pacls{}~a)} constraint.
(Our concrete implementation is more complex with more operations, but the code
shown here conveys the basic idea.)  

An instance declaration fills in an implementation for each of these 
elements.  For example, the instance declaration for integers takes the
following form:
%
\begin{alltt}
  class \pacls Int where
    data \pa{Int} = AInt ByteArray
    indexPA (AInt ba) i = indexIntArray ba i
    lengthPA (AInt ba) = lengthIntArray ba
    replicatePA n i = AInt (replicateIntArray n i)
    ...more operations...
\end{alltt}
% 
We represent the array by a contiguous region of bytes (aka @ByteArray@) with
primitives such as @indexIntArray@ that operate on individual 32-bit integers
from a @ByteArray@.  (The code again simplifies the concrete
implementation by omitting the use of unboxed types.)

\subsection{Arrays of structured data} \label{sec:arr-alg}

The \ccode{\pacls} instance for @Float@, and other primitive types, follows
the same pattern.  But what about more complex data structures, such
as an array of pairs? It is quite unacceptable to represent it by an
array of pointers to (heap-allocated) records, because the indirection
costs would be too heavy. Instead, we represent it by a \emph{pair of arrays}:
%
\begin{alltt}
  class (\pacls a, \pacls b) => \pacls (a, b) where
    data \pa{(a,b)} = ATup\sub{2} Int \paB{a} \paB{b}
    indexPA (ATup\sub{2} _ arr1 arr2) i = (indexPA arr1 i, indexPA arr2 iv)
    lengthPA (ATup\sub{2} n _ _) = n
\end{alltt}
%
Thus, a \ccode{\pa{(Float,Float)}} is represented by a pair of unboxed arrays,
each storing a vector of floating point values.  Crucially, the two arrays
must have the same length; and we record that length in the
@Int@ field of the \ccode{ATup\sub{2}} constructor.  This length field is
convenient, but usually redundant --- but not always! Consider an array of
@()@ elements:
\begin{alltt}
  class \pacls () where
    data \pa{()} = ATup\sub{0} Int
    indexPA (ATup\sub{0} _) i = ()
    lengthPA (ATup\sub{0} n) = n
\end{alltt}
We need no data storage to store a vector of @()@ values, but we must still remember
its length.

Notice that the representation is
\emph{compositional}; that is, the representation of an array of pairs
is given by combining the representations of an array of the first and
second elements of the pair, and so on recursively.

The representation also allows us to combine two arrays element-wise
into an array of pairs \emph{in constant time}, with unzipping being
equally easy:
\begin{alltt}
  zipPA :: PAElem a => PA a -> PA b -> PA (a,b)
  zipPA as bs = ATup\sub{2} (lengthPA as) as bs

  unzipPA :: PA (a,b) -> (PA a, PA b)
  unzipPA (ATup\sub_2 _ as bs) = (as,bs)
\end{alltt}
This stands in contrast to lists, where zipping and unzipping take linear time.

Lastly, since records are converted into product types by the
desugarer, the @Particle@ arrays in Barnes-Hut are
represented by tuples of arrays.  

\subsection{Nested arrays}

Even more interesting is the representation of nested arrays. A classic example is
that of \emph{sparse matrices}, in which we represent a sparse matrix as a vector
of rows, each row consisting of a vector of (index,value) pairs, where only the
non-zero values in the row are represented.  Thus
\begin{code}
  type SparseMatrix a = [:[:(Int,a):]:]
\end{code}
Since our ultimate goal is to eliminate nested parallelism, it is not
surprising that we also want to represent nested arrays in terms of
flat ones. Indeed, a nested array \ccode{\pa{\paB{a}}} can be encoded by
\begin{itemize}
\item a flat \emph{data array} of type \ccode{\pa{a}} which contains the data elements and
\item a \emph{segment descriptor} of type \ccode{\pa{(Int, Int)}} which stores the
  starting position and length of the subarrays embedded in the flat
  data array.
\end{itemize}
%
This is captured by the following instance:
\begin{alltt}
  class \pacls a => \pacls \paB{a} where
    data \pa{\paB{a}} = AArr \paB{a} \paB{(Int, Int)} 
    indexPA (AArr arr segd) i = slicePA arr (indexPA segd i)
    lengthPA (AArr _ seg) = lengthPA seg
\end{alltt}
%
where @sliceP@ extracts a subarray from a larger array in constant time.
% NOTE: the constant time remark is referenced at the end of 5.2
Thus, the sparse matrix
%
\begin{code}
  [:[:(0,15),(2,9),(3,20):], [::], [:(3,46):]:]
\end{code}
%
will be represented as
%\pagebreak[4] % !!!Urgh!
%
\begin{alltt}
  AArr (ATup\sub{2} [#0,2,3,3#] [#15,9,20,46#]) -- Data
       (ATup\sub{2} [#0,3,3#] [#3,0,1#])        -- Segment descriptor
\end{alltt}
%
where we write @[#...#]@ for a literal @ByteArray@. 
The first @ByteArray@ contains all the column indexes, the second one all the
@Float@s, and the third and fourth the start indexes and lengths of the
segments, respectively.  Since all four @ByteArray@ are unboxed, programs which
process such matrices can be compiled to highly efficient code.

Remarkably, we can now give constant-time implementations of the two functions
@concatPA@ and @unconcatPA@, as promised in Section~\ref{sec:vect-overview}:
\begin{alltt}
  concatPA :: \pa{\paB{a}} -> \pa{a}
  concatPA (AArr cts _) = cts
  
  unconcatPA :: \pa{\paB{a}} -> \pa{b} -> \pa{\paB{b}}
  unconcatPA (AArr _ shape) cts = AArr cts shape
\end{alltt}

\begin{figure}[tbp]
  \centering
  \begin{pgfpicture}{0cm}{0cm}{11cm}{3cm}
%     \mpgfline{0cm}{1cm}{12cm}{1cm}
%     \mpgfline{0cm}{2cm}{12cm}{2cm}
%     \mpgfrect{0cm}{0cm}{12cm}{3cm}

%     \mpgfline{1cm}{0cm}{1cm}{3cm}
%     \mpgfline{2cm}{0cm}{2cm}{3cm}
%     \mpgfline{3cm}{0cm}{3cm}{3cm}
%     \mpgfline{4cm}{0cm}{4cm}{3cm}
%     \mpgfline{5cm}{0cm}{5cm}{3cm}

    % boxes 

    \mpgfbox{1cm}{2.5cm}
    \mpgfbox{1.6cm}{2.5cm}
    \mpgfbox{2.3cm}{2.5cm}
    \mpgfbox{3.1cm}{2.5cm}

    % arrows 
    \pgfsetendarrow{\pgfarrowsingle}

    \mpgfline{1cm}{2.5cm}{0.9cm}{2cm}
    \mpgfline{1.2cm}{2.5cm}{1.3cm}{2cm}

    \mpgfline{1.7cm}{2.5cm}{1.7cm}{2cm}

    \mpgfline{2.3cm}{2.5cm}{2.1cm}{2cm}
    \mpgfline{2.4cm}{2.5cm}{2.4cm}{2cm}
    \mpgfline{2.5cm}{2.5cm}{2.7cm}{2cm}
    
    \mpgfline{3.1cm}{2.5cm}{3.0cm}{2cm}
    \mpgfline{3.3cm}{2.5cm}{3.4cm}{2cm}

    % boxes

    \mpgfbox{0.8cm}{1.8cm}
    \mpgfbox{1.2cm}{1.8cm}
    \mpgfbox{1.6cm}{1.8cm}
    \mpgfbox{2.0cm}{1.8cm}
    \mpgfbox{2.3cm}{1.8cm}
    \mpgfbox{2.6cm}{1.8cm}
    \mpgfbox{2.9cm}{1.8cm}
    \mpgfbox{3.3cm}{1.8cm}

    % arrows

    \mpgfline{0.9cm}{1.8cm}{0.9cm}{1.3cm}    
    \mpgfline{1.6cm}{1.8cm}{1.5cm}{1.3cm}    
    \mpgfline{1.8cm}{1.8cm}{1.9cm}{1.3cm}    
    \mpgfline{2.4cm}{1.8cm}{2.4cm}{1.3cm}    
    \mpgfline{2.7cm}{1.8cm}{2.7cm}{1.3cm}    
    \mpgfline{3.4cm}{1.8cm}{3.4cm}{1.3cm}    

    % boxes
    \mpgfbox{0.8cm}{1.1cm}
    \mpgfbox{1.4cm}{1.1cm}
    \mpgfbox{1.8cm}{1.1cm}
    \mpgfbox{2.3cm}{1.1cm}
    \mpgfbox{2.6cm}{1.1cm}
    \mpgfbox{3.3cm}{1.1cm}


    \pgfputat{\pgfxy(0.7,2.5)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(1.3,2.5)}{\pgfbox[center,base]{\texttt{,}}}
    \pgfputat{\pgfxy(2.0,2.5)}{\pgfbox[center,base]{\texttt{,}}}
    \pgfputat{\pgfxy(2.8,2.5)}{\pgfbox[center,base]{\texttt{,}}}
    \pgfputat{\pgfxy(3.6,2.5)}{\pgfbox[center,base]{\texttt{]}}}

    % end tree 



    % flat structure
    % l1
    \mpgfbox{6.5cm}{2.5cm}
    \mpgfbox{7.4cm}{2.5cm}
    \mpgfbox{8.3cm}{2.5cm}
    \mpgfbox{9.2cm}{2.5cm}

    \pgfputat{\pgfxy(6.3,2.5)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(7.1,2.5)}{\pgfbox[center,base]{\texttt{,}}}
    \pgfputat{\pgfxy(8.0,2.5)}{\pgfbox[center,base]{\texttt{,}}}
    \pgfputat{\pgfxy(8.9,2.5)}{\pgfbox[center,base]{\texttt{,}}}
    \pgfputat{\pgfxy(9.7,2.5)}{\pgfbox[center,base]{\texttt{]}}}

    \mpgfline{7.9cm}{2.37cm}{7.9cm}{2.1cm}    

    %l2
    \pgfputat{\pgfxy(6.15,1.8)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(6.25,1.8)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(6.95,1.8)}{\pgfbox[center,base]{\texttt{]}}}
    \pgfputat{\pgfxy(7.1,1.8)}{\pgfbox[center,base]{\texttt{,}}}

    \mpgfbox{6.35cm}{1.8cm}
    \mpgfbox{6.65cm}{1.8cm}
        
    \pgfputat{\pgfxy(7.3,1.8)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(7.7,1.8)}{\pgfbox[center,base]{\texttt{]}}}
    \pgfputat{\pgfxy(7.8,1.8)}{\pgfbox[center,base]{\texttt{,}}}

    \mpgfbox{7.4cm}{1.8cm}

    \pgfputat{\pgfxy(7.9,1.8)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(8.9,1.8)}{\pgfbox[center,base]{\texttt{]}}}
    \pgfputat{\pgfxy(8.975,1.8)}{\pgfbox[center,base]{\texttt{,}}}
    \mpgfbox{8cm}{1.8cm}
    \mpgfbox{8.3cm}{1.8cm}
    \mpgfbox{8.6cm}{1.8cm}

    \pgfputat{\pgfxy(9.05,1.8)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(9.75,1.8)}{\pgfbox[center,base]{\texttt{]}}}
    \pgfputat{\pgfxy(9.85,1.8)}{\pgfbox[center,base]{\texttt{]}}}
    \mpgfbox{9.15cm}{1.8cm}
    \mpgfbox{9.45cm}{1.8cm}

    \mpgfline{7.9cm}{1.67cm}{7.9cm}{1.4cm}    

    %l3
    \pgfputat{\pgfxy(-0.4,0)}{
    \pgfputat{\pgfxy(6.15,1.1)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(6.25,1.1)}{\pgfbox[center,base]{\texttt{[}}}
    
    \mpgfbox{6.35cm}{1.1cm}
     \pgfputat{\pgfxy(6.65,1.1)}{\pgfbox[center,base]{\texttt{]}}}
     \pgfputat{\pgfxy(6.75,1.1)}{\pgfbox[center,base] {\texttt{,}}}
     \pgfputat{\pgfxy(6.85,1.1)}{\pgfbox[center,base]{\texttt{[}}}
     \pgfputat{\pgfxy(7.05,1.1)}{\pgfbox[center,base]{\texttt{]}}}
     \pgfputat{\pgfxy(7.15,1.1)}{\pgfbox[center,base] {\texttt{,}}}
     \pgfputat{\pgfxy(7.25,1.1)}{\pgfbox[center,base] {\texttt{[}}}
    \mpgfbox{7.35cm}{1.1cm}    
     \pgfputat{\pgfxy(7.65,1.1)}{\pgfbox[center,base] {\texttt{,}}}
    \mpgfbox{7.75cm}{1.1cm}    
     \pgfputat{\pgfxy(8.05,1.1)}{\pgfbox[center,base] {\texttt{]}}}
     \pgfputat{\pgfxy(8.15,1.1)}{\pgfbox[center,base] {\texttt{,}}}

     \pgfputat{\pgfxy(8.25,1.1)}{\pgfbox[center,base]{\texttt{[}}}
     \pgfputat{\pgfxy(8.45,1.1)}{\pgfbox[center,base]{\texttt{]}}}

     \pgfputat{\pgfxy(8.55,1.1)}{\pgfbox[center,base] {\texttt{,}}}
     \pgfputat{\pgfxy(8.65,1.1)}{\pgfbox[center,base] {\texttt{[}}}
    \mpgfbox{8.75cm}{1.1cm}    
     \pgfputat{\pgfxy(9.05,1.1)}{\pgfbox[center,base] {\texttt{]}}}


     \pgfputat{\pgfxy(9.15,1.1)}{\pgfbox[center,base] {\texttt{,}}}
     \pgfputat{\pgfxy(9.25,1.1)}{\pgfbox[center,base] {\texttt{[}}}
    \mpgfbox{9.35cm}{1.1cm}    
     \pgfputat{\pgfxy(9.65,1.1)}{\pgfbox[center,base] {\texttt{]}}}
     \pgfputat{\pgfxy(9.75,1.1)}{\pgfbox[center,base] {\texttt{,}}}

     \pgfputat{\pgfxy(9.85,1.1)}{\pgfbox[center,base]{\texttt{[}}}
     \pgfputat{\pgfxy(10.05,1.1)}{\pgfbox[center,base]{\texttt{]}}}

     \pgfputat{\pgfxy(10.15,1.1)}{\pgfbox[center,base] {\texttt{,}}}
     \pgfputat{\pgfxy(10.25,1.1)}{\pgfbox[center,base] {\texttt{[}}}
    \mpgfbox{10.35cm}{1.1cm}    
     \pgfputat{\pgfxy(10.65,1.1)}{\pgfbox[center,base] {\texttt{]}}}
     \pgfputat{\pgfxy(10.75,1.1)}{\pgfbox[center,base] {\texttt{]}}}}

    \mpgfline{7.9cm}{0.97cm}{7.9cm}{0.7cm}    

    %l4
    \pgfputat{\pgfxy(0.5,0)}{
    \pgfputat{\pgfxy(6.15,0.4)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(6.25,0.4)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(6.45,0.4)}{\pgfbox[center,base]{\texttt{]}}}
    \pgfputat{\pgfxy(6.55,0.4)}{\pgfbox[center,base]{\texttt{,}}}

    \pgfputat{\pgfxy(6.65,0.4)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(6.85,0.4)}{\pgfbox[center,base]{\texttt{]}}}
    \pgfputat{\pgfxy(6.95,0.4)}{\pgfbox[center,base]{\texttt{,}}}

    \pgfputat{\pgfxy(7.05,0.4)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(7.25,0.4)}{\pgfbox[center,base]{\texttt{]}}}
    \pgfputat{\pgfxy(7.35,0.4)}{\pgfbox[center,base]{\texttt{,}}}

    \pgfputat{\pgfxy(0.4,0)}{
    \pgfputat{\pgfxy(7.05,0.4)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(7.25,0.4)}{\pgfbox[center,base]{\texttt{]}}}
    \pgfputat{\pgfxy(7.35,0.4)}{\pgfbox[center,base]{\texttt{,}}}
}

    \pgfputat{\pgfxy(0.8,0)}{
    \pgfputat{\pgfxy(7.05,0.4)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(7.25,0.4)}{\pgfbox[center,base]{\texttt{]}}}
    \pgfputat{\pgfxy(7.35,0.4)}{\pgfbox[center,base]{\texttt{,}}}
  }
    \pgfputat{\pgfxy(1.2,0)}{
    \pgfputat{\pgfxy(7.05,0.4)}{\pgfbox[center,base]{\texttt{[}}}
    \pgfputat{\pgfxy(7.25,0.4)}{\pgfbox[center,base]{\texttt{]}}}
    \pgfputat{\pgfxy(7.35,0.4)}{\pgfbox[center,base]{\texttt{]}}}
  }}

\pgfsetlinewidth{0.8pt}
\mpgfline{4cm}{1.75cm}{5.5cm}{1.75cm}    
  \end{pgfpicture}
\caption{Value of type \ccode{[:Tree:]} and its vectorised representation}
  \label{fig:tree-vect}
\end{figure}

\subsection{Recursive types}\label{sec:rec-types}

If the array elements are recursive, the non-parametric representation of the
array has to be recursive, too. For instance, arrays of \ccode{Tree} from
Section~\ref{sec:bh-dph} are represented as follows:
%
\begin{alltt}
  instance PAElem Tree where
    data \pa{Tree} = ATree Int \paB{Mass} \paB{Location} \paB{\paB{Tree}}
\end{alltt}
%
Since \ccode{Tree} is a product, the representation is similar to
arrays of tuples. It stores the masses and locations of the centroids and a
nested array containing the subtrees of each node. As
described in the previous section, the latter is encoded by a flat array
of trees together with a segment descriptor. In effect, this means that an
array of trees is represented by a list with each element containing the
centroids and segmentation information for one tree level. This allows all
data in one level to be processed in parallel, although the levels
have to be processed one after another. Figure~\ref{fig:tree-vect} illustrates
this representation. User-defined types are discussed in more detail in
Section~\ref{sec:user-types}.

\subsection{Polymorphism}

If we were only interested in \emph{monomorphic} code, or if we would use a
whole-program compiler that specialises a polymorphic to a monomorphic
program, as was the case in NESL, then life would have been much easier. 
We could implement the non-parametric array type by \emph{statically replacing}
\ccode{\pa{Int}} by \ccode{PAInt}, say, where the latter is a perfectly
ordinary data type, defined as
\begin{alltt}
  data PAInt = AInt ByteArray
\end{alltt}
Now we do not need non-parametric types; in the vectorised code, original
types @[:Int:]@, @[:(Int,Float):]@, etc, are simply replaced by @PAInt@,
@PAPair PAInt PAFloat@, and so on, where all these are ordinary data types.

Alas, this does not work for \emph{polymorphic} functions.  For example, how
could we translate the type of this function?  What would we statically
replace @[:a:]@ by?
\begin{alltt}
  firstRow :: [:[:a:]:] -> [:a:]
\end{alltt}
We also cannot turn polymorphic functions into families of monomorphic
functions, as we support separate compilation and polymorphic recursion.  No ---
if we want polymorphism, we must use something akin to type classes, as we
have described in this section.  A key component of our work is the extension
of the non-parametric representation idea to work in a polymorphic setting.
In particular, our Core language regards \ccode{\pacon} as a type-level
function from types to types \cite{sulzmann-et-al:fc}.

\section{Representing functions in vectorised code} \label{sec:functions}

Haskell is a higher order language, so we have to consider how to vectorise
programs that manipulate functions.  Vectorising higher-order programs raises
two distinct problems.

\subsection{Functions are pairs} \label{sec:fun-pairs}

Consider this (contrived) definition:
\begin{alltt}
  ho :: (Int->Bool) -> (Bool, [:Bool:])
  ho f = (f 2, mapP f [:1,2,3:])
\end{alltt}
In our overview (Section~\ref{sec:vect-overview}), we said that we
should replace @(mapP f)@ with a call to \ccode{\LV{f}}, the lifted
version of @f@.  But since @f@ is lambda-bound, it is not so easy to
call ``the lifted version of @f@''.  Clearly the
caller must pass the lifted version of @f@ as a parameter to @ho@.
But there is an ordinary,
scalar call @(f x)@ in the body of @ho@, so we can't pass \emph{only} the 
lifted version.  The
obvious alternative is to pass a pair that gives \emph{both} the
lifted \emph{and} unlifted versions. With such a representation, @mapP@
can just extract the lifted variant of its argument (a pair), 
while the vanilla application of @f@ extracts the unlifted variant.

\subsection{Functions are closures}

There is a second challenge.  In Section~\ref{sec:non-parametric} 
we discussed the efficient, non-parametric representation of data-parallel 
arrays.  A higher order language forces us to confront the question of
how to represent an array of \emph{functions}.  For example:
\begin{alltt}
  distance :: [:Float:] -> [: Float->Float :]
  distance xs = mapP (\hslam x y. sqrt (x*x + y*y)) xs

  distY :: [:Float:] -> Float -> Float
  distY xs y = sumP [: d y | d <- distance xs :]
\end{alltt}
There are more direct ways to write @distY@, of course, but arrays
of functions also arise inevitably when we think about lifting. If we
start with
\begin{alltt}
  f :: (Int->Int) -> Int
\end{alltt}
then the lifted version of @f@ has the type
\begin{alltt}
  \LV{f} :: [: Int->Int :] -> [:Int:]
\end{alltt}
How should we represent such an array of functions?
A possible answer is ``as an array of pointers to function closures''.  
Bad answer!  In data-parallel computations, all the processors are supposed
to execute the same code in parallel, but if every function closure in the
array is potentially different, they clearly cannot do that.  Furthermore,
a pointer-based representation destroys locality.  

Happily, there is no need for the generality of a distinct function
pointer for each array element.  Consider the result of @distance@, for example.
Every element of this array is a function \emph{with the same code}, but
a different value for the function's free variable @x@.  This makes the
solution obvious: we must represent an array of functions by a pair of
a \emph{single} code pointer and an \emph{array} of environment records, 
which give the per-element bindings for the free variables.

\subsection{Putting it together}

Putting our two solutions together, we see that in vectorised code a function
must be represented by a triple:
\begin{enumerate}
\item The scalar version of the function
\item The lifted version of the function
\item An environment record of the free variables of the function
\end{enumerate}
To be concrete, here is the data type declaration for vectorised functions:
\begin{alltt}
  data (a :-> b) = forall e. \pacls e => 
                   Clo \{ env  :: e
                       , clo\sub{s} :: e -> a -> b
                       , clo\sub{l} :: \pa{e} -> \pa{a} -> \pa{b} \}
\end{alltt}
This declaration says that @(:->)@ is an algebraic data type (written infix),
with a single constructor @Clo@.  The constructor @Clo@ has
an existentially-quantified type variable
@e@, and three fields, @env@, \ccode{clo\sub{s}}, and \ccode{clo\sub{l}}.  
The vectorisation transformation will transform every function type
$\tau_1 @->@ \tau_2$ to $\tau_1' @:->@ \tau_2'$, where $\tau_1'$ is the transformed version of $\tau_1$
and similarly for $\tau_2$.  In effect, the vectorisation transform performs
closure conversion~\cite{appel-jim:CPS-CPS}.

With this definition in hand, we can now explain how arrays of 
values of type @(a :-> b)@ are represented: 
\begin{alltt}
  instance \pacls (a :-> b) where
    data \pa{(a :-> b)} = forall e. \pacls e 
                       => AClo \{ aenv  :: \pa{e}
                               , aclo\sub{s} :: e -> a -> b
                               , aclo\sub{l} :: \pa{e} -> \pa{a} -> \pa{b} \}
    lengthPA (AClo env f\sub{s} f\sub{l}) = lengthPA env
    indexPA (AClo env f\sub{s} f\sub{l}) n = Clo (indexPA env n) f\sub{s} f\sub{l}
    replicatePA n (Clo env f\sub{s} f\sub{l}) = AClo (replicatePA n env) f\sub{s} f\sub{l}
\end{alltt}
To represent an array of functions, we keep a \emph{single} code pointer for
each of the 
scalar and lifted code, but have an \emph{array} of environment records.
Notice that @Clo@ and @AClo@ differ \emph{only} in the type of the environment field; their
\ccode{clo\sub{s}} and \ccode{clo\sub{l}} fields are identical.

As in the case of other types, it is worth noting that this representation supports very simple and
direct implementations of indexing, replication, and so on.
It does \emph{not} efficiently support literal arrays of various different
functions, such as @[:sin,cos:]@. This is quite deliberate: in a data-parallel
computation all the processors should be performing the same computation at the
same time.  Nevertheless, such arrays can be handled, essentially using
conditionals which ensure that different functions are executed one after
another.

In Section~\ref{sec:ho-vect}, we will see how @Clo@ and @AClo@ are used
in the transformation of the @ho@ example.  Before we can do so, we must
first specify the vectorisation transformation more precisely.

% However, we do not want to
% represent it exactly in the same way as, for example, the array of functions
% \begin{alltt}
%   [: (+1) ,(*2), (-3):]
% \end{alltt}
% which are a collection of distinct functions, which have to be applied
% sequentially in a data parallel execution model. The first type of
% values appears quite frequently: namely, whenever a partial application is
% lifted. The second only whenever such an array is constructed
% explicitely.  We therefore need gear our representation of 
% functions towards the first case (while still being able to handle the
% second case, of course).
% With this in mind, we no longer need to say
% ``replace @(mapP f)@ by \ccode{\LV{f}}''.  Instead we can define
% @\VV{mapP}@ to select the second component of its function argument:
% \begin{alltt}
%   \VV{mapP} f as = AClo \SE{f} \LV{f} (lengthP as)  \LV{\capp} as
% \end{alltt}
% which reduces to @\LV{f} (lengthP as) as@.
% 
% The partial application in our previous example, @mapP (+) [:1,2,3:]@,
% is now represented by @AClo (+) (\LV{+}) [:1,2,3:]@. The closure
% environment in this expression also fixes the length of arrays
% acceptable as second argument. 



% =================================================================
\section{Vectorisation}
\label{sec:vectorisation}

\begin{figure}
$$\begin{array}{c}
\mbox{$\TT{\tau}:: \mathit{Type} \to \mathit{Type}$ 
is the vectorisation transformation on types} \\[2mm]
\begin{array}{rcll}
\TT{\tau_1 @->@ \tau_2} & = & \TT{\tau_1}~ @:->@~\TT{\tau_2}  & \mbox{Functions} \\
\TT{@[:@\tau@:]@}    & = & \TTL{\tau} & \mbox{Parallel arrays} \\
\TT{@Int@} & = & @Int@  & \mbox{Primitive scalar types}\\
\TT{@Float@} & = & @Float@ \\
\TT{T~\tau_1 \ldots \tau_n} & = & \VV{T}~\TT{\tau_1}\ldots \TT{\tau_n} & \mbox{Algebraic data types (e.g. lists)}
\\ \\
\TTL{\tau} & = & \pa{\TT{\tau}} 
\end{array} \\ \\[-2mm]
\hline \\[-2mm]
\mbox{$\VT{e}::\mathit{Term} \to \mathit{Term}$ 
is the full vectorisation transformation on terms} \\
\mbox{Invariant: if $\tup{x_i : \sigma_i} \vdash e:\tau$ then 
      $\tup{x_i:\TT{\sigma_i}} \vdash \VT{e}:\TT{\tau}$}\\[1em]
\begin{array}{rll}
\VT{k}       & = k & \mbox{$k$ is a literal} \\       
\VT{f}       & = \VV{f} & \mbox{$f$ is bound at top level} \\
\VT{x}       & = x     & \mbox{$x$ is locally bound (lambda, let, etc)} \\
\VT{e_1 \; e_2} & = \VT{e_1} ~\capp~ \VT{e_2} \\
\VT{\lambda{x}.e} 
  & \multicolumn{2}{l}{= @Clo@\begin{array}[t]{l} <
              \{ @env@ = (y_1,\dots,y_k) \\
              , @clo@_s = \lambda e\,x.\,
                 @case@\,e\,@of@\;(y_1,\dots,y_k) \to \VT{e} \\
              , @clo@_l = \lambda e\,x.\,
                 @case@\, e\, @of@\; @ATup@_k\; n\; y_1 \dots y_k \to \LT{n}{e} \}
           \end{array}} \\
 & \multicolumn{2}{l}{\mbox{where} \{y_1, \ldots, y_k\} = \mbox{free variables of $\lambda x.e$}} \\
\VT{\begin{array}{l}
        @if@\,e_1 \\
      @then@\,e_2 \\
      @else@\, e_3
  \end{array}} & 
\multicolumn{2}{l}{= @if@\,\VT{e_1}\,@then@\,\VT{e_2}\,@else@\,\VT{e_3}} \\
\end{array} \\ \\[-2mm]
 \hline \\[-2mm]
\mbox{$ \LT{n}{e}:: \mathit{Term} \to \mathit{Term} \to \mathit{Term}$
is the lifting transformation on terms} \\
\mbox{Invariant:} \begin{array}[t]{l}
  \mbox{if $\tup{x_i : \sigma_i} \vdash e : \tau$ 
        then $\tup{x_i : \TTL{\sigma_i}} \vdash \LT{n}{e}: \TTL{\tau}$} \\
  \mbox{where $n$ is the length of the result array}
  \end{array} \\[2em]
\begin{array}{rll}
\LT{n}{k} & = @replicatePA@~ n~ k & \mbox{$k$ is a literal} \\
\LT{n}{f} & = @replicatePA@~ n~ \VV{f} & \mbox{$f$ is bound at top level} \\
\LT{n}{x} & = x  & \mbox{$x$ is locally bound (lambda, let, etc)} \\
\LT{n}{e_1\; e_2} & = \LT{n}{e_1}~~ \LV{\capp}~~ \LT{n}{e_2} \\
\LT{n}{\lambda{x}.e} & \multicolumn{2}{l}{= @AClo@ \begin{array}[t]{l}
        \{ @aenv@ = @ATup@_k\; n\; y_1 \dots y_k, \\
        , @aclo@_s = \lambda e\,x.\,
           @case@\,e\,@of@\; (y_1,\dots,y_k) \to \VT{e} \\
        , @aclo@_l = \lambda e\,x.\,
           @case@\,e\,@of@\; @ATup@_k\; n'\; y_1 \dots y_k \to \LT{n'}{e} \}
     \end{array}} \\
 & \multicolumn{2}{l}{\mbox{where} \{y_1, \ldots, y_k\} = \mbox{free variables of $\lambda x.e$}} \\
\raisebox{0pt}[0.9\height][0.2\depth]{\LT{n}{\begin{array}{l}
        @if@\,e_1 \\
      @then@\,e_2 \\
      @else@\, e_3
  \end{array}}} & = @combinePA@~~e_1'~~e_2'~~e_3' \\
& \multicolumn{2}{l}{\mbox{where}\begin{array}[t]{l}
 e_1' = \LT{n}{e_1} \\
 e_2' = @case@ ~\mathit{ys_2}~@of@ \; @ATup@_k\; n_2\; y_1 \dots y_k \to \LTprime{n_2}{e_2} \\
 e_3' = @case@ ~\mathit{ys_3}~@of@ \; @ATup@_k\; n_3\; y_1 \dots y_k \to \LTprime{n_3}{e_3} \\
 (\mathit{ys_2},\mathit{ys_3}) = @splitPA@~e_1'~ (@ATup@_k\; n\; y_1 \dots y_k) \\
 \{y_1, \ldots, y_k\} = \mbox{free variables of $e_2, e_3$}
 \end{array}}\\
\LTprime{n}{e} & \multicolumn{2}{l}{= ~ @if@~n@==@0~@then@~~@emptyPA@~~@else@~~\LT{n}{e}}
\end{array} \\
\end{array}$$
\caption{The vectorisation transformation} \label{fig:vect}
\end{figure}

We are finally ready to discuss the vectorisation transformation itself.
Consider a top-level function definition $f::\tau = e$,
where $\tau$ is the type of \ccode{\var f}.  
The \emph{full vectorisation} transformation produces a definition for
the vectorised version of \var{f} called $\VV{\var{f}}$, thus:
%
$$\VV{\var{f}} :: \TT{\tau} = \VT{\var{e}}$$
%
Here, \ccode{\VV{\var f}} is the \emph{fully vectorised} variant of \ccode{\var f}, 
whose right-hand side is generated by the \emph{full vectorisation
  transform} \VT{\cdot}.  As we have already seen, vectorisation returns an expression of a
different type to the input, so the type of \ccode{\VV{\var f}} is obtained by
vectorising the type $\tau$, thus $\TT{\tau}$.  In general, if
$e::\tau$ then $\VT{\var e}::\TT{\tau}$.
Figure~\ref{fig:vect} gives the functions for both type and term vectorisation.
In our compiler, the transformation applies to an explicitly-typed program,
but we omit all type information in Figure~\ref{fig:vect}, in order to
concentrate on the essentials.

Vectorisation is applied separately to each top level function in the
program, so it is a whole-program transformation.  In real programs,
only a part will be data-parallel, while much of it is not (e.g. input/output,
user interaction etc).  We ignore this issue here, but in reality our
compiler performs \emph{selective} vectorisation -- see
Section~\ref{sec:selective} and \cite{chakravarty-et-al:partial-vect}.

The type transformation $\TT{\tau}$ transforms a source-program type
to the corresponding type in the vectorised program.  As can be seen
in Figure~\ref{fig:vect}, its effect is simple: it transforms every
function arrow @(->)@ to a vectorised function arrow @(:->)@, and
every parametric array constructor @[::]@ to a non-parametric parallel
array constructor \ccode{\pacon}. A user-defined algebraic data type might
have nested uses of @(->)@ or @[::]@ --- for example, @Tree@ does so ---
and for these we must generate a vectorised variant ($\VV{@Tree@}$)
of the data type itself.  We elaborate this point in Section~\ref{sec:user-DT}.

This type transformation forces the 
vectorised program to differ quite radically from the original.
In particular, since a ``function''
is now a triple constructed with @Clo@, we need an infix application 
operator \capp{} to extract the scalar copy:
\begin{alltt}
  (\capp) :: (a :-> b) -> a -> b
  (\capp) (Clo env fs fl) = fs env
\end{alltt}
and a lifted version
\begin{alltt}
  (\LV{\capp}) :: \pa{(a:->b)} -> \pa{a} -> \pa{b}
  (\LV{\capp}) (AClo env fs fl) = fl env
\end{alltt}

The transformation rules in Figure~\ref{fig:vect} are given with their type
invariants, which make the rules much more comprehensible.  For example,
consider the rule for $\VT{e_1\;e_2}$.  Since $e_1 : \tau_1 @->@ \tau_2$, we
know that $\VT{e_1} : \TT{\tau_1} @:->@ \TT{\tau_2}$; that is why we need the
application function $(\capp)$ to transform the application to an expression
of type $\TT{\tau_2}$.

Similarly, the rule for $\VT{\lambda x.e}$ must produce a value of
type $\TT{\tau_1} @:->@ \TT{\tau_2}$, and that in turn must be built
with a @Clo@ constructor.  We build an environment tuple $(y_1, \ldots, y_k)$, of the
free variables of $(\lambda x.e)$.  Now the type of the arguments of @Clo@ tell us
what functions we must build.  The scalar function simply requires a recursive
use of $\VT{e}$, while the lifted function requires us to generate a lifted version of
the code for $e$, $\LT{n}{e}$.  

The rules for $\LT{n}{e}$ can be read in the same way.
The main new complication is with conditionals.  First we compute in parallel $e_1'$, the vector
of booleans (of length $n$) for the discriminant of the conditional.  Then
we use that vector to split the a vector of environment tuples into two parts,
$ys_2$ (for which corresponding elements of $e_1'$ is true), and $ys_3$ (for which
$e_1'$ is false). The lengths $n_2, n_3$ of these vectors will sum to $n$.
Then we compute each of $e_2'$ and $e_3'$ in parallel, and finally interleave them
together with @combinePA@. 

Why do we need to pack and split the free variables in the conditional
rule?  Each free variable $y_i$ is bound to an $n$-vector; but in
the @then@ branch we need a (shorter) $n_2$-vector (namely $ys_2$) of
the elements of $y_i$ for which $e$ is True; and dually for the @else@
branch. We must also test for $n_2$ or $n_3$ being
zero (done by $\LTprime{n}{e}$), otherwise when transforming a recursive
function we would generate a program that recurses infinitely deep.
The \emph{operational} behaviour of the translated function will 
compute $e_1'$, $e_2'$ and $e_3'$ \emph{in sequence}; as in any data-parallel 
machine, the ``then'' and ``else'' branches of a conditional are computed
separately.

Figure~\ref{fig:vect} is the core of this paper.  Our real system
handles @let@ expressions, @case@ expressions, and constructors, and
hence is a bit more complicated.  But Figure~\ref{fig:vect} describes all the 
essential ideas.

Since $\cal V$ invokes both $\cal V$ and $\cal L$ (as does $\cal L$) you
might worry about a code explosion.  But notice that the \ccode{clo\sub{s}} field in
$\cal V$ is identical to the \ccode{aclo\sub{s}} field in $\cal L$, and both are closed
functions that can be named, and bound at top level; and 
similarly for the \ccode{clo\sub{l}} and \ccode{aclo\sub{l}} fields.
Hence, as we will see in the examples that follow, we can avoid the code explosion simply by
naming and sharing these functions.

\subsection{A simple example}

Here is the simplest possible example:
\begin{alltt}
  inc :: Float -> Float
  inc = \hslam{}x. x + 1
\end{alltt}
The full vectorisation transformation in Figure~\ref{fig:vect}
gives us this:
\begin{alltt}
  \VV{inc} :: Float :-> Float
  \VV{inc} = Clo () inc\sub{S} inc\sub{L}

  \SE{inc} :: () -> Float -> Float
  \SE{inc} = \hslam{}e x. case e of () -> \VV{(+)} \$: x \$: 1 

  \LV{inc} :: \pa{()} -> \pa{Float} -> \pa{Float}
  \LV{inc} = \hslam{}e x. case e of ATup\sub{0} n -> \VV{(+)} \LV{\$:} x \LV{\$:} (replicatePA n 1)
\end{alltt}
To aid explanation we have named \ccode{\SE{inc}} and \ccode{\LV{inc}},
but otherwise we have simply applied Figure~\ref{fig:vect} blindly.
Notice the way we have systematically transformed @inc@'s type, replacing @(->)@
by @(:->)@.  Notice too that this transformation neatly embodies the idea that
we need two versions of every top-level function @inc@, a \emph{scalar
  version} \ccode{\SE{inc}} and a \emph{lifted version} \ccode{\LV{inc}}.  These two
versions paired together to form the \emph{fully vectorised version}
\VV{@inc@}.

The vectorised code makes use of vectorised addition \ccode{(+)}, which is part
of a fixed, hand-written library of vectorised primitives:
\begin{alltt}
  \VV{(+)} :: Float :-> Float :-> Float
  \VV{(+)} = Clo () \SE{(+)} \LV{(+)}

  \SE{(+)} :: () -> Float -> Float :-> Float
  \SE{(+)} = \hslam{}e x. Clo x \SE{addFloat} \LV{addFloat}
 
  \LV{(+)} :: \pa{()} -> \pa{Float} -> \pa{(Float:->Float)}
  \LV{(+)} = \hslam{}e xs. AClo xs \SE{addFloat} \LV{addFloat}

  -- Implemented in the back end
  addFloat\sub{S} :: Float -> Float -> Float
  addFloat\sub{L} :: \pa{Float} -> \pa{Float} -> \pa{Float}
\end{alltt}
The intermediate functions \ccode{\SE{(+)}} and \ccode{\LV{(+)}} deal with
partial applications of @(+)@.  Finally we reach ground truth: invocations of
\ccode{\SE{addFloat}} and \ccode{\LV{addFloat}}, which are implemented by
the back end.  The former is the ordinary floating point addition instruction; the latter is 
a ``vector instruction'', which will be implemented differently on different targets.
On a sequential machine it will be implemented as a loop; on a GPU it will be implemented
using vector hardware; on a cluster it will be implemented using a loop on each
CPU with barrier synchronisation at the end.   Section~\ref{sec:post-vect} elaborates.

These functions look grotesquely inefficient, especially considering how trivial
the original function @inc@ was.  Fortunately, most of the clutter is introduced
to account for the \emph{possibility} of higher order programming, and can
be removed by straightforward optimisations.

For example, consider the sub-term \ccode{\VV{(+)} \capp{}~x \capp{}~1} in
the definition of \ccode{\SE{inc}}. We can simplify it in the following way:
% $$
% \begin{array}{rl}
%  & \ccode{\VV{(+)} \capp{}~x \capp{}~1} \\
% \mbox{(Inline \ccode{\VV{(+)}})} & \implies 
%   \ccode{(Clo () \SE{(+)} \LV{(+)}) \capp{}~x \capp{}~1} \\
% \mbox{(Definition of \ccode{\capp})} & \implies 
%   \ccode{\SE{(+)} () x  \capp{}~1}  \\
% \mbox{(Inline \ccode{(+)\sub{S}})} & \implies 
%   \ccode{(Clo x \SE{addFloat} \LV{addFloat})  \capp{}~1} \\
% \mbox{(Definition of \ccode{\capp})} & \implies 
%   \ccode{\SE{addFloat} x 1}  \\
% \end{array}
% $$
\begin{alltt}
    \VV{(+)} \capp{} x \capp{} 1)
  \(\implies\sub{\textrm{Inline \ccode{\VV{(+)}}}}\)
    (Clo () \SE{(+)} \LV{(+)}) \capp{} x \capp{} 1
  \(\implies\sub{\textrm{Definition of \ccode{\capp}}}\)
    \SE{(+)} () x \capp{} 1
  \(\implies\sub{\textrm{Inline \ccode{\SE{(+)}}}}\)
    (Clo x \SE{addFloat} \LV{addFloat}) \capp{} 1
  \(\implies\sub{\textrm{Definition of \ccode{\capp}}}\)
    \SE{addFloat} x 1
\end{alltt}
All the intermediate closure data structures are removed.  (To save
generating huge intermediates during compilation, we are exploring
whether the vectorisation transformation could have special cases to
avoid introducing them in the first place.)

\subsection{The higher order example again} \label{sec:ho-vect}

It is instructive to see a case where the use of higher order functions
prevents complete removal of intermediate closures.  Let us return to the
@ho@ example of Section~\ref{sec:fun-pairs}:
\begin{alltt}
  ho :: (Int->Bool) -> (Bool, [:Bool:])
  ho f = (f 2, mapP f [:1,2,3:])
\end{alltt}
Again applying the vectorisation transformation blindly we get this:
\begin{alltt}
  \VV{ho} :: (Int :-> Bool) :-> (Bool, \pa{Bool})
  \VV{ho} = Clo () \SE{ho} \LV{ho}

  \SE{ho} :: () -> (Int :-> Bool) -> (Bool, \pa{Bool})
  \SE{ho} () f = (f \capp 2, \VV{mapP} \capp f \capp [:1,2,3:])

  \LV{ho} :: \pa{()} -> \pa{(Int :-> Bool)} -> \pa{(Bool, \pa{Bool})}
  \LV{ho} (ATup_0 n) fs 
    = \LV{(,)} (fs \LV{\capp} replicatePA n 2)
           (replicatePA n \VV{mapP} \LV{\capp} fs \LV{\capp} replicatePA n [:1,2,3:])  
\end{alltt}
We have taken a short-cut here by using optimised transformation rules for pairs:
$$\begin{array}{rcl} 
\VT{(e_1,e_2)}    & = & (\VT{e_1}, \VT{e_2}) \\
\LT{n}{(e_1,e_2)} & = & \LV{(,)}\;\;n\;\;(\LT{n}{e_1})\;\;(\LT{n}{e_2})
\end{array}
$$
The reader may verify the correctness of this 
optimised rule by seeing what happens instead if we
use the normal translation $\VV{\ccode{(,)}}\;\capp\;  \VT{e_1}\; \capp \; \VT{e_1}$,
and the definition of $\VV{\ccode{(,)}}$, which in turn is very like that for $\ccode{(+)}$.
Because of our array representation, the lifted pairing 
function $\LV{\ccode{(,)}}$ is a constant-time operation:
\begin{alltt}
  \LV{\ccode{(,)}} :: Int -> \pa{a} -> \pa{b} -> \pa{(a,b)}
  \LV{\ccode{(,)}} n xs ys = ATup\sub{2} n xs ys
\end{alltt}

\subsection{How flattening happens}

In our informal overview (Section~\ref{sec:vect-overview}) we said that we
``replace a call (@mapP f@) by \ccode{\LV{f}}''.   Higher order flattening
takes that \emph{static} decision and makes it \emph{dynamic}, by representing @f@ by a pair of
functions, thereby allowing @mapP@ to select at runtime.  (With the usual compile-time
optimisations when @f@ is known, of course.)  The code for @mapP@ itself is
therefore the heart of the way in which nested data parallelism is transformed
to flat data parallelism.  Here it is:
\begin{alltt}
  \VV{mapP} :: (a :-> b) :-> \pa{a} :-> \pa{b}
  \VV{mapP} = Clo () mapP\sub{1} mapP\sub{2}
  
  mapP\sub{1} :: () -> (a :-> b) -> \pa{a} :-> \pa{b}
  mapP\sub{1} _ f = Clo f \SE{mapP} \LV{mapP}

  mapP\sub{2} :: \pa{()} -> \pa{(a :-> b)} -> \pa{(\pa{a} :-> \pa{b})}
  mapP\sub{2} _ fs = AClo fs \SE{mapP} \LV{mapP}

  \SE{mapP} :: (a :-> b) -> \pa{a} -> \pa{b}
  \SE{mapP} (Clo env fs fl) xss 
    = fl (replicatePA (lengthPA xss) env) xss

  \LV{mapP} :: \pa{(a :-> b)} -> \pa{\paB{a}} -> \pa{\paB{b}}
  \LV{mapP} (AClo env _ fl) xss
    = unconcatPA xss (fl (expandPA xss env) (concatPA xss))
        -- xss :: \pa{\paB{a}}
        -- env :: \pa{e}
        -- fl  :: \pa{e} -> \pa{a} -> \pa{b}
\end{alltt}
In \SE{@mapP@} we exploit the key observation from
Section~\ref{sec:vect-overview}, namely that we can define the doubly-lifted
function using the singly-lifted one @fl@, using constant-time reshaping
operations on the data.  Unfortunately, to account for free variables, we face
a small complication: the environment \ccode{env} contains one element for
each \emph{subarray} of \ccode{xss}. Thus, before applying \ccode{fl} we must
expand \ccode{env}, i.e., repeat each element as many times as the
corresponding subarray of \ccode{xss} has elements.  For top-level functions,
the environment will be empty and @expandPA@ performs no work.

Of course, @mapP@ is not the only function that the library must implement.
All of (the @PA@ versions of) the functions in Figure~\ref{fig:parr-ops} must be
provided in vectorised form. For example, here is the lifted version of @zipP@
(the definition of @zipPA@ is given in Section~\ref{sec:arr-alg}):
\begin{alltt}
  zipP\sub{L} :: \pa{\paB{a}} -> \pa{\paB{b}} -> \pa{\paB{a,b}}
  zipP\sub{L} (AArr segd xs) (AArr \_ ys) = AArr segd (zipPA xs ys)
\end{alltt}
These library functions are the heart of flattening: they make nested data 
parallelism ``go''.  Everything is
organised to make their implementation, especially their lifted variants,
work efficiently.

\subsection{User-defined data types}\label{sec:user-types}

One of Haskell's strengths is the ease with which programmers can declare new
algebraic data types, and process them using pattern matching.  DPH allows all of
this expressiveness in fully-vectorised code as well.
There are two main complications: occurrences of @(->)@ and @[::]@ in
user-defined data types; and representing arrays of values drawn from 
such types.  We discuss each in turn.

\subsubsection{Vectorising user-defined data types} \label{sec:user-DT}

In Figure~\ref{fig:vect}, the type transform $\TT{\tau}$ replaces a user-defined
data type $T$ by its vectorised counterpart $\VV{T}$.  But what exactly is $\VV{T}$?
Consider
\begin{code}
  data Fun = MkFun (Int -> Int) 
\end{code}
Remember that in the vectorised program, each function arrow @(->)@ must
be replaced by a function closure @(:->)@ --- and of course that must also
happen 
inside data types.  So we must generate a vectorised version of @Fun@, thus:
\begin{alltt}
  data \VV{Fun} = \VV{MkFun} (Int :-> Int)
\end{alltt}
This must be done recursively: if a constructor of data type @T@ mentions @Fun@, 
then @T@ too must have a vectorised version. So the vectorised variant of each
data type obtained by simply applying the $\TT{}$ transform to every type in the
data type declaration.  The @Tree@ type of Section~\ref{sec:bh-dph} 
is another good example, because we must replace @[::]@ by \ccode{\pacon}:
\begin{alltt}
  data \VV{Tree} = \VV{Node} Mass Location \paB{Tree}
\end{alltt}

While we \emph{can} generate a vectorised version of every data type,
it is \emph{unnecessary} to do so for data types that do not mention
functions or parallel arrays.  Happily, almost all data types fall
into this category; for example @Bool@, @Maybe@, lists, tuples, and so
on.  We quietly took advantage of this in the @Tree@ example, by not
transforming @Mass@ to \ccode{\VV{Mass}} (and similarly @Location@)
because @Mass@ = \ccode{\VV{Mass}}.
In Section~\ref{sec:selective} we will see a second reason to avoid
vectorising a data type unless it is absolutely necessary to do so.


\subsubsection{Arrays of user-defined data types}

The ideas of Section~\ref{sec:arr-alg} can readily be extended to work for
arbitrary user-defined algebraic data types. We have already seen how this
works for \ccode{Tree} in Section~\ref{sec:rec-types}. Here is another
example, a sum type:
\begin{code}
  data Maybe a = Nothing | Just a
\end{code}
How can we represent an array of @Maybe Float@ values?  The natural
dense representation is as a pair of (a) an array of booleans (@True@
for @Nothing@, and @False@ for @Just@), and (b) an array of @Float@
containing only the @Just@ values:
\begin{alltt}
  instance \pacls a => \pacls (Maybe a) where
    data \pa{Maybe a} = AMaybe \paB{Bool} \paB{a}
    indexPA (AMaybe bs vs) i 
      | indexPA bs i = Nothing
      | otherwise    = indexPA vs (indexPA just\_indices i)
      where
        just\_indices = scanPA (+) 0 (mapPA boolToInt bs)
    lengthPA (AMaybe bs _) = lengthPA bs
\end{alltt}
In practice, to avoid computing @just_indices@ on each indexing operation we precompute the
index vector, and cache it in an extra field of the @AMaybe@ constructor.

In our real implementation, we avoid generating a big instance declaration for 
every such user-defined data type, by instead generating code to convert it
to a simple sum-of-products representation, and then using a set of fixed 
instances for \ccode{\pacls} at those representation types.


\subsection{What we have swept under the carpet}

Vectorisation is a complicated transformation, and to keep it comprehensible
we have simplified several aspects.  In this section we briefly mention some of them.

\subsubsection{Types and dictionaries}

The alert and Haskell-savvy reader will have noticed the following
discontinuity in our presentation.  We described \ccode{\pacon} type
in association with a \emph{type class}, \ccode{\pacls}.  However, type classes
are dealt with by the type inference system, right at the front end of
the compiler, and are completely translated out in the passage to the
Core intermediate language. In this desugaring, a function with an overloaded
type, such as @nub :: Eq a => [a] -> [a]@ is given a second parameter
which is a record, or ``dictionary'', of the functions that implement the
operations of the @Eq@ class.

In the desugared program, @nub@ has type @EqD a -> [a] -> [a]@, where
@EqD@ is an ordinary data type, thus:
\begin{alltt}
  data EqD a = EqD \{ (==) :: a -> a -> Bool
                   , (/=) :: a -> a -> Bool \}
\end{alltt}
Correspondingly, the desugarer injects an extra argument at every call to @nub@,
namely the correct method suite for that particular call site.

The vectoriser generates many calls to 
@replicatePA@, @splitPA@, etc, which have type-class-constrained types,
\emph{yet the vectoriser runs after typechecking and desugaring are complete}.
So the vectoriser cannot take advantage of the implicit injection of extra arguments;
instead it must insert them itself.
In the real implementation of Figure~\ref{fig:vect}, the
vectoriser therefore adds appropriate dictionary abstractions and
applications.  (In fact, since GHC's Core language is an
explicitly-typed variant of System F, we also inject type abstractions
and applications.)  All this is tiresome but routine; showing
the implicit abstractions and applications 
in Figure~\ref{fig:vect} would have dramatically obfuscated an
already-dense figure.

\subsubsection{Selective vectorisation} \label{sec:selective}

As mentioned earlier, we do not really vectorise the \emph{whole}
program; rather, we selectively vectorise parts of it.  We must also
generate marshaling code to allow us to ``cross the border'' between
vectorised and unvectorised code.  For example, in Barnes-Hut, we
presumably want to vectorise the @oneStep@ function, which will give
us
\begin{alltt}
  \VV{oneStep} :: \pa{Particle} :-> \pa{Particle}
\end{alltt}
If we want to be able to call @oneStep@ from ordinary scalar code, 
we must generate the following marshalling code:
\begin{alltt}
  oneStep :: [:Particle:] -> [:Particle:]
  oneStep ps = fromPA (\VV{oneStep} \capp (toPA ps))

  toPA   :: \pacls a => [:a:] -> PA a
  fromPA :: \pacls a => PA a -> [:a:]
\end{alltt}
Marshaling may also be necessary for user-defined data types..  
For example, suppose we vectorise a function @f :: Int -> Fun@, 
so that \ccode{\VV{f} ::\ Int :-> \VV{Fun}} (cf. Section~\ref{sec:user-types}
for the definition of @Fun@).  If we want to call @f@ from
normal scalar code, we must generate: 
\begin{alltt}
  f :: Int -> Fun
  f n = case \VV{f} n of
          \VV{MkFun} tf -> MkFun ((\capp) tf)
\end{alltt}
Of course, it gets worse if the data type is recursive, because the marshaling
code has to traverse the whole structure.  On the other hand, no marshaling is
needed for types that have have no functions or arrays inside them, which
is a strong reason for exploiting that special case (Section~\ref{sec:user-DT}).

Marshaling has a run-time cost.
In particular, the calls to @toPA@ and @fromPA@ change the data representation for 
parallel arrays, and so are potentially \emph{very} expensive. In fact, it is
possible to choose a representation for @[:a:]@ that mitigates these costs
somewhat but in general, marshaling data across the border should be avoided.  

The question of just which parts of the program to vectorise is therefore
an interesting one.  We want to vectorise code that can run in parallel; we want
to reduce marshaling to a minimum; and we do not want to vectorise code where there
is little or no benefit.  We suggest automatically<ic approaches in \cite{chakravarty-et-al:partial-vect},
but it may also be reasonable to seek help from the programmer (e.g. ``vectorise module
X but not module Y'').

\subsubsection{Vectorisation of Scalar Functions}

For some scalar functions, it is not necessary to apply the vectorisation transformation to get a flat data parallel version, as a simple parallel map would do the job. For example, the increment function we used to demonstrate the effect of lifting earlier:
\begin{code}
  f :: Float -> Float
  f x = x*x + 1
\end{code}
and it's lifted version:
\begin{alltt}
  \LV{f} :: [:Float:] -> [:Float:]
  \LV{f} x = (x \LV{*} x) \LV{+} (replicateP n 1)
         where
           n = lengthP x
\end{alltt}
But, as the argument as well as the result type is a type whose vectorised version is identical to the original type and there are no parallel computations in the function, we could have just  used parallel map to get the liften version:
\begin{alltt}
  \LV{f} :: [:Float:] -> [:Float:]
  \LV{f} = mapP (\hslam z. z + 1) x
\end{alltt}
This definition is much simpler, and also more efficient, since we have a single parallel loop with no intermediate data structure (as would be produced by replicate). Now, we use fusion aggressively to address the problem of additional traversals and intermediate structures, but  there is no much point in decomposing the functions, replacing every operation with the parallel version, and then trying to reassemble it via fusion into a single parallel loop if we can statically identify some of these functions, and avoid the transformation and need for successive optimisations completely.


So, the rules for lambda-expressions are as follows

$$
\begin{array}{l}
\begin{array}{rll}
\VT{\lambda{x}.e}   
  & \multicolumn{2}{l}{
              | isScalar (\lambda{x}.e)   = @Clo@\begin{array}[t]{l} <
              \{ @env@ = (y_1,\dots,y_k) \\
              , @clo@_s = \lambda e\,x.\,
                 @case@\,e\,@of@\;(y_1,\dots,y_k) \to e \\
              , @clo@_l = \lambda e\,x.\,
                 @case@\, e\, @of@\; @ATup@_k\; n\; y_1 \dots y_k \to @map@\, \lambda{x}.{e} \}
           \end{array}} \\
  & \multicolumn{2}{l}{
              | otherwise  = @Clo@\begin{array}[t]{l} <
              \{ @env@ = (y_1,\dots,y_k) \\
              , @clo@_s = \lambda e\,x.\,
                 @case@\,e\,@of@\;(y_1,\dots,y_k) \to \VT{e} \\
              , @clo@_l = \lambda e\,x.\,
                 @case@\, e\, @of@\; @ATup@_k\; n\; y_1 \dots y_k \to \LT{n}{e} \}
           \end{array}} \\
 & \multicolumn{2}{l}{\mbox{where} \{y_1, \ldots, y_k\} = \mbox{free variables of $\lambda x.e$}} \\
\end{array}
\\ \\
% \begin{array}{rll}
% 	isScalar (\lambda{x}.e :: \tau_1 \to \tau_2) = &\VT{\tau_1}  == \tau_1\; and\; \VT{\tau_2} == \tau_2 &\mbox{Primitive scalar type}\\
% 	&
% \end{array}
\end{array}
$$
where \textit{isScalar} holds if the following conditions are true:
\begin{itemize}
	\item both the argument and result types are of primitive scalar type
	\item subexpression of \textit{e} cannot be of array type or compound type which may contain an array type. This excludes values abstract type, as we can't guarantee that there might not be a hidden array type. User defined types are  not a problem - we traverse the type and make sure that it can't possibly contain an array type.
	\item any identifier occurring in \textit{e} must either refer to a  global scalar value or function, a locally bound value, or a scalar, possibly user defined data constructor.
	\item the expression \textit{e} has to contain a reference to a a global scalar. This requirement rules out functions which don't do anything. For example, the identity function gets lifted to the identity function over arrays, which is more efficient than mapping the identity function over scalars over an array.
\end{itemize}

The way \textit{isScalar} is currently implemented means that a lambda-expression is sometimes lifted, even if it's not actually necessary, for example, if it contains a call to another scalar function. To catch these cases, we need to make information about values accessible, so that we can identify all scalar values, not only the fixed set of global scalar functions.

\subsubsection{Laziness}

Consider this function:
\begin{alltt}
  f :: Int -> Int
  f x = h x (1/x)
\end{alltt}
Although @x@ might be zero, let us assume that @h@ only evaluates its second argument
if its first argument is non-zero.  Haskell's lazy evaluation therefore ensures that
no divide-by-zero exception is raised.

The lifted version will look something like this:
\begin{alltt}
  \LV{f} :: \pa{Int} -> \pa{Int} -> \pa{Int}
  \LV{f} xs = \LV{h} xs (replicatePA (lengthPA xs) 1 \LV{/} xs)
\end{alltt}
The trouble is that a demand for \emph{any} element of \ccode{\LV{h}}'s second
argument will force \emph{all} the elements to be evaluated, including the 
divisions by zero.  Something very similar arises in a more local context when
we have @let@ expressions:
\begin{alltt}
  f x = let y = 1/x in
        if x==0 then 0 else y+1
\end{alltt}
Although this is something of a corner case, we do not yet have a very
satisfying solution.  We currently simply ignore the problem, and accept
the slight change in semantics.  A better solution might be to reify
the exception into an exceptional value (like a IEEE NaN); but that
carries an efficiency cost.  Lastly, we might treat the argument as a
nullary function, accepting the loss of sharing that would result.



% ---------------------------------------------------------------------

\section{Multicore execution model} \label{sec:post-vect}

The vectorisation transformation turns all nested data parallelism into
parallel operations on flat arrays, as used by the instances of the @PAElem@
class.  The transformation is crucial to express parallel algorithms on a
high-level of abstraction and in a modular fashion.  However, purely
functional array operations, even if restricted to flat arrays, are still a
far cry from the hardware model of multicore CPUs.

In particular, the vectorised code uses many superfluous intermediate arrays,
which increase the overhead of memory management and whose creation involves
extra synchronisation between parallel CPUs.  Even worse, the repeated
traversal of large structures compromises locality of reference, and so, has a
very negative effect on execution performance.  Finally, we need to map the
data parallel array operations onto the multi-threaded execution model of
multicore CPUs by way of the Single Program Multiple Data (SPMD) model \cite{darema:SPMD}.

In contrast to the vectorisation transformation, we can implement the mapping
from flat data-parallel code to SPMD code using existing transformation and
optimisation phases of GHC; in particular, we make heavy use of GHC's inliner
and rewrite rules~\cite{ghc-inliner,pj-et-al:playing-by-the-rules}, which enable
library-specific optimisations as part of library source code, in the form of
compiler pragmas.  Consequently, we can implement these transformations
without altering GHC's source code, which greatly simplifies experimentation
with different transformations.

In the reminder of this section, we illustrate the transformation of flat
data-parallel code into SPMD code by way of an example.  Further details are
in~\cite{keller-chak:dtypes,chak-keller:afp-arrays,CSL06,CLS07,chak-etal:DPH}.

\subsection{Running example}

As an example, we consider the computation of the value @far_forces@, in the
function @calcForces@ of Section~\ref{sec:bh-dph}, by way of the array
comprehension,
%
\begin{code}
  [: forceOn p m l | p <- ps, isFar len l p :]
\end{code}
%
After vectorisation and simplification to remove intermediate closure data
structures, we have
%
\begin{alltt}
  \LV{forceOn'} (\SE{filterP} (isFar len l) ps) m l
\end{alltt}
%
The code performs two collective operations on the input array @ps@ in
sequence.  Firstly, the application of \SE{@filterP@} to remove all particles
that are not far, and secondly, a computation that corresponds to lifting
@forceOn@ only on its first argument (here called, \LV{@forceOn'@}):
\begin{alltt}
  \LV{forceOn'} :: \pa{Particle} -> Mass -> Location -> \pa{Force}
\end{alltt}
Such pipelines of collective operations are typical for data-parallel code.

\subsection{The SPMD execution model}

The implementation of collective array operations, such as @mapP@ and
@filterP@, needs to distribute the workload evenly across the the available
processing elements (PEs), such as multiple cores and CPUs.  In the
data-parallel model, the workload of a PE is dependent on the number of array
elements residing on that PE.  Hence, we balance work by suitably distributing
the array elements.  By default, we choose an even distribution; i.e., given
$p$ PEs and an array of length $n$, each PE gets about $n/p$ array elements.

In the SPMD model, the individual PEs process local array
elements until they arrive at a point in the computation where they require
non-local data, and need to cooperate with other PEs.  In
our example, the result of \SE{@filterP@} is such a point.  Even if the
\emph{input} to \SE{@filterP@} is an array that is evenly spread across the
PEs, the \emph{output} of \SE{@filterP@} might be wildly unbalanced,
depending on which elements of the array are selected by the predicate.  If
so, any further processing of that array would have an equally unbalanced work
distribution.

To avoid a work imbalance, arrays need to be re-distributed when their size
changes.  Redistribution is a cooperative process in which all PEs need to
coordinate.  However, redistribution is not the only such operation in an SPMD
implementation of data parallelism.  Other prominent cooperative operations
are reductions (such as @foldP@), pre-scans (such as @scanlP@), and
permutation operations.  Overall, a parallel program executing in SPMD-style
alternates between \emph{processing phases,} where the PEs operate
independently on local data, and \emph{communication phases,} where the
processing elements interact and exchange data.

Communication phases are typically expensive because they include data
exchange and blocking to allow any slower PEs to catch up.  Hence, compiler
optimisations that remove communication phases in favour of longer-running
processing phases are often worthwhile.  In particular, the redistribution of
arrays after operations that change the array length, such as \SE{@filterP@},
does not necessarily improve overall runtime.  An inexpensive, purely local
operation may be faster, even if work is not ideally balanced, than an
expensive redistribution followed by the same local operation with a perfectly
balanced workload.

\subsection{Gang parallelism}

Our implementation of the SPMD model for data parallelism is based on the
coordinated execution of a \emph{gang of threads,} with one thread per PE.
GHC includes a Haskell library for concurrent programming with explicit
thread forking and thread communication primitives.  It forms the lowest level
of abstraction in our data-parallel array framework and enables us to implement
the entire library in Haskell without any special compiler support or the
need to resort to C code.

We need to make the distributed nature of computations in the SPMD model
explicit to further compile the code resulting from vectorisation, such as
%
\begin{alltt}
  \LV{forceOn'} (\SE{filterP} (isFar len l) ps) m l
\end{alltt}
%
In this context, distribution does not imply that the data is necessarily
located on physically distinct memory banks, but that different threads are
responsible for the processing of different portions of parallel arrays.  By
being explicit about distribution, we are automatically also explicit about
the distinction of processing phases versus communication phases.

Our main vehicle for distinguishing between these two phases and making
distribution explicit is the type @Dist a@ of \emph{distributed values.} For
instance, @Dist Int@, pronounced ``distributed @Int@'', denotes a collection
of \emph{local} integers, such that there is one local integer value per gang
thread.  Arrays can be distributed, too: @Dist [:Float:]@ is a collection of
local array \emph{chunks}, again one per gang member, which together make up
the array.  Arrays are distributed across gang members and joined back
together by the following functions:
%
\begin{quote}
\begin{code}
splitD :: PA a -> Dist (PA a)
joinD  :: Dist (PA a) -> PA a
\end{code}
\end{quote}
%
Distributed values support a number of operations, most importantly mapping:
%
\begin{quote}
\begin{code}
mapD :: (a -> b) -> Dist a -> Dist b
\end{code}
\end{quote}
%
While @splitD@ and @joinD@ denotes communication, @mapD@ is the main means of
implementing parallel processing phases: the gang members
concurrently apply the (purely sequential) function to their respective local
values. 

\subsection{Inlining of gang code}

Given a sequential filter function operating on a \emph{single chunk} of a
parallel array
%
\begin{quote}
\begin{alltt}
filterS :: (a -> Bool) -> PA a -> PA b
\end{alltt}
\end{quote}
%
we can define \ccode{\SE{filterP}} as a distributed gang computation as follows:
%
\begin{quote}
\begin{alltt}
\SE{filterP} p arr = joinD (mapD (filterS p) (splitD arr))
\end{alltt}
\end{quote}
%
Given a global parallel array, which is not distributed, @splitD@ distributes
the array across the gang, @mapD (filterS p)@ applies the sequential filter
function in parallel to all chunks of the distributed array, and finally,
@joinD@ combines the various chunks, which may now be of varying length, into
one global array.

Similarly, \LV{@forceOn'@} internally consists of @mapD@s that compute the force
for each particle.  The force computations for the individual particles are
entirely independent, so we can assume \LV{@forceOn'@} to have the following
structure:
%
\begin{quote}
\begin{alltt}
\LV{forceOn'} ps m l 
  = joinD (mapD (mapS (\hslam p. forceOn p m l)) (splitD ps))
\end{alltt}
\end{quote}
%
where @forceOn@ is the original, sequential function from the source of our
Barnes-Hut implementation and @mapS@ is a purely sequential array mapping
function.

GHC's inliner will inline the definition of both \SE{@filterP@} and
\LV{@forceOn'@}; i.e., it will perform the following rewriting:
%
\begin{alltt}
    \LV{forceOn'} (\SE{filterP} (isFar len l) ps) m l
  \(\Longrightarrow\sub{\textrm{Inlining}}\)
    joinD (mapD (mapS (\hslam p. forceOn p m l)) (
      splitD (joinD (mapD (filterS (isFar len l)) (splitD ps)))))
\end{alltt}

Of special interest here is the function @splitD@ which is applied to the
immediate result of
@joinD@ (in the second line of the resulting expression).  This turns a
distributed array into a global array and distributes it again.  In contrast
to the original array, the newly distributed one is guaranteed to be
distributed evenly; hence, a @splitD/joinD@ combination performs load
balancing. 

However, as we remarked earlier, it is often an advantage to accept some load
imbalance in favour of avoiding communication phases in an SPMD computation.
In GHC, we easily achieve that by specifying the following \emph{rewrite rule}:
%
\begin{quote}
\begin{code}
"splitD/joinD" forall xs. splitD (joinD xs) = xs
\end{code}
\end{quote}
%
GHC has support for specifying such rewrite rules directly in the library source
code as compiler pragmas~\cite{pj-et-al:playing-by-the-rules}.  Applications
of the @splitD/joinD@ rule frequently produce two adjacent applications of
@mapD@, which signal two adjacent purely sequential and thread-local
computations.  We can combine them, and hence eliminate a synchronisation
point, using the well known map fusion law:
%
\begin{quote}
\begin{code}
"mapD/mapD" forall f g xs. 
  mapD f (mapD g xs) = mapD (f . g) xs

\end{code}
\end{quote}

Applying both rules to our example, we get
%
\begin{alltt}
    joinD (mapD (mapS (\hslam p. forceOn p m l)) (
      splitD (joinD (mapD (filterS (isFar len l)) (splitD ps)))))
  \(\Longrightarrow\sub{\textrm{Apply \texttt{splitD/joinD}}}\)
    joinD (mapD (mapS (\hslam p. forceOn p m l)) (
      mapD (filterS (isFar len l)) (splitD ps)))
  \(\Longrightarrow\sub{\textrm{Apply \texttt{mapD/mapD}}}\)
    joinD (mapD (mapS (\hslam p. forceOn p m l) . filterS (isFar len l))
                (splitD ps)))
\end{alltt}
%
At this point, the question arises whether we can combine adjacent
sequential array combinators, such as @mapS@ and @filterS@, to reduce the
number of array traversals and intermediate data structures.  Indeed, we
aggressively remove such inefficiencies using a fusion framework known as
\emph{stream fusion}~\cite{CSL06,CLS07,chak-etal:DPH}, but we will refrain
from discussing this in detail.

This concludes our brief overview of the post-vectorisation aspects
of Data Parallel Haskell.  A somewhat more detailed discussion can be found
in \cite{chak-etal:DPH}.

\section{Related work}

We discussed prior work on the implementation of language support
for nested data parallelism in detail in~\cite{chak-etal:DPH}. In this
paper, we will only give a brief overview of existing work, and how
they compare to our approach. 

The starting point for our work was the nested data parallel programming model
of NESL~\cite{blelloch:diss,blelloch-etal:impl-nested-parallel}, which we
extended and implemented in the context of a general-purpose language and GHC,
a state-of-the-art compiler. Consequently, we have to deal with a multitude of
issues not previously addressed, as for example the combination of
user-defined and parallel data structures, selective vectorisation,
higher-order functions, separate compilation, and aggressive cross-function
optimisation.

Prins et al.\ worked on various aspects of the vectorisation of nested data
parallel programs; see,
e.g.,~\cite{prins-etal:nested-trafo,palmer-etal:efficient-trafo}.  Most of
their work was also in the context of a functional language, but one that like
NESL lacks many of Haskell's features.  Their work is largely orthogonal to
ours.

The Proteus system~\cite{mills-etal:proteus} promised a combination of data
and control parallelism, but Proteus had a particular focus on manual
refinement of algorithms, where data parallel components were automatically
vectorised, this again was a complete whole-program transformation.  Moreover,
the system was never fully implemented.

Manticore~\cite{fluet-etal:manticore-status} supports a range of forms
of parallelism including nested data parallelism.  Manticore employs
some of the same techniques that we use, but does not implement
flattening yet~\cite{fluet-etal:implicit-threaded}.  According
to the project web page, a preliminary implementation of the Manticore
system should be available around the time when this paper is published.

So et al.~\cite{so-etal:dp-many-core} developed a parallel library of
immutable arrays for C/C++ supporting what they call \emph{sub-primitive
  fusion}.  Their choice of immutable arrays, despite working with imperative
languages, is to enable aggressive program transformations, much like in our
approach.  However, where we apply transformations statically during compile
time, their library builds a representation of the to be executed computation
at runtime.  Consequently, they require less compiler support and do not have
to worry about inlining and similar optimisations.  However, they incur a
runtime penalty by performing optimisations at runtime and need to amortise
that penalty by further optimisations.  Like us, they also strive for a
seamless integration of data parallelism and explicit concurrency within a
single program.


\section{Conclusion}

We are excited about Data Parallel Haskell because it gives us 
some chance of writing parallel programs that can in principle
efficiently exploit very large parallel machines working on large data sets.

In this paper we have outlined solutions to the challenges of 
polymorphism, higher order functions, and user-defined data types.
There is much to do, however, before we can declare victory.
The very generality of Data Parallel Haskell makes it an ambitious undertaking.
Many components have to work together smoothly to generate efficient code
--- and that is before we start to consider matters such as using SSE vector
instructions or GPUs, or mapping to a distributed memory architecture.
Nevertheless, we regard nested data parallelism general, and Data Parallel Haskell
in particular, as a very promising and exciting approach to harnessing the multicores.

\section*{Acknowledgements}

We gratefully acknowledge the help of Max Bolingbroke, 
Ryan Ingram and John Reppy, whose comments
led to real improvements in the paper.  Thank you.


\bibliographystyle{alpha}
\bibliography{fsttcs2008,slpj-author,simon}
\end{document}
