\documentclass[preprint]{sigplanconf}

%\usepackage[british]{babel}
\usepackage{version}
\usepackage[pdftex]{graphicx}
\usepackage{amsmath}
\usepackage{amsfonts,amssymb}
\usepackage{url}
\usepackage{alltt}
\usepackage{code}


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

\newcommand{\simon}[1]{\emph{\textbf{SLPJ:} #1}}
\newcommand{\slpj}[1]{\emph{\textbf{SLPJ:} #1}}
\newcommand{\rl}[1]{\emph{\textbf{RL:} #1}}
\newcommand{\gck}[1]{\emph{\textbf{GCK:} #1}}
\newcommand{\manuel}[1]{\emph{\textbf{MMTC:} #1}}
\newcommand{\benl}[1]{\emph{\textbf{BL:} #1}}
\newcommand{\TODO}[1]{\emph{\textbf{TODO:} #1}}

% \makeatletter
% \def \ps@plain {%
%   \let \@mkboth = \@gobbletwo
%   \let \@evenhead = \@empty
%   \def \@evenfoot { \thepage \hfil \scriptsize \textit{\@preprintfooter}\hfil
%                    \textit{\@formatyear}}%
%   \let \@oddhead = \@empty
%   \let \@oddfoot = \@evenfoot}
% \makeatother

\begin{document}

\preprintfooter{\textbf{--- DRAFT --- DRAFT --- DRAFT --- DRAFT ---}}

%\title{Parallel Divide-and-conquer Algorithms for Data-intensive Problems%
\title{Parallel Divide-and-conquer Algorithms as Array Code in Haskell%
  \\--- DRAFT ---}
\authorinfo{
  \shortstack{%
    Ben$^\dagger$ 
    \and
    Gabriele$^\dagger$ 
    \and 
    Manuel$^\dagger$ 
    \and 
    Roman$^\dagger$ 
    \and
    Simon$^\ddagger$ 
%    \\[4pt]
%    -- second line
  }
}{
  \vspace{5pt}
  \shortstack{
    $^\dagger$Computer Science and Engineering \\
    University of New South Wales, Australia \\[2pt]
    \textsf{\{keller,chak,rl,benl\}@cse.unsw.edu.au}
  }
  \and
  \shortstack{
    $^\ddagger$Microsoft Research Ltd \\
    Cambridge, England \\[2pt]
    \textsf{\{simonpj\}@microsoft.com}
  }
}{%
}
\maketitle

\makeatactive

\begin{abstract}
  Divide-and-conquer algorithms are an attractive category of algorithms for parallel programming.  However, the performance and scalability of list and tree-based versions is limited, especially for data-intensive problems.  Array-based versions have the potential to be more efficient, but are often harder to implement and parallelise.  We illustrate at \TODO{HOWMANY} examples that \emph{nested data parallelism}, as implemented in Data Parallel Haskell, might provide a middle ground.  It supports a purely functional, list-like array abstraction and implicit parallelism for all array operations.  We compare both the source code as well as the sequential and parallel performance of divide-and-conquer algorithms implemented with lists, the vector package, and Data Parallel Haskell.
\end{abstract}


\section{Introduction} % (fold)
\label{sec:introduction}

Divide-and-conquer algorithms are obvious candidates for parallel programming.  They split a large problem into a number of smaller subproblems that can be solved independently, possibly in parallel.  Nevertheless, the effective parallelisation of these algorithms comes with its own set of challenges.  Some of these challenges are related to load balancing and granularity: problems may be split into subproblems of varying complexity and recursive problem decomposition can lead to a huge number of fine-grained work items.  Other challenges are related to the splitting of a problem into subproblems and the merging of subsolutions into overall solutions.  Both tasks can themselves be rather work intensive and require parallelisation to achieve good overall parallel scalability --- Amdahl's law tells us that even a small sequential fraction of the overall runtime can seriously limit scalability~\cite{amdahl}.

Haskell, and specifically the Glasgow Haskell Compiler (GHC), has grown three major extension to facilitate parallel programming: control parallel annotations~\cite{trinder-etal:GUM}, explicit multi-threading~\cite{peyton-jones:concurrent-haskell,harris-etal:haskell-shared-memory}, and nested data parallelism~\cite{blelloch:ppa,peyton-jones-etal:harnessing-multicores}.  All three extensions come with their own trade offs and their development is at different stages of maturity.

In this paper, we investigate the implementation of divide-and-conquer algorithms in all three extensions for parallel programming.  We compare implementation trade offs as well as sequential and parallel performance.  Our focus is on algorithms that solve data-intensive problems --- that is, problems where we have to process large datasets, but the number of operations on individual data items are comparatively few.  


In summary, we make the following contributions:
%
\begin{itemize}
  \item We present and compare methods to parallelise divide-and-conquer algorithms in all three of the parallel programming frameworks supported in the Glasgow Haskell Compiler, that is parallel annotations (Section~\ref{sec:parallelising_quickhull}), explicit multi-threading (Section~\ref{sec:an_imperative_version}),
   and Data Parallel Haskell (Section~\ref{sec:nested_data_parallelism}).
  \item We implement and discuss a parallel divide-and-conquer skeleton.  \TODO{Will we do that?}
  \item We present a quantitative evaluation of the various approaches on two different multicore architectures (Section~\ref{sec:benchmarks}).
\end{itemize}

% section introduction (end)


\section{Parallelising Quickhull} % (fold)
\label{sec:parallelising_quickhull}

As an example of fairly simple, but practically useful divide-and-conquer algorithm, let us consider Quickhull~\cite{barber-etal:quickhull} ---an algorithm to determine the planar convex hull of a set of points--- and it's parallelisation.  Quickhull gets its name because of its similarity to the Quicksort algorithm.   In line with Quicksort, the algorithmic pattern is to (1) pick a \emph{pivot} element, (2) split the data set based on the pivot, and (3) recurse on each of the split sets.  Also as with Quicksort, the pivot element is not guaranteed to split the data set into equally sized sets, and in the worst case the algorithm requires \(\mathcal{O}(n^2)\) work. Nevertheless, in practice, the algorithm is often very efficient.

The canonical approach to parallelising divide-and-conquer algorithms is to execute the two (or more) recursive calls in parallel.  This suggests the use of \emph{control parallelism} to execute the recursive calls in separate threads of control.  In the reminder of this section, after briefly illustrating the Quickhull algorithm, we shall discuss this form of parallelising the algorithm, using the @Control.Parallel.par@ combinator.  We will also consider a variant of the algorithm that uses a balanced tree structure instead of lists as the central data structure.

\TODO{Is there literature on turning d-n-c algorithms into DP algorithms other than that in the context of NESL (and Prins' paper(s) on irregular parallelism in Fortran)?}

\subsection{QuickHull in a nutshell} % (fold)
\label{sub:quickhull_in_a_nutshell}

\begin{figure*}
\begin{alltt}
upper :: (a -> a -> Bool) -> [(a, b)] -> b
upper above = snd . foldl1 pick
  where
    pick left@(kl, _) right@(kr, _) | kl `above` kr = left
                                    | otherwise     = right

distance :: Point -> Line -> Double
distance (Point xo yo) (Line (Point x1 y1) (Point x2 y2))
  = (x1-xo) * (y2 - yo) - (y1 - yo) * (x2 - xo)

hsplit :: [Point] -> Line -> [Point]
hsplit points line@(Line p1 p2)
  | length packed < 2 = p1:packed
  | otherwise         = hsplit packed (Line p1 pm) ++ hsplit packed (Line pm p2)
  where
    cross  = [ (distance p line, p) | p <- points ]
    packed = [ p | (p, (c, _)) <- zip points cross, c > 0.0 ]

    pm     = upper (>) cross

quickHull :: [Point] -> [Point]
quickHull [] = []
quickHull points
  = hsplit points (Line minx maxx) ++ hsplit points (Line maxx minx)
  where
    xs   = [ (x, p) | p@(Point x y) <- points ]
    minx = upper (<) xs
    maxx = upper (>) xs
\end{alltt}
\caption{The Quickhull algorithm in plain (sequential) Haskell} 
\label{fig:quickhull}
\end{figure*}
%
Figure~\ref{fig:quickhull} displays a plain Haskell version of Quickhull using lists.
\TODO{Explain the algorithm using a small example and a figure that graphical illustrates the splitting and how the remaining points form the convex hull.}

% subsection quickhull_in_a_nutshell (end)

\subsection{Annotations} % (fold)
\label{sub:control_parallelism}

The standard strategy to parallelise divide-and-conquer algorithms is to execute the two or more recursive calls simultaneously.  In the case of our Haskell implementation of Quickhull, we might like to execute the two calls to @hsplit@ in @quickHull@ and the recursive calls in @hsplit@ in parallel.  The module @Control.Parallel@ provides a function @par :: a -> b -> b@, which sparks the evaluation of its first argument before evaluating, and returning, the second.  In other words, the expression @e1 `par` e2@ yields the same result as just @e2@, but creates in addition a \emph{spark} for @e1@.  A spark is like a marker that identifies a subcomputation that \emph{may} be executed in parallel to the main evaluation.  Parallel evaluation is not compulsory; it is merely a hint to the runtime system where it can find worthwhile parallel work if idle resource are available.  Sparks a very lightweight and it is usually a good idea to create many.  However, it is important that, if a spark is turned into an actual parallel task, that task should be sufficiently expensive to justify the overhead of task creation, synchronisation, and so on.

A standard idiom for using @par@ is the following:
%
\begin{quote}
\begin{alltt}
let y = e2 in y `par` (e1 + y)
\end{alltt} 
\end{quote}
%
Assuming that @(+)@ starts by evaluating its first argument, @e2@ may be evaluated in parallel with @e1@, but the addition can only be performed after both @e1@ and @e2@ have been fully evaluated.

The dependence on the evaluation order of the arguments of @(+)@ leads to fragile code.  We can avoid that by using a second annotation, @pseq@, that has the same type as @par@, but evaluates it's first argument to WHNF before evaluating and returning the second:
%
\begin{quote}
\begin{alltt}
let x = e1
    y = e2
in x `par` (y `pseq` x + y)
\end{alltt} 
\end{quote}

% subsection control_parallelism (end)

\subsection{Splitting in parallel} % (fold)
\label{sub:splitting_in_parallel}

The simplest form of using @par@ to parallelise @hsplit@ is to create a spark for the right hand-side argument of @(++)@.  We achieve this as follows:
%
\begin{quote}
\begin{alltt}
hsplit :: [Point] -> Line -> [Point]
hsplit points line@(Line p1 p2)
  | length packed < 2 = p1:packed
  | otherwise         
  = let left  = hsplit packed (Line p1 pm)
        right = hsplit packed (Line pm p2)
    in
    right `par` (left ++ right)
  where
    \(\langle\)\textrm{as before}\(\rangle\)
\end{alltt}
\end{quote}
%
In addition, we apply the same sort of annotation to the two calls of @hsplit@ in the function @quickHull@.  That is important as the initial splitting of the dataset gives us the two tasks with the biggest granularity.  As we descent into the recursion, we eliminate points quickly from the input set by means of the list comprehension bound to @packed@.  In fact, for typical inputs, the set of points shrinks rather quickly.

It turns out that we get good speedup for two cores with this approach, but in our experiments, we weren't able to utilise any more cores.  Why is that?  Quickhull processes very long lists; we benchmarked data sets with 1 million to 8 million points.  Lists, however, are a sequential data structure and although annotations have been shown to be useful in parallelising list-based algorithms~\cite{?}, the approach appears to be limited when the input lists grow very large.  \TODO{Also mention that matters are even worse without -H500M or similar.}

In fact, much of the speed up on two cores is due to parallel garbage collection.  We also experimented with different uses of @par@ and @pseq@, but the only configuration that showed any speed up at all was the one described above.

% subsection splitting_in_parallel (end)

\subsection{Lists or trees} % (fold)
\label{sub:lists_or_trees}

Given the problems inherent in using large lists for data storage, we might want to try using a tree structure instead.  In fact, the use of a tree structure instead of a linear structure to expose more of the parallelism of an algorithm is a standard technique in parallel programming.  Specifically, we use finger trees~\cite{hinze-paterson:finger-trees} as exported by the module @Data.Sequence@ defined in the package @containers@~\cite{containers}.  Finger trees have been shown to be an efficient tree structure to represent sequences.

By replacing list operations by the corresponding operations on @Data.Sequence.Seq@ and replacing the list comprehensions by the explicit use of @filter@ and @map@, rewriting the parallel @hsplit@ to uses sequences instead of lists is straight forward.
%
\begin{quote}
\begin{alltt}
hsplit :: Seq Point -> Line -> Seq Point
hsplit points line@(Line p1 p2)
  | length packed < 2 = p1<|packed
  | otherwise             
  = let left  = hsplit packed (Line p1 pm)
        right = hsplit packed (Line pm p2)
    in
    right `par`
    (left >< right)
  where
    cross  = fmap (\string\p -> (distance p line, p)) 
                  points
    packed = fmap fst $ 
               filter (\string\(p, (c, _)) -> c > 0.0) 
                      (zip points cross)
    pm     = upper (>) cross
\end{alltt}
\end{quote}
%
Unfortunately, the use of sequences leads to a large performance penalty when compared to lists.  On 1 million points, the list version is about \TODO{?10?} times faster than the version based on sequences.  This is due to the very efficient implementation of list combinator by GHC.  In particular, it uses shortcut fusion to eliminate intermediate lists~\cite{gill-etal:build-fold}.  \TODO{Double check explanation by looking at the generated Core code.}  Consequently, tree-based sequences do not seem to provide a solution to the problem of data strorage for divide-and-conquer algorithms for data-intensive problems.

% subsection lists_or_trees (end)

% section parallelising_quickhull (end)


\section{The baseline} % (fold)
\label{sec:an_imperative_version}

The detailed benchmarks in Section~\ref{sec:benchmarks} include an optimised implementation of Quickhull in C to establish a baseline for the absolute performance of the sequential and parallel Haskell code.  However, as GHC's lack of low-level optimisations adds an inherent penalty to all Haskell code, we also establish a baseline of what is possible in Haskell by version of Quickhull implemented with the high-performance, mutable arrays of the @vector@ package~\cite{leshchinskiy:vector}.  By comparing the performance of the sequential C code with the sequential version of the Haskell code using @vector@, we deduce that the penalty due to the code generator is about \TODO{?} on x86 and about \TODO{?} on UltraSPARC T2.

After discussing a sequential version of the imperative algorithm, we will have a look at its parallelisation.  The imperative nature of the algorithm requires explicit concurrency control (thread spawning and synchronisation) to avoid races.

\subsection{Package vector} % (fold)
\label{sub:package_vector}

\TODO{Explain the imperative vector version of Quickhull (relative to the list version).}

% subsection package_vector (end)

\subsection{Explicit multi-threading} % (fold)
\label{sub:explicit_multi_threading}

\TODO{Explain the version of Quickhull based on @forkIO@ and @MVar@s.  Explain how this affords more control, but is also a more low-level approach to parallel programming by managing concurrency explicitly.  (For example, we can have new types of bugs, such as deadlocks and races.)  We also sacrifice purely functional programming for a state-based, imperative model.  The question is whether we can get better efficiency that way.  If so, how much better?}

% subsection explicit_multi_threading (end)

% section an_imperative_version (end)


\section{Nested data parallelism} % (fold)
\label{sec:nested_data_parallelism}

We can draw the following conclusions from our experience with the list, tree, and vector-based versions of Quickhull:
%
\begin{itemize}
  \item Sparse, boxed representations, such as trees, by default incur a high overhead on data-intensive problems.  Such representations stress the storage manager and are not cache friendly.
  \item Deforestation combined with GHC's unboxing transformation removes much of the overhead for lists.  However, very long lists are expensive during garbage collection and hinder parallel garbage collection.
  \item Dense, unboxed representations ---the classic data structures used in data-intensive applications--- lead to very efficient Haskell code, and can be parallelised.
\end{itemize}
%
Unfortunately, the latter option is tedious, especially with the very imperative programming style of Section~\ref{sec:an_imperative_version}.

% =======

briefly describes our extension of Haskell with \emph{nested data parallelism,} also known as Data Parallel Haskell, and how it enables the data-parallel implementation of divide-and-conquer algorithms, such as Quickhull, following the approach introduced by \citet{blelloch:ppa} with the language NESL.

\TODO{A lesson that we can learn from the list and tree implementation (and also the successful vector implementation) is that data-intensive code (a) needs a linear storage structure (trees take too much space and are not cache friendly) and (b) huge lists get the GC into trouble.  That is, we need to use large unboxed arrays.  We can manage them manually (e.g., with vector) using an imperative programming style or rely on sophisticated compiler-driven program transformations as in DPH.}

\TODO{Explain the DPH version of the algorithm.  Use this to introduce DPH as well.}

\TODO{At least in our experience so far, a practical advantage of DPH over @par@ for divide-and-conquer algorithms is that it is fairly clear how to parallelise them.  Using @par@ and @pseq@, there are usually a few different approaches, which can vary wildly in performance.  And the choice can be counterintuitive --- e.g., how using @nf@ in Quickhull kills scalability.}

% section nested_data_parallelism (end)


\section{Delaunay triangulation} % (fold)
\label{sec:delaunay_triangulation}

% section delaunay_triangulation (end)


\section{A divide-and-conquer skeleton} % (fold)
\label{sec:a_divide_and_conquer_skeleton}

\TODO{Play with a d-n-c skeleton and see whether it makes sense.}

% section a_divide_and_conquer_skeleton (end)


\section{Implementation challenges} % (fold)
\label{sec:implementation_challenges}

\TODO{Outline recent lessons learnt in implementing DPH, especially wrt to d-n-c algorithms: prealloc, ?what else?}

% section implementation_challenges (end)


\section{Benchmarks} % (fold)
\label{sec:benchmarks}

\TODO{Do we want to try to compare to qhull.org (need to convert the data first, to compare the same point sets.)}

% section benchmarks (end)


\section{Related work} % (fold)
\label{sec:related_work}

\TODO{Have a look at classic parallel implementations of divide-and-conquer algorithms.}
%
\begin{verbatim}
Shared memory:
* A parallel algorithm for computing the singular 
  value decomposition of a matrix (1994)
  by E R Jessup, D C Sorensen
  SIAM J. Matrix Anal. Appl
* A fully parallel algorithm for the symmetric 
  eigenvalue problem (1987)
  by J J DONGARRA, D C SORENSEN
  Tel-Aviv University
Distributed memory:
* Parallelizing the Divide and Conquer Algorithm 
  for the Symmetric Tridiagonal Eigenvalue Problem 
  on Distributed Memory Architectures (1998)
  by Françoise Tisseur ,  Jack Dongarra
  SIAM J. Sci. Comput
  http://citeseerx.ist.psu.edu/viewdoc/summary?doi=10.1.1.57.606
\end{verbatim}

% section related_work (end)


\paragraph{Acknowledgements.}
%
This research was funded in part by the Australian Research Council under grant number LP0989507.

\bibliographystyle{abbrvnat}
\bibliography{divide-conquer}

\end{document}
