dph-examples (empty) → 0.5.1.1

raw patch · 46 files changed

+3031/−0 lines, 46 filesdep +basedep +dph-basedep +dph-parsetup-changed

Dependencies added: base, dph-base, dph-par, dph-prim-par, dph-prim-seq, dph-seq, gloss, old-time, parseargs, random, vector

Files

• LICENSE +31/−0
• Setup.hs +2/−0
• dph-examples.cabal +240/−0
• imaginary/DotProduct/dph/DotProductVector.hs +7/−0
• imaginary/DotProduct/dph/DotProductVectorised.hs +15/−0
• imaginary/DotProduct/dph/Main.hs +45/−0
• imaginary/Evens/dph/EvensVector.hs +8/−0
• imaginary/Evens/dph/EvensVectorised.hs +16/−0
• imaginary/Evens/dph/Main.hs +45/−0
• imaginary/Reverse/dph/Main.hs +58/−0
• imaginary/Reverse/dph/ReverseVector.hs +21/−0
• imaginary/Reverse/dph/ReverseVectorised.hs +29/−0
• imaginary/SumSquares/dph/Main.hs +34/−0
• imaginary/SumSquares/dph/SumSquaresVector.hs +11/−0
• imaginary/SumSquares/dph/SumSquaresVectorised.hs +15/−0
• lib/Points2D/Generate.hs +98/−0
• lib/Points2D/Types.hs +49/−0
• lib/Randomish.hs +82/−0
• lib/Timing.hs +70/−0
• lib/Util.hs +25/−0
• real/NBody/Batch/Config.hs +57/−0
• real/NBody/Batch/MainArgs.hs +87/−0
• real/NBody/Common/Body.hs +103/−0
• real/NBody/Common/Dump.hs +32/−0
• real/NBody/Common/Util.hs +18/−0
• real/NBody/Common/World.hs +48/−0
• real/NBody/Gloss/Config.hs +74/−0
• real/NBody/Gloss/Draw.hs +104/−0
• real/NBody/Gloss/Main.hs +151/−0
• real/NBody/Gloss/MainArgs.hs +113/−0
• real/NBody/Solver.hs +59/−0
• real/NBody/Solver/ListBH/Solver.hs +161/−0
• real/NBody/Solver/NestedBH/Solver.hs +154/−0
• real/NBody/Solver/VectorBH/Solver.hs +163/−0
• real/NBody/Solver/VectorNaive/Solver.hs +18/−0
• spectral/QuickHull/dph/Main.hs +59/−0
• spectral/QuickHull/dph/QuickHullVect.hs +42/−0
• spectral/QuickHull/lib/SVG.hs +34/−0
• spectral/QuickHull/vector/Main.hs +71/−0
• spectral/QuickHull/vector/QuickHullIO.hs +269/−0
• spectral/QuickHull/vector/QuickHullSplit.hs +38/−0
• spectral/QuickHull/vector/QuickHullVector.hs +80/−0
• spectral/QuickSort/dph/Main.hs +50/−0
• spectral/QuickSort/dph/QuickSortVect.hs +25/−0
• spectral/SMVM/dph/Main.hs +102/−0
• spectral/SMVM/dph/SMVMVectorised.hs +18/−0

@@ -0,0 +1,31 @@+Copyright (c) 2001-2011, The DPH Team+All rights reserved.++The DPH Team is:+  Manuel M T Chakravarty+  Gabriele Keller+  Roman Leshchinskiy+  Ben Lippmeier+  George Roldugin++Redistribution and use in source and binary forms, with or without+modification, are permitted provided that the following conditions are met:+    * Redistributions of source code must retain the above copyright+      notice, this list of conditions and the following disclaimer.+    * Redistributions in binary form must reproduce the above copyright+      notice, this list of conditions and the following disclaimer in the+      documentation and/or other materials provided with the distribution.+    * Neither the name of the University of New South Wales nor the+      names of its contributors may be used to endorse or promote products+      derived from this software without specific prior written permission.++THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS ''AS IS'' AND ANY+EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE+DISCLAIMED. IN NO EVENT SHALL COPYRIGHT HOLDERS BE LIABLE FOR ANY+DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.

@@ -0,0 +1,2 @@+import Distribution.Simple+main = defaultMain

@@ -0,0 +1,240 @@++++Name:                dph-examples+Version:             0.5.1.1+License:             BSD3+License-file:        LICENSE+Author:              The DPH Team+Maintainer:          Ben Lippmeier <benl@ouroborus.net>+Build-Type:          Simple+Cabal-Version:       >=1.8+Stability:           experimental+Category:            Data Structures+Homepage:            http://www.cse.unsw.edu.au/~chak/project/dph/+Description:         Examples using the DPH library.+Synopsis:            Examples using the DPH library.++Flag gloss+ Description:        Enable graphical front ends that use the gloss library.+ Default:            False++-- With these examples, +--   the plain dph-NAME versions are parallel versions built against dph-par+--   while the dph_NAME-seq versions are sequential ones built against dph-seq++-- Sum of Squares -------------------------------------------------------------+Executable dph-sumsq+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  SumSquaresVector+                  SumSquaresVectorised+                  Timing Randomish+  hs-source-dirs: lib imaginary/SumSquares/dph+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3++Executable dph-sumsq-seq+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  SumSquaresVector+                  SumSquaresVectorised+                  Timing Randomish+  hs-source-dirs: lib imaginary/SumSquares/dph+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+++-- Dot Product ----------------------------------------------------------------+Executable dph-dotp+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  DotProductVector+                  DotProductVectorised+                  Timing Randomish+  hs-source-dirs: lib imaginary/DotProduct/dph+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3++Executable dph-dotp-seq+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  DotProductVector+                  DotProductVectorised+                  Timing Randomish+  hs-source-dirs: lib imaginary/DotProduct/dph+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+++-- Evens ----------------------------------------------------------------------+Executable dph-evens+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  EvensVector+                  EvensVectorised+                  Timing Randomish+  hs-source-dirs: imaginary/Evens/dph lib+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3++Executable dph-evens-seq+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  EvensVector+                  EvensVectorised+                  Timing Randomish+  hs-source-dirs: imaginary/Evens/dph lib+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+++-- Primes ---------------------------------------------------------------------+-- Primes dies with a bounds check error+-- Executable dph-primes+--  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*+--  Main-is:        Main.hs+--  hs-source-dirs: imaginary/Primes lib+--  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3++-- Executable dph-primes-seq+--  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*+--  Main-is:        Main.hs+--  hs-source-dirs: imaginary/Primes lib+--  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+++-- Reverse --------------------------------------------------------------------+Executable dph-reverse+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  ReverseVector+                  ReverseVectorised+                  Util Timing Randomish+  hs-source-dirs: imaginary/Reverse/dph lib+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+++Executable dph-reverse-seq+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  ReverseVector+                  ReverseVectorised+                  Util Timing Randomish+  hs-source-dirs: imaginary/Reverse/dph lib+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+++-- Words ---------------------------------------------------------------------+-- This builds and runs, but takes ~5min. +-- Executable dph-words+--  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*+--  Main-is:        Main.hs+--  hs-source-dirs: imaginary/Words lib+--  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3 -fno-liberate-case++-- This one also stalls in SpecConstr, like quicksort-seq+-- Executable dph-words-seq+--   Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*+--   Main-is:        Main.hs+--   hs-source-dirs: imaginary/Words lib+--   ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+++-- SMVM ---------------------------------------------------------------------+Executable dph-smvm+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  SMVMVectorised+                  Timing+  hs-source-dirs: spectral/SMVM/dph lib+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3++Executable dph-smvm-seq+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  SMVMVectorised+                  Timing+  hs-source-dirs: spectral/SMVM/dph lib+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3++++-- QuickHull ------------------------------------------------------------------+Executable dph-quickhull+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  QuickHullVect+                  Timing Points2D.Types SVG+  hs-source-dirs: lib spectral/QuickHull/dph spectral/QuickHull/lib+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3++Executable dph-quickhull-seq+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  QuickHullVect+                  Timing Points2D.Types SVG+  hs-source-dirs: lib spectral/QuickHull/dph spectral/QuickHull/lib+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3++Executable dph-quickhull-vector+  Build-depends: +        base          == 4.*,+        vector        == 0.7.*,    +        old-time      == 1.0.*,+        random        == 1.0.*+  Main-is:        Main.hs+  other-modules:  QuickHullIO+                  QuickHullSplit+                  QuickHullVector+                  Timing Points2D.Types SVG+  hs-source-dirs: lib spectral/QuickHull/vector spectral/QuickHull/lib+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+++-- QuickSort ------------------------------------------------------------------+Executable dph-quicksort+  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*+  Main-is:        Main.hs+  other-modules:  QuickSortVect+                  Timing+  hs-source-dirs: lib spectral/QuickSort/dph+  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3++-- This one stalls in SpecConstr then eventually runs out of memory.+-- Executable dph-quicksort-seq+--  Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*+--  Main-is:        Main.hs+--  other-modules:  QuickSortVect+--                  Timing+--  hs-source-dirs: lib spectral/QuickSort/dph+--  ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+++-- NBody ----------------------------------------------------------------------+Executable dph-nbody+  Main-is:        Main.hs+  other-modules:  Common.Dump Common.World Common.Body Common.Util +                  Solver Solver.ListBH.Solver+                         Solver.NestedBH.Solver+                         Solver.VectorBH.Solver+                         Solver.VectorNaive.Solver+                  Timing Points2D.Types Points2D.Generate++  if flag(gloss)+    other-modules:  Gloss.MainArgs Gloss.Draw Gloss.Config+    Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*, parseargs == 0.1.*, gloss == 1.3.*+    hs-source-dirs: lib real/NBody real/NBody/Gloss+    ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3 -fno-liberate-case+  else+    other-modules:  Batch.MainArgs Batch.Config+    Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-par == 0.5.*, dph-par == 0.5.*, random == 1.0.*, old-time == 1.0.*, parseargs == 0.1.*+    hs-source-dirs: lib real/NBody real/NBody/Batch+    ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-par -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3 -fno-liberate-case+  endif++-- Executable dph-nbody-seq+--  Main-is:        Main.hs+--  if flag(gloss)+--    Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*, parseargs == 0.1.*, gloss == 1.3.*+--    hs-source-dirs: lib real/NBody real/NBody/Gloss+--    ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+--  else+--    Build-depends:  base == 4.4.*, vector == 0.7.*, dph-base == 0.5.*, dph-prim-seq == 0.5.*, dph-seq == 0.5.*, random == 1.0.*, old-time == 1.0.*, parseargs == 0.1.*+--    hs-source-dirs: lib real/NBody real/NBody/Batch+--    ghc-options:    -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3+--  endif

@@ -0,0 +1,7 @@++module DotProductVector (dotV) where+import Data.Vector.Unboxed		(Vector)+import qualified Data.Vector.Unboxed	as V++dotV :: Vector Double -> Vector Double -> Double+dotV vec1 vec2 = V.sum $ V.zipWith (*) vec1 vec2

@@ -0,0 +1,15 @@+{-# LANGUAGE ParallelArrays #-}+{-# OPTIONS -fvectorise #-}+module DotProductVectorised ( dotPA ) where++import Data.Array.Parallel+import Data.Array.Parallel.Prelude.Double as D++import qualified Prelude++dotPA :: PArray Double -> PArray Double -> Double+{-# NOINLINE dotPA #-}+dotPA v w = dotp' (fromPArrayP v) (fromPArrayP w)++dotp' :: [:Double:] -> [:Double:] -> Double+dotp' v w = D.sumP (zipWithP (*) v w)

@@ -0,0 +1,45 @@++import Timing+import Randomish+import System.Environment+import Data.Array.Parallel	as P+import Data.Array.Parallel.PArray	as P+import qualified Data.Vector.Unboxed	as V+import qualified DotProductVector	as V+import qualified DotProductVectorised	as Z++main :: IO ()+main + = do	args	<- getArgs+	case args of+	  [alg, len] 	-> run alg (read len) +	  _		-> usage++usage+ = putStr $ unlines+ 	[ "usage: dotp <alg> <length>"+ 	, "  alg one of " ++ show ["vectorised", "vector"] ]+	+run alg len+ = do	let vec1 = randomishDoubles len 0 1 1234+	let vec2 = randomishDoubles len 0 1 12345++	(result, tElapsed) <- runAlg alg vec1 vec2++	putStr	$ prettyTime tElapsed+	putStr	$ (take 12 $ show result) ++ "\n"++runAlg "vectorised" vec1 vec2 + = do	let arr1 = fromUArrPA' vec1+	let arr2 = fromUArrPA' vec2+	arr1 `seq` arr2 `seq` return ()++	time	$ let result	= Z.dotPA arr1 arr2+		  in  result `seq` return result++runAlg "vector" vec1 vec2+ = do	vec1 `seq` vec2 `seq` return ()++	time	$ let result	= V.dotV vec1 vec2 +		  in  result `seq` return result+		

@@ -0,0 +1,8 @@++module EvensVector (evensV) where+import Data.Vector.Unboxed		(Vector)+import qualified Data.Vector.Unboxed	as V++evensV :: Vector Int -> Vector Int+evensV ints = V.filter (\x -> x `mod` 2 == 0) ints+

@@ -0,0 +1,16 @@+{-# LANGUAGE ParallelArrays #-}+{-# OPTIONS -fvectorise #-}++module EvensVectorised (evensPA) where+import Data.Array.Parallel+import Data.Array.Parallel.Prelude.Int+import qualified Prelude as P++-- | DPH filter opereations are reasonably involved because they use packByTag.+evens :: [:Int:] -> [:Int:]+evens ints = filterP (\x -> x `mod` 2 == 0) ints++evensPA :: PArray Int -> PArray Int+{-# NOINLINE evensPA #-}+evensPA arr = toPArrayP (evens (fromPArrayP arr))+

@@ -0,0 +1,45 @@++import EvensVectorised			as Z+import EvensVector			as V+import qualified Data.Vector.Unboxed	as V+import Timing+import System.Environment+import Data.Array.Parallel.PArray	as P++main :: IO ()+main + = do	args	<- getArgs+	case args of+	  [alg, len]	-> run alg (read len)+	  _ 		-> usage++usage+ = putStr $ unlines+ 	[ "usage: evens <alg> <length>"+ 	, "  alg one of " ++ show ["vectorised", "vector"] ]++run alg len+ = do	let vec	= V.fromList [0..len - 1]+	runAlg alg vec+	+runAlg "vectorised" vec+ = do	let arr	= P.fromUArrPA' vec+	arr `seq` return ()++	(arr', tElapsed)+	 <- time $ let	arr' = evensPA arr+		   in	arr' `seq` return arr'+					+	putStr   $ prettyTime tElapsed+	putStrLn $ show $ P.length arr'++runAlg "vector" vec+ = do	vec `seq` return ()+	+	(vec', tElapsed)+	 <- time $ let	vec' = evensV vec+		   in	vec' `seq` return vec'+		+	putStr   $ prettyTime tElapsed+	putStrLn $ show $ V.length vec'+	

@@ -0,0 +1,58 @@++import Util+import Timing+import Randomish+import System.Environment+import System.Random+import Control.Exception+import qualified ReverseVector                  as RV+import qualified ReverseVectorised              as RD+import qualified Data.Array.Parallel.PArray     as P+import qualified Data.Vector.Unboxed            as V+++main + = do   args    <- getArgs+        +        case args of+         [alg, count]   +           |  not $ isPowerOfTwo (read count)+           -> error "reverse: length of array must be a power of two."++           | otherwise -> run alg (read count)++         _ -> usage+++run "vectorised" count+ = do   let arr = P.fromList [0 .. count - 1]+        arr `seq` return ()     +                +        (arrReversed, tElapsed)+         <- time+         $  let  arr'    = RD.treeReversePA arr+            in   arr' `seq` return arr'++        print arrReversed+        putStr  $ prettyTime tElapsed++run "vector" count+ = do   let arr = V.fromList [0 .. count - 1]+        arr `seq` return ()+        +        (arrReversed, tElapsed)+         <- time+         $  let arr'    = RV.treeReverse arr+            in  arr' `seq` return arr'+            +        print arrReversed+        putStr  $ prettyTime tElapsed++run _ _+ = usage+++usage   = putStr $ unlines+        [ "usage: reverse <algorithm> <count>\n"+        , "  algorithm one of " ++ show ["vectorised"]+        , ""]

@@ -0,0 +1,21 @@++module ReverseVector+        (treeReverse)+where+import qualified Data.Vector.Unboxed as V+import Data.Vector.Unboxed              (Vector)++-- | Reverse the elements in an array using a tree.+treeReverse :: Vector Int -> Vector Int+{-# NOINLINE treeReverse #-}+treeReverse xx+        | V.length xx == 1+        = xx+        +        | otherwise+        = let   len     = V.length xx+                half    = len `div` 2+                s1      = V.slice 0    half  xx+                s2      = V.slice half half  xx         +          in    treeReverse s2 V.++ treeReverse s1+

@@ -0,0 +1,29 @@+{-# LANGUAGE ParallelArrays #-}+{-# OPTIONS -fvectorise #-}+module ReverseVectorised        +        (treeReversePA)+where+import Data.Array.Parallel+import Data.Array.Parallel.Prelude.Int+import qualified Prelude as P+++treeReversePA :: PArray Int -> PArray Int+{-# NOINLINE treeReversePA #-}+treeReversePA ps+        = toPArrayP (treeReverse (fromPArrayP ps))+++-- | Reverse the elements in an array using a tree.+treeReverse :: [:Int:] -> [:Int:]+{-# NOINLINE treeReverse #-}+treeReverse xx+        | lengthP xx == 1+        = xx+        +        | otherwise+        = let   len     = lengthP xx+                half    = len `div` 2+                s1      = sliceP 0    half xx+                s2      = sliceP half half  xx           +          in    concatP (mapP treeReverse [: s2, s1 :])

@@ -0,0 +1,34 @@++import Timing+import Randomish+import System.Environment+import qualified Data.Vector.Unboxed	as V+import qualified SumSquaresVector	as V+import qualified SumSquaresVectorised	as Z++main :: IO ()+main + = do	args	<- getArgs+	case args of+	  [alg, len] 	-> run alg (read len) +	  _		-> usage++usage+ = putStr $ unlines+ 	[ "usage: sumsq <alg> <length>"+ 	, "  alg one of " ++ show ["vectorised", "vector"] ]+	+run alg num+ = do	(result, tElapsed) <- runAlg alg num++	putStr	$ prettyTime tElapsed+	print result++runAlg "vectorised" num+ =	time	$ let result	= Z.sumSq num+		  in  result `seq` return result++runAlg "vector" num+ =	time	$ let result	= V.sumSq num+		  in  result `seq` return result+		

@@ -0,0 +1,11 @@++module SumSquaresVector where+import qualified Data.Vector.Unboxed	 as V+	+sumSq :: Int -> Double+sumSq num	+	= V.sum +	$ V.map (\x -> x * x) +	$ V.map fromIntegral +	$ V.enumFromTo 1 num+	

@@ -0,0 +1,15 @@+{-# LANGUAGE ParallelArrays #-}+{-# OPTIONS -fvectorise #-}+module SumSquaresVectorised (sumSq) where+import Data.Array.Parallel+import Data.Array.Parallel.Prelude.Int 		as I+import Data.Array.Parallel.Prelude.Double	as D+import qualified Prelude++sumSq :: Int -> Double+{-# NOINLINE sumSq #-}+sumSq n	+ = 	D.sumP +	(mapP (\x -> x D.* x) +	(mapP D.fromInt +	(enumFromToP 1 n)))

@@ -0,0 +1,98 @@+{-# LANGUAGE TypeOperators #-}+module Points2D.Generate +	( genPointsUniform+	, genPointsUniformWithSeed+	, genPointsDisc+	, genPointsCombo+	, pointsPArrayOfUArray )+where+import Points2D.Types+import Randomish+import qualified Data.Array.Parallel.Unlifted 	    as U+import qualified Data.Array.Parallel 	    as P+import qualified Data.Array.Parallel.PArray         as P+import Data.Array.Parallel.PArray		    (PArray)+import Control.Exception++-- Random points generation+-- IMPORTANT: We use the same seed with the same random generator in all+--            quickhull codes.  The asymptotic work complexity of quickhull+--            is between O (N) and O (N^2) depending on the input.+--            To compare benchmark results, they always need to use the same+--            input.+seed :: Int+seed 	= 42742++-- | Some uniformly distributed points+genPointsUniform +	:: Int			-- ^ number of points+	-> Double		-- ^ minimum coordinate+	-> Double		-- ^ maximum coordinate+        -> U.Array (Double, Double)++genPointsUniform n pointMin pointMax+ = let  pts             = randomishDoubles (n*2) pointMin pointMax seed+        xs              = U.extract pts 0 n+        ys              = U.extract pts n n+   in   U.zip xs ys+++-- | Some uniformly distributed points+genPointsUniformWithSeed+	:: Int			-- ^ seed+	-> Int			-- ^ number of points+	-> Double		-- ^ minimum coordinate+	-> Double		-- ^ maximum coordinate+        -> U.Array (Double, Double)++genPointsUniformWithSeed seed' n pointMin pointMax+ = let  pts             = randomishDoubles (n*2) pointMin pointMax seed'+        xs              = U.extract pts 0 n+        ys              = U.extract pts n n+   in   U.zip xs ys+++-- | Some points distributed as a disc+genPointsDisc+	:: Int			-- ^ number of points+	-> (Double, Double) 	-- ^ center of disc+	-> Double 		-- ^ radius of disc+        -> U.Array (Double, Double)++genPointsDisc n (originX, originY) radiusMax+ = let	radius = randomishDoubles n 0     radiusMax seed+        angle  = randomishDoubles n (-pi) pi        (seed + 1234)++	makeXY r a	+	 	= ( originX + r * cos a+		  , originY + r * sin a)	++    in	originX `seq` originY `seq` U.zipWith makeXY radius angle+++-- | A point cloud with areas of high an low density+genPointsCombo +	:: Int 			-- ^ number of points+        -> U.Array (Double, Double)++genPointsCombo n+ 	=  genPointsDisc    (n `div` 5) (250, 250) 200+	U.+:+ genPointsDisc (n `div` 5) (100, 100) 80 +	U.+:+ genPointsDisc (n `div` 5) (150, 300) 30 +	U.+:+ genPointsDisc (n `div` 5) (500, 120) 30 +	U.+:+ genPointsDisc (n `div` 5) (300, 200) 150+++-- | Convert a list of points to a PArray+pointsPArrayOfUArray +	:: U.Array (Double, Double)+	-> IO (PArray Point)++pointsPArrayOfUArray ps+  = do+      let pts = makePointsPA +			(P.fromUArrPA' (U.fsts ps))+ 			(P.fromUArrPA' (U.snds ps))+      evaluate $ P.nf pts+      return pts+

@@ -0,0 +1,49 @@+{-# LANGUAGE ParallelArrays #-}+{-# OPTIONS -fvectorise #-}++module Points2D.Types +	( Point+	, Line+	, makePoints, makePointsPA+	, xsOf, xsOfPA+	, ysOf, ysOfPA)+where+import Data.Array.Parallel++type Point = (Double, Double)+type Line  = (Point, Point)++-- | Make some points from their components.+makePoints :: [:Double:] -> [:Double:] -> [:Point:]+makePoints = zipP+++-- | Make some points from their components, as a `PArray`.+makePointsPA :: PArray Double -> PArray Double -> PArray Point+{-# NOINLINE makePointsPA #-}+makePointsPA xs ys+	= toPArrayP (makePoints (fromPArrayP xs) (fromPArrayP ys))+++-- | Take the x values of some points.+xsOf :: [:Point:] -> [:Double:]+xsOf ps = [: x | (x, _) <- ps :]+++-- | Take the x values of some points as a `PArray`.+xsOfPA :: PArray Point -> PArray Double+{-# NOINLINE xsOfPA #-}+xsOfPA ps = toPArrayP (xsOf (fromPArrayP ps))+++-- | Take the y values of some points.+ysOf :: [:Point:] -> [:Double:]+ysOf ps = [: y | (_, y) <- ps :]+++-- | Take the y values of some points as a `PArray`.+ysOfPA :: PArray Point -> PArray Double+{-# NOINLINE ysOfPA #-}+ysOfPA ps = toPArrayP (ysOf (fromPArrayP ps))++

@@ -0,0 +1,82 @@+{-# LANGUAGE BangPatterns #-}++module Randomish+ 	( randomishInts+	, randomishDoubles)+where+import Data.Word+import Data.Vector.Unboxed			(Vector)+import qualified Data.Vector.Unboxed.Mutable	as MV+import qualified Data.Vector.Unboxed		as V+import qualified Data.Vector.Generic		as G		+++-- | Use the "minimal standard" Lehmer generator to quickly generate some random+--   numbers with reasonable statistical properties. By "reasonable" we mean good+--   enough for games and test data, but not cryptography or anything where the+--   quality of the randomness really matters.+-- +--   From "Random Number Generators: Good ones are hard to find"+--   Stephen K. Park and Keith W. Miller.+--   Communications of the ACM, Oct 1988, Volume 31, Number 10.+--+randomishInts +	:: Int 			-- Length of vector.+	-> Int 			-- Minumum value in output.+	-> Int 			-- Maximum value in output.+	-> Int 			-- Random seed.	+	-> Vector Int		-- Vector of random numbers.++randomishInts !len !valMin' !valMax' !seed'+	+ = let	-- a magic number (don't change it)+	multiplier :: Word64+	multiplier = 16807++	-- a merzenne prime (don't change it)+	modulus	:: Word64+	modulus	= 2^(31 :: Integer) - 1++	-- if the seed is 0 all the numbers in the sequence are also 0.+	seed	+	 | seed' == 0	= 1+	 | otherwise	= seed'++	!valMin	= fromIntegral valMin'+	!valMax	= fromIntegral valMax' + 1+	!range	= valMax - valMin++	{-# INLINE f #-}+	f x		= multiplier * x `mod` modulus+ in G.create +     $ do	+	vec		<- MV.new len++	let go !ix !x +	  	| ix == len	= return ()+		| otherwise+		= do	let x'	= f x+			MV.write vec ix $ fromIntegral $ (x `mod` range) + valMin+			go (ix + 1) x'++	go 0 (f $ f $ f $ fromIntegral seed)+	return vec+++-- | Generate some randomish doubles with terrible statistical properties.+--   This is good enough for test data, but not much else.+randomishDoubles +	:: Int			-- Length of vector+	-> Double		-- Minimum value in output+	-> Double		-- Maximum value in output+	-> Int			-- Random seed.+	-> Vector Double	-- Vector of randomish doubles.++randomishDoubles !len !valMin !valMax !seed+ = let	range	= valMax - valMin++	mx	= 2^(30 :: Integer) - 1+	mxf	= fromIntegral mx+	ints	= randomishInts len 0 mx seed+	+   in	V.map (\n -> valMin + (fromIntegral n / mxf) * range) ints

@@ -0,0 +1,70 @@+module Timing+	(time, showTime, prettyTime)+where+import GHC.Exts	(traceEvent)+import System.CPUTime+import System.Time+++-- Time -----------------------------------------------------------------------+data Time +	= Time +	{ cpu_time  :: Integer+        , wall_time :: Integer+        }++zipT :: (Integer -> Integer -> Integer) -> Time -> Time -> Time+zipT f (Time cpu1 wall1) (Time cpu2 wall2) +	= Time (f cpu1 cpu2) (f wall1 wall2)++minus :: Time -> Time -> Time+minus = zipT (-)+++-- TimeUnit -------------------------------------------------------------------+type TimeUnit +	= Integer -> Integer++milliseconds :: TimeUnit+milliseconds n = n `div` 1000000000++cpuTime :: TimeUnit -> Time -> Integer+cpuTime f = f . cpu_time++wallTime :: TimeUnit -> Time -> Integer+wallTime f = f . wall_time+++-- | Get the current time.+getTime :: IO Time+getTime =+  do+    cpu          <- getCPUTime+    TOD sec pico <- getClockTime+    return $ Time cpu (pico + sec * 1000000000000)+++-- | Show a time as a string, in milliseconds.+showTime :: Time -> String+showTime t = (show $ wallTime milliseconds t)+          ++ "/"+          ++ (show $ cpuTime  milliseconds t)++-- | Pretty print the times.+prettyTime :: Time -> String+prettyTime t+	= unlines+	[ "elapsedTimeMS   = " ++ (show $ wallTime milliseconds t)+	, "cpuTimeMS       = " ++ (show $ cpuTime  milliseconds t) ]+++-- Timing benchmarks ----------------------------------------------------------+time :: IO a -> IO (a, Time)+{-# NOINLINE time #-}+time p = do+           start <- getTime+           traceEvent "dph-examples: start timing"+           x     <- p+           traceEvent "dph-examples: finished timing"+           end   <- getTime+           return (x, end `minus` start)

@@ -0,0 +1,25 @@+{-# LANGUAGE PatternGuards #-}+module Util +	( isPowerOfTwo+	, isSorted)+where+	++-- | Check if an integer is a power of two.+isPowerOfTwo :: Int -> Bool+isPowerOfTwo n+	| 0	<- n		= True+	| 2	<- n		= True+	| n `mod` 2 == 0	= isPowerOfTwo (n `div` 2)+	| otherwise		= False+++isSorted :: Ord a => [a] -> Bool+isSorted xx+ = case xx of+	[]		-> True+	[x]		-> True+	(x1:x2:rest)+	 | x1 <= x2	-> isSorted (x2 : rest)+	 | otherwise	-> False+

@@ -0,0 +1,57 @@++module Batch.Config+	( Config(..)+	, loadConfig)+where+import Batch.MainArgs+import System.Console.ParseArgs+import Data.Maybe+++-- | Program config+data Config+	= Config {+	-- What solver to use+	  configSolverName	:: String++	-- System setup+	, configBodyCount	:: Int+	, configBodyMass	:: Double+	, configTimeStep	:: Double+	, configEpsilon		:: Double++	-- Initial conditions.+	, configStartDiscSize	:: Double+	, configStartSpeed	:: Double++	-- Terminating conditions.+	, configMaxSteps	:: Int+	+	-- Dump points to file+	, configDumpFinal	:: Maybe FilePath }+	++-- | Load program config from its command line arguments.	+loadConfig :: Args MainArg -> Config+loadConfig args+ = let	Just solverName	= getArgString	args ArgSolver+	Just timeStep	= getArgDouble	args ArgTimeStep+	Just bodyCount	= getArgInt	args ArgBodyCount+	Just bodyMass	= getArgDouble  args ArgBodyMass+	Just epsilon	= getArgDouble	args ArgEpsilon+	Just discSize	= getArgDouble	args ArgDiscSize+	Just startSpeed	= getArgDouble	args ArgStartSpeed++	Just maxSteps	= getArgInt	args ArgMaxSteps+	mFilePath	= getArgString	args ArgDumpFinal++   in	Config+	{ configSolverName	= solverName+	, configBodyCount	= bodyCount+	, configBodyMass	= bodyMass+	, configTimeStep	= timeStep+	, configEpsilon		= epsilon+	, configStartDiscSize	= discSize+	, configStartSpeed	= startSpeed+	, configMaxSteps	= maxSteps+	, configDumpFinal	= mFilePath }

@@ -0,0 +1,87 @@++module Batch.MainArgs+	( MainArg(..)+	, mainArgs)+where+import System.Console.ParseArgs++data MainArg+	= ArgHelp+	| ArgSolver+	| ArgMaxSteps+	| ArgTimeStep+	| ArgBodyCount+	| ArgBodyMass+	| ArgEpsilon+	| ArgDiscSize+	| ArgStartSpeed+		+	-- Dump output+	| ArgDumpFinal+	deriving (Eq, Ord, Show)+	+mainArgs :: [Arg MainArg]+mainArgs+ = 	[ Arg	{ argIndex	= ArgHelp+		, argAbbr	= Just 'h'+		, argName	= Just "help"+		, argData	= Nothing+		, argDesc	= "Print this usage help." }++	-- Solver selection.+	, Arg	{ argIndex	= ArgSolver+		, argAbbr	= Just 's'+		, argName	= Just "solver"+		, argData	= argDataDefaulted "name" ArgtypeString "vector-naive"+		, argDesc	= "One of: list-bh, vector-naive, vector-bh, nested-bh. (default vector-naive)" }++	-- Simulation setup.+	, Arg	{ argIndex	= ArgTimeStep+		, argAbbr	= Just 't'+		, argName	= Just "timestep"+		, argData	= argDataDefaulted "Double" ArgtypeDouble 1+		, argDesc	= "Time step between states (default 1)" }++	, Arg	{ argIndex	= ArgBodyCount+		, argAbbr	= Just 'b'+		, argName	= Just "bodies"+		, argData	= argDataDefaulted "Int" ArgtypeInt 200 +		, argDesc	= "Number of bodies in simulation (default 200)" }++	, Arg	{ argIndex	= ArgBodyMass+		, argAbbr	= Just 'm'+		, argName	= Just "mass"+		, argData	= argDataDefaulted "Double" ArgtypeDouble 10+		, argDesc	= "Mass of each body (default 10)" }++	, Arg	{ argIndex	= ArgEpsilon+		, argAbbr	= Just 'e'+		, argName	= Just "epsilon"+		, argData	= argDataDefaulted "Double" ArgtypeDouble 100+		, argDesc	= "Smoothing parameter (default 100)" }+		+	, Arg	{ argIndex	= ArgDiscSize+		, argAbbr	= Just 'd'+		, argName	= Just "disc"+		, argData	= argDataDefaulted "Double" ArgtypeDouble 50+		, argDesc	= "Starting size of disc containing bodies (default 50)" }++	, Arg	{ argIndex	= ArgStartSpeed+		, argAbbr	= Just 'p'+		, argName	= Just "speed"+		, argData	= argDataDefaulted "Double" ArgtypeDouble 0.5+		, argDesc	= "Starting rotation speed of bodies (default 0.5)" }++	, Arg	{ argIndex	= ArgMaxSteps+		, argAbbr	= Nothing+		, argName	= Just "max-steps"+		, argData	= argDataDefaulted "steps" ArgtypeInt 1000+		, argDesc	= "Exit simulation after this many steps (default 1000)" }++	, Arg	{ argIndex	= ArgDumpFinal+		, argAbbr	= Nothing+		, argName	= Just "dump-final"+		, argData	= argDataOptional "FilePath" ArgtypeString+		, argDesc	= "Dump final body positions and masses to file" }++	]

@@ -0,0 +1,103 @@+{-# LANGUAGE BangPatterns #-}++-- | Massful bodies in the simulation.+module Common.Body+	( Velocity+	, Accel+	, MassPoint+	, Body+	+	, accel++	, unitBody+	, massPointOfBody+	, setMassOfBody+	, setAccelOfBody+	, setStartVelOfBody+	, advanceBody)+where+import Common.Util+++-- Types ----------------------------------------------------------------------+-- We're using tuples instead of ADTs so we can put them in unboxed vectors.++-- | The velocity of a point.+type Velocity	= (Double, Double)++-- | The acceleration of a point.+type Accel	= (Double, Double)++-- | A point in 2D space with its mass.+type MassPoint	= (Double, Double, Double)++-- | Bodies consist of a MassPoint, but also carry their velocity+--   and acceleration between steps of the simulation.+type Body	= (MassPoint, Velocity, Accel)+++-- Acceleration ---------------------------------------------------------------+-- | Calculate the acceleration on a point due to some other point.+accel 	:: Double 	-- ^ Smoothing parameter.+	-> MassPoint	-- ^ The point being acclerated.+	-> MassPoint	-- ^ Neibouring point.+	-> Accel++{-# INLINE accel #-}+accel epsilon (!x1, !y1, _) (!x2, !y2, !m)  + = (aabs * dx / r , aabs * dy / r)  + where	dx   = x1 - x2 +	dy   = y1 - y2 +	rsqr = (dx * dx) + (dy * dy) + epsilon * epsilon+	r    = sqrt rsqr +	aabs = m / rsqr +++-- Body -----------------------------------------------------------------------+-- | Make a body with unit mass and zero vel and acc.+unitBody :: Double -> Double -> Body+unitBody x y+	= ((x, y, 1), (0, 0), (0, 0))+++-- | Take the MassPoint of a body.+massPointOfBody :: Body -> MassPoint+massPointOfBody (mp, vel, acc)	+	= mp+++-- | Set the mass of a body.+setMassOfBody :: Double -> Body -> Body+setMassOfBody mass ((x, y, _), vel, acc)+	= ((x, y, mass), vel, acc)+++-- | Set the acceleration of a body.+setAccelOfBody :: Accel -> Body -> Body+setAccelOfBody acc' (mp, vel, _)	+	= (mp, vel, acc')++	+-- | Set the starting velocity of a body.+--   It is set to rotate around the origin, with the speed proportional+--   to the sqrt of the distance from it. This seems to make nice simulations.+setStartVelOfBody :: Double -> Body -> Body+setStartVelOfBody startVel (mp@(x, y, mass), vel, acc)+ = let	pos		= (x, y)+	(x', y')	= normaliseV (x, y)+   	vel'		= (y', -x')+	vel''		= mulSV (sqrt (magV pos) * startVel) vel'+	+   in	(mp, vel'', acc)+++-- | Advance a body forwards in time.+advanceBody :: Double -> Body -> Body+advanceBody time +	( (px, py, mass) +	, (vx, vy) +	, acc@(ax, ay))++  =	( (px + time * vx, py + time * vy, mass)+	, (vx + time * ax, vy + time * ay)+	, acc)

@@ -0,0 +1,32 @@+{-# LANGUAGE ScopedTypeVariables #-}+module Common.Dump+	(dumpWorld)+where+import Common.Body+import Common.World+import System.IO+import Data.List+import qualified Data.Vector.Unboxed	as V+++-- | Dump the bodies in a world to a file.+dumpWorld :: World -> FilePath -> IO ()+dumpWorld world filePath+ = do	h	<- openFile filePath WriteMode+	mapM_ (hWriteBody h)+		$ V.toList +		$ worldBodies world++-- | Write a single body to a file.+hWriteBody :: Handle -> Body -> IO ()+hWriteBody h ((px, py, mass), (vx, vy), (ax, ay))+ 	= hPutStrLn h +	$ concat +	$ (  (padRc 8 ' ' $ show mass)+	  :  " " : (intersperse " " $ map (padRc 22 ' ' . show) [ px, py, vx, vy, ax, ay ]))+		++-- | Right justify a doc, padding with a given character.+padRc :: Int -> Char -> String -> String+padRc n c str+	= replicate (n - length str) c ++ str

@@ -0,0 +1,18 @@++module Common.Util+	( magV+	, mulSV+	, normaliseV)+where+	+-- | The magnitude of a vector.+magV :: (Double, Double) -> Double+magV (x, y)	= sqrt (x * x + y * y)+	+-- | Multiply a vector by a scalar.+mulSV :: Double -> (Double, Double) -> (Double, Double)+mulSV s (x, y)	= (s * x, s * y)+	+-- | Normalise a vector, so it has a magnitude of 1.+normaliseV :: (Double, Double) -> (Double, Double)+normaliseV v	= mulSV (1 / magV v) v

@@ -0,0 +1,48 @@+{-# LANGUAGE BangPatterns #-}++module Common.World+	( World(..)+	, advanceWorld)+where+import Common.Body+import qualified Data.Vector.Unboxed		as V++data World +	= World+	{ -- | Bodies in the simulation.+	  worldBodies	:: !(V.Vector Body)++	  -- | Number of steps taken in the simulation so far.+	, worldSteps	:: !Int }+++-- | Advance the world forward in time.+advanceWorld +	:: (V.Vector MassPoint	-> V.Vector Accel)+				-- ^ Fn to compute accelerations of each point.+	-> Double		-- ^ Time step.+	-> World+	-> World++advanceWorld calcAccels timeStep world+ = let	+	-- Calculate the accelerations on each body.+	accels	= calcAccels +		$ V.map massPointOfBody +		$ worldBodies world++	-- Apply the accelerations to the bodies and advance them.+	bodies'	= V.zipWith +		(\body (ax, ay) +			-> advanceBody timeStep+				(setAccelOfBody (-ax, -ay) body))+		(worldBodies world)+		accels++	-- Update the world.+	steps'	= worldSteps world + 1++   in	world	{ worldBodies	= bodies'+		, worldSteps	= steps' }++	

@@ -0,0 +1,74 @@++module Gloss.Config+	( Config (..)+	, loadConfig)+where+import Gloss.MainArgs+import System.Console.ParseArgs+import Data.Maybe++-- | Program config.+data Config+	= Config {+	 +	-- How to present the output.+	  configWindowSize	:: Maybe Int+	, configShouldDrawTree	:: Bool+	, configRate		:: Int++	-- What solver to use+	, configSolverName	:: String++	-- System setup+	, configBodyCount	:: Int+	, configBodyMass	:: Double+	, configTimeStep	:: Double+	, configEpsilon		:: Double++	-- Initial conditions.+	, configStartDiscSize	:: Double+	, configStartSpeed	:: Double++	-- Terminating conditions.+	, configMaxSteps	:: Maybe Int+	+	-- dump points to file+	, configDumpFinal	:: Maybe FilePath }+	++-- | Load program config from its command line arguments.	+loadConfig :: Args MainArg -> Config+loadConfig args+ = let	mWindowSize	= getArgInt	args ArgGloss+	Just solverName	= getArgString	args ArgSolver+	shouldDrawTree	= gotArg  	args ArgDrawTree+	Just timeStep	= getArgDouble	args ArgTimeStep+	Just rate	= getArgInt	args ArgRate+	Just bodyCount	= getArgInt	args ArgBodyCount+	Just bodyMass	= getArgDouble  args ArgBodyMass+	Just epsilon	= getArgDouble	args ArgEpsilon+	Just discSize	= getArgDouble	args ArgDiscSize+	Just startSpeed	= getArgDouble	args ArgStartSpeed++	mMaxSteps	= getArgInt	args ArgMaxSteps+	mFilePath	= getArgString	args ArgDumpFinal++	checkMode x+	 = if not (isJust mWindowSize || isJust mMaxSteps)+		then error "you must specify either --max-steps or --gloss <window size :: Int>"+		else x	++   in	checkMode $+	Config+	{ configWindowSize	= mWindowSize+	, configShouldDrawTree	= shouldDrawTree+	, configRate		= rate+	, configSolverName	= solverName+	, configBodyCount	= bodyCount+	, configBodyMass	= bodyMass+	, configTimeStep	= timeStep+	, configEpsilon		= epsilon+	, configStartDiscSize	= discSize+	, configStartSpeed	= startSpeed+	, configMaxSteps	= mMaxSteps +	, configDumpFinal	= mFilePath }

@@ -0,0 +1,104 @@+{-# LANGUAGE PatternGuards #-}+-- | Drawing the world as a gloss picture.+module Gloss.Draw+	(drawWorld)+where+import Common.World+import Common.Body+import Graphics.Gloss+import qualified Solver.ListBH.Solver		as L+import qualified Data.Vector.Unboxed		as V+++-- | Radius of the circle representing each body.+pointSize :: Float+pointSize	= 4+++-- | Draw the world, and optionally show the Barnes-Hut tree.+--   NOTE: We always show the list version tree, which might not be the one+--         that's actually being used to calculate the accelerations.+--         To display the other trees we'd have to write draw functions for them,+--         or functions to convert them to the list version.+--+drawWorld :: Bool -> World -> Picture+drawWorld shouldDrawTree world+ = let	picPoints	= Color (makeColor 1 1 1 0.4)+			$ Pictures +			$ map drawBody+			$ V.toList +			$ worldBodies world++   	picTree		= drawBHTree+			$ L.buildTree +			$ map massPointOfBody+			$ V.toList +			$ worldBodies world++   in	Pictures +		[ if shouldDrawTree +			then Color (makeColor 0.5 1.0 0.5 0.2) $ picTree+			else Blank+			+		, picPoints ]+++-- | Draw a list version Barnes-Hut tree.+drawBHTree :: L.BHTree -> Picture+drawBHTree bht+ = drawBHTree' 0 bht++drawBHTree' depth bht+ = let	+	-- The bounding box+	L.Box left down right up	= L.bhTreeBox bht+	[left', down', right', up']	= map realToFrac [left, down, right, up]++	picCell		= lineLoop [(left', down'), (left', up'), (right', up'), (right', down')]+++	-- Draw a circle with an area equal to the mass of the centroid.+	centroidX	= realToFrac $ L.bhTreeCenterX bht+	centroidY	= realToFrac $ L.bhTreeCenterY bht+	+	centroidMass	= L.bhTreeMass bht+	circleRadius	= realToFrac $ sqrt (centroidMass / pi)++	midX		= (left' + right') / 2+	midY		= (up'   + down')  / 2++	picCentroid	+	 | _:_	<- L.bhTreeBranch bht+	 , depth >= 1+	 = Color (makeColor 0.5 0.5 1.0 0.4)+		$  Pictures+			[ Line [(midX, midY), (centroidX, centroidY)]+			, Translate centroidX centroidY +			$ ThickCircle+				(circleRadius * 4 / 2) +				(circleRadius * 4) ]+			+	 | otherwise+	 = Blank++	-- The complete picture for this cell.+	picHere		= Pictures [picCentroid, picCell]+		+	-- Pictures of children.+	picSubs		= map (drawBHTree' (depth + 1))+			$ L.bhTreeBranch bht++   in	Pictures (picHere : picSubs)+++-- | Draw a single body.+drawBody :: Body -> Picture+drawBody ((x, y, _), _, _)+	= drawPoint (x, y)+++-- | Draw a point using a filled circle.+drawPoint :: (Double, Double) -> Picture+drawPoint (x, y)+	= Translate (realToFrac x) (realToFrac y) +	$ ThickCircle (pointSize / 2) pointSize

@@ -0,0 +1,151 @@+{-# LANGUAGE ParallelListComp, BangPatterns #-}++import Gloss.MainArgs+import Gloss.Draw+import Gloss.Config++import Common.Dump+import Common.World+import Common.Body+import Common.Util++import Solver+import Timing+import Points2D.Generate++import Graphics.Gloss+import Graphics.Gloss.Interface.Simulate++import System.Environment+import System.Console.ParseArgs+import System.IO.Unsafe+import Control.Monad+import Data.Maybe+import qualified Data.Vector.Unboxed		as V+++main :: IO ()+main  + = do	args	<- parseArgsIO ArgsComplete mainArgs+	+	when (gotArg args ArgHelp)+	 $ usageError args ""++	mainWithArgs args+	++mainWithArgs :: Args MainArg -> IO ()+mainWithArgs args+ = let	config		= loadConfig args++	-- The solver we're using to calculate the acclerations.+	solverName	= configSolverName config+	calcAccels	= fromMaybe (error $ unlines+					[ "unknown solver " ++ show solverName+					, "choose one of "  ++ (show $ map fst solvers) ])+			$ lookup solverName solvers+	+	-- Setup initial world+	vPoints 	= genPointsDisc +				(configBodyCount config)+	 			(0, 0) +				(configStartDiscSize config)++	vBodies		= V.map (setStartVelOfBody $ configStartSpeed config)+			$ V.map (setMassOfBody     $ configBodyMass   config)+			$ V.map (uncurry unitBody) +			$ vPoints++	worldStart	= World+			{ worldBodies	= vBodies+			, worldSteps	= 0 }+				+    in	case configWindowSize config of+	 Just windowSize	-> mainGloss config calcAccels worldStart windowSize+	 Nothing		-> mainBatch config calcAccels worldStart +++-- | Run the simulation in a gloss window.+mainGloss +	:: Config+	-> Solver	-- ^ Fn to calculate accels of each point.+	-> World	-- ^ Initial world.+	-> Int		-- ^ Size of window.+	-> IO ()+	+mainGloss config calcAccels worldStart windowSize+ = let	draw	= drawWorld (configShouldDrawTree config)++	advance _viewport time world	+	 = let	world'	= advanceWorld +				(calcAccels $ configEpsilon config)+				(configTimeStep config)+				world++		-- if we've done enough steps then bail out now.+	   in	case configMaxSteps config of+	 	 Nothing		-> world'+	 	 Just maxSteps+	  	   | worldSteps world' < maxSteps	-> world'+	  +		   -- Gloss doesn't provide a clean way to end the animation...+	  	   | otherwise	+	  	   -> unsafePerformIO (mainEnd (configDumpFinal config) world' >> (error $ "done")) +			`seq` error "advanceWorld: we're finished, stop calling me."++   in	simulateInWindow+		"Barnes-Hut"		-- window name+		(windowSize, windowSize)-- window size+		(10, 10)		-- window position+		black			-- background color+		(configRate config)	-- number of iterations per second+		worldStart		-- initial world+		draw			-- fn to convert a world to a picture+		advance			-- fn to advance the world+++-- | Run the simulation in batch mode, not displaying anything to the screen.+mainBatch+	:: Config+	-> Solver		-- ^ Fn to calculate accels of each point.+	-> World		-- ^ Initial world.+	-> IO ()+	+mainBatch config calcAccels worldStart+ = do+	worldStart `seq` return ()++	(world', tElapsed)+		<- time +		$  let 	world	= mainBatchRun config calcAccels worldStart+		   in	world `seq` return world+					+	putStr $ prettyTime tElapsed+	mainEnd (configDumpFinal config) world'+	++mainBatchRun config calcAccels worldStart + = go worldStart+ where	go !world+ 	  = let world' = advanceWorld+				(calcAccels $ configEpsilon config)+				(configTimeStep config)+				world+	    in case configMaxSteps config of+		Nothing	-> go world'+		Just maxSteps+ 		  | worldSteps world' < maxSteps -> go world'+		  | otherwise 	-> world'+++-- | Called at end of run to dump final world state.+mainEnd +	:: Maybe FilePath	-- ^ Write final bodies to this file.+	-> World		-- ^ Final world state.+	-> IO ()++mainEnd mDumpFinal world+ = do	-- Dump the final world state to file if requested.+	maybe 	(return ()) (dumpWorld world) mDumpFinal++

@@ -0,0 +1,113 @@++module Gloss.MainArgs+	( MainArg(..)+	, mainArgs)+where+import System.Console.ParseArgs++data MainArg+	= ArgHelp+	| ArgSolver+	| ArgMaxSteps+	| ArgTimeStep+	| ArgBodyCount+	| ArgBodyMass+	| ArgEpsilon+	| ArgDiscSize+	| ArgStartSpeed+	+	-- Gloss output+	| ArgGloss+	| ArgRate+	| ArgDrawTree+	+	-- Dump output+	| ArgDumpFinal+	deriving (Eq, Ord, Show)+	+mainArgs :: [Arg MainArg]+mainArgs+ = 	[ Arg	{ argIndex	= ArgHelp+		, argAbbr	= Just 'h'+		, argName	= Just "help"+		, argData	= Nothing+		, argDesc	= "Print this usage help." }++	-- Solver selection.+	, Arg	{ argIndex	= ArgSolver+		, argAbbr	= Just 's'+		, argName	= Just "solver"+		, argData	= argDataDefaulted "name" ArgtypeString "vector-naive"+		, argDesc	= "One of: list-bh, vector-naive, vector-bh, nested-bh. (default vector-naive)" }++	-- Simulation setup.+	, Arg	{ argIndex	= ArgTimeStep+		, argAbbr	= Just 't'+		, argName	= Just "timestep"+		, argData	= argDataDefaulted "Double" ArgtypeDouble 1+		, argDesc	= "Time step between states (default 1)" }++	, Arg	{ argIndex	= ArgBodyCount+		, argAbbr	= Just 'b'+		, argName	= Just "bodies"+		, argData	= argDataDefaulted "Int" ArgtypeInt 200 +		, argDesc	= "Number of bodies in simulation (default 200)" }++	, Arg	{ argIndex	= ArgBodyMass+		, argAbbr	= Just 'm'+		, argName	= Just "mass"+		, argData	= argDataDefaulted "Double" ArgtypeDouble 10+		, argDesc	= "Mass of each body (default 10)" }++	, Arg	{ argIndex	= ArgEpsilon+		, argAbbr	= Just 'e'+		, argName	= Just "epsilon"+		, argData	= argDataDefaulted "Double" ArgtypeDouble 100+		, argDesc	= "Smoothing parameter (default 100)" }+		+	, Arg	{ argIndex	= ArgDiscSize+		, argAbbr	= Just 'd'+		, argName	= Just "disc"+		, argData	= argDataDefaulted "Double" ArgtypeDouble 50+		, argDesc	= "Starting size of disc containing bodies (default 50)" }++	, Arg	{ argIndex	= ArgStartSpeed+		, argAbbr	= Just 'p'+		, argName	= Just "speed"+		, argData	= argDataDefaulted "Double" ArgtypeDouble 0.5+		, argDesc	= "Starting rotation speed of bodies (default 0.5)" }++	-- Termination conditions.+	, Arg	{ argIndex	= ArgMaxSteps+		, argAbbr	= Nothing+		, argName	= Just "max-steps"+		, argData	= argDataOptional "steps" ArgtypeInt+		, argDesc	= "Exit simulation after this many steps." }++	-- Gloss animation output.+	, Arg	{ argIndex	= ArgGloss+		, argAbbr	= Nothing+		, argName	= Just "gloss"+		, argData	= argDataOptional "Int" ArgtypeInt+		, argDesc	= "Animate simulation in window of this size" }++	, Arg	{ argIndex	= ArgRate+		, argAbbr	= Nothing+		, argName	= Just "gloss-rate"+		, argData	= argDataDefaulted "Double" ArgtypeInt 50+		, argDesc	= "(opt. for gloss) Number of steps per second of real time (default 50)" }++	, Arg	{ argIndex	= ArgDrawTree+		, argAbbr	= Nothing+		, argName	= Just "gloss-tree"+		, argData	= Nothing+		, argDesc	= "(opt. for gloss) Draw the Barnes-Hut quad tree"}++	-- Dump points to file+	, Arg	{ argIndex	= ArgDumpFinal+		, argAbbr	= Nothing+		, argName	= Just "dump-final"+		, argData	= argDataOptional "FilePath" ArgtypeString+		, argDesc	= "Dump final body positions and masses to file" }++	]

@@ -0,0 +1,59 @@++-- | Wrappers for the various solvers.+module Solver+	( Solver+	, solvers)+where+import Common.Body++import qualified Solver.ListBH.Solver		as SolverLB++import qualified Data.Vector.Unboxed		as V+import qualified Solver.VectorBH.Solver		as SolverVB+import qualified Solver.VectorNaive.Solver	as SolverVN++import qualified Data.Array.Parallel	as P+import qualified Data.Array.Parallel.PArray	as P+import qualified Solver.NestedBH.Solver		as SolverNB++type Solver	= Double -> V.Vector MassPoint -> V.Vector Accel++solvers :: [(String, Solver)]+solvers+ = 	[ ("list-bh",		calcAccels_lb)+	, ("vector-naive",	calcAccels_vn)+	, ("vector-bh",		calcAccels_vb)+	, ("nested-bh",		calcAccels_nb) ]+++-- | Lists + Barnes-Hut algorithm.+calcAccels_lb	:: Solver+calcAccels_lb epsilon mpts+	= V.fromList+	$ SolverLB.calcAccels epsilon+	$ V.toList mpts+++-- | Vector + Naive algorithm.+calcAccels_vn	:: Solver+calcAccels_vn epsilon+	= SolverVN.calcAccels epsilon +	++-- | Vector + Barnes-Hut algorithm.+calcAccels_vb 	:: Solver+calcAccels_vb epsilon mpts+	= SolverVB.calcAccels epsilon mpts+++-- | Nested Data Parallelism + Barnes-Hut algorithm.+calcAccels_nb	:: Solver+calcAccels_nb epsilon mpts+ = let	+	-- bounds finding isn't vectorised yet.+	(llx, lly, rux, ruy)	= SolverVB.findBounds mpts++	mpts'	= P.fromList $ V.toList mpts+	accels'	= SolverNB.calcAccelsWithBoxPA epsilon llx lly rux ruy mpts'+	+   in	V.fromList $ P.toList accels'

@@ -0,0 +1,161 @@+{-# LANGUAGE BangPatterns, PatternGuards #-}++-- | The list version of the solver also builds the bounding box at every+--   node of the tree, which is good for visualisation.+module Solver.ListBH.Solver+	( MassPoint	(..)+	, BoundingBox	(..)+	, BHTree	(..)+	, calcAccels+	, buildTree)+where+import Common.Body++eClose :: Double+eClose  = square 500++square x = x * x++-- | A rectangular region in 2D space.+data BoundingBox+	= Box+	{ boxLowerLeftX	 :: {-# UNPACK #-} !Double+	, boxLowerLeftY	 :: {-# UNPACK #-} !Double+	, boxUpperRightX :: {-# UNPACK #-} !Double+	, boxUpperRightY :: {-# UNPACK #-} !Double }+	deriving Show+	+-- | The Barnes-Hut tree we use to organise the points.+data BHTree+	= BHT+	{ bhTreeBox	:: {-# UNPACK #-} !BoundingBox+	, bhTreeCenterX	:: {-# UNPACK #-} !Double+	, bhTreeCenterY	:: {-# UNPACK #-} !Double+	, bhTreeMass	:: {-# UNPACK #-} !Double+	, bhTreeBranch	:: ![BHTree] }+	deriving Show+++-- | Compute the acclerations on all these points.+calcAccels :: Double -> [MassPoint] -> [Accel]+calcAccels epsilon mpts+	= map (calcAccel epsilon (buildTree mpts)) mpts+	++-- | Build a Barnes-Hut tree from these points.+buildTree :: [MassPoint] -> BHTree+buildTree mpts+ = let	(llx, lly, rux, ruy)	= findBounds mpts+	box			= Box llx lly rux ruy+   in	buildTreeWithBox box mpts+++-- | Find the coordinates of the bounding box that contains these points.+findBounds :: [MassPoint] -> (Double, Double, Double, Double)+{-# INLINE findBounds #-}+findBounds ((x1, y1, _) : rest1)+ = go x1 y1 x1 y1 rest1	+ where	go !left !right !down !up pts+	 = case pts of+		[]	-> (left, down, right, up)+		(x, y, _) : rest+		 -> let	left'	= min left  x+			right'	= max right x+			down'	= min down  y+			up'	= max up    y+	   	    in	go left' right' down' up' rest+++-- | Given a bounding box that contains all the points, +--   build the Barnes-Hut tree for them.+buildTreeWithBox+	:: BoundingBox		-- ^ bounding box containing all the points.+	-> [MassPoint]		-- ^ points in the box.+	-> BHTree++buildTreeWithBox bb particles+  | length particles <= 1	= BHT bb x y m []+  | otherwise			= BHT bb x y m subTrees+  where	(x, y, m)		= calcCentroid particles+    	(boxes, splitPnts)	= splitPoints bb particles +    	subTrees		= [buildTreeWithBox bb' ps | (bb', ps) <- zip boxes splitPnts]++  +-- | Split massPoints according to their locations in the quadrants.+splitPoints+	:: BoundingBox		-- ^ bounding box containing all the points.+	-> [MassPoint]		-- ^ points in the box.+	-> ( [BoundingBox]	-- +	   , [[MassPoint]])++splitPoints b@(Box llx lly rux  ruy) particles +  | noOfPoints <= 1 = ([b], [particles])+  | otherwise         +  = unzip [ (b,p) | (b,p) <- zip boxes splitPars, length p > 0]+  where+        noOfPoints	= length particles++	-- The midpoint of the parent bounding box.+        (midx,  midy)	= ((llx + rux) / 2.0 , (lly + ruy) / 2.0) ++	-- Split the parent bounding box into four quadrants.+        b1		= Box llx  lly  midx midy+        b2		= Box llx  midy midx  ruy+        b3		= Box midx midy rux   ruy+        b4		= Box midx lly  rux  midy+        boxes		= [b1,   b2,  b3,  b4]++	-- Sort the particles into the smaller boxes.+        lls		= [ p | p <- particles, inBox b1 p ]+        lus		= [ p | p <- particles, inBox b2 p ]+        rus		= [ p | p <- particles, inBox b3 p ]+        rls		= [ p | p <- particles, inBox b4 p ]+        splitPars	= [lls, lus, rus, rls]+++-- | Check if a particle is in box (excluding left and lower border)+inBox:: BoundingBox -> MassPoint -> Bool+{-# INLINE inBox #-}+inBox (Box llx  lly rux  ruy) (px, py, _) +	= (px > llx) && (px <= rux) && (py > lly) && (py <= ruy)+++-- | Calculate the centroid of some points.+calcCentroid :: [MassPoint] -> MassPoint+{-# INLINE calcCentroid #-}+calcCentroid mpts = (sum xs / mass, sum ys / mass, mass)+  where+    mass     = sum   [ m | (_, _, m)  <- mpts ]+    (xs, ys) = unzip [ (m * x, m * y) | (x, y,  m) <- mpts ]   +++-- | Calculate the accelleration of a point due to the points in the given tree.+--   If the distance between the points is less then some small number+--   we set the accel to zero to avoid the acceleration going to infinity+--   and the points escaping the simulation. +--+--   We also use this behavior as a hacky way to discard the acceleration+--   of a point due to interaction with itself.+--+calcAccel:: Double -> BHTree -> MassPoint -> (Double, Double)	+calcAccel !epsilon (BHT _ x y m subtrees) mpt+	| []	<- subtrees+	= accel epsilon mpt (x, y, m)+	+	| not $ isClose mpt x y+	= accel epsilon mpt (x, y, m)++	| otherwise+	= let	(xs, ys)  = unzip [ calcAccel epsilon st mpt | st <- subtrees]+	  in	(sum xs, sum ys) +++-- | If the a point is "close" to a region in the Barnes-Hut tree then we compute+--   the "real" acceleration on it due to all the points in the region, otherwise+--   we just use the centroid as an approximation of all the points in the region.+--+isClose :: MassPoint -> Double -> Double -> Bool+{-# INLINE isClose #-}+isClose (x1, y1, m) x2 y2 +	= (x1-x2) * (x1-x2) + (y1-y2) * (y1-y2) < eClose+

@@ -0,0 +1,154 @@+{-# LANGUAGE ParallelArrays #-}+{-# OPTIONS -fvectorise #-}+module Solver.NestedBH.Solver +	(calcAccelsWithBoxPA)+where+import Data.Array.Parallel+import Data.Array.Parallel.Prelude.Double+import qualified Data.Array.Parallel.Prelude.Int as I+import qualified Prelude++data BoundingBox	+	= Box	Double 		-- lower left  X+		Double 		-- lower left  Y+		Double		-- upper right X+		Double		-- upper right Y++data MassPoint+	= MP 	Double 		-- pos X+		Double		-- pos Y+		Double		-- mass++type Accel	+	= (Double, Double)++data BHTree+	= BHT	Double		-- size of cell+		Double		-- centroid X+		Double		-- centroid Y+		Double		-- centroid mass+		[:BHTree:]	-- children++calcAccelsWithBoxPA+	:: Double+	-> Double -> Double -> Double -> Double+	-> PArray (Double, Double, Double)+	-> PArray (Double, Double)++calcAccelsWithBoxPA epsilon llx lly rux ruy mpts+ = let	mpts'	= [: MP x y m | (x, y, m) <- fromPArrayP mpts :]+	accs'	= calcAccelsWithBox epsilon llx lly rux ruy mpts'+   in	toPArrayP accs'+	++-- | Given the extend of a bounding box containing all the points,+--   calculate the accelerations on all of them.+calcAccelsWithBox+	:: Double+	-> Double -> Double -> Double -> Double+	-> [: MassPoint :]+	-> [: Accel :]++calcAccelsWithBox epsilon llx lly rux ruy mspts+ = accs+ where	accs = [: calcAccel epsilon m tree | m <- mspts :]+	tree = buildTree (Box llx lly rux ruy) mspts+++-- | Build the Barnes-Hut quadtree tree.+buildTree :: BoundingBox -> [: MassPoint :] -> BHTree+buildTree bb particles+ | lengthP particles I.<= 1	= BHT s x y m emptyP+ | otherwise			= BHT s x y m subTrees+ where	(MP x y m)		= calcCentroid particles+	(boxes, splitPnts)	= splitPoints bb particles +    	subTrees		= [:buildTree bb' ps | (bb', ps) <- zipP boxes splitPnts:]+  +	(Box llx lly rux ruy)	= bb+	sx			= rux - llx+	sy			= ruy - lly+	s			= if sx < sy then sx else sy+++-- | Split massPoints according to their locations in the quadrants.+splitPoints+	:: BoundingBox+	-> [: MassPoint :]+	-> ([:BoundingBox:], [:[: MassPoint :]:])++splitPoints b@(Box llx lly rux  ruy) particles +  | noOfPoints I.<= 1 = (singletonP b, singletonP particles)+  | otherwise         +  = unzipP [: (b,p) | (b,p) <- zipP boxes splitPars, lengthP p I.> 0:]+  where	noOfPoints	= lengthP particles+	lls		= [: p | p <- particles, inBox b1 p :]+	lus		= [: p | p <- particles, inBox b2 p :]+	rus		= [: p | p <- particles, inBox b3 p :]+	rls		= [: p | p <- particles, inBox b4 p :]+	b1		= Box llx  lly  midx midy+	b2		= Box llx  midy midx  ruy+	b3		= Box midx midy rux   ruy+ 	b4		= Box midx lly  rux  midy+	boxes		= singletonP b1  +:+ singletonP b2  +:+ singletonP b3 +:+ singletonP b4 +	splitPars	= singletonP lls +:+ singletonP lus +:+ singletonP rus +:+ singletonP rls+	(midx,  midy)	= ((llx + rux) / 2.0 , (lly + ruy) / 2.0) +++-- | Checks if particle is in box (excluding left and lower border)+inBox :: BoundingBox -> MassPoint -> Bool+inBox (Box llx  lly rux  ruy) (MP px  py  _) + 	= (px > llx) && (px <= rux) && (py > lly) && (py <= ruy)+++-- | Calculate the centroid of some points.+calcCentroid:: [:MassPoint:] -> MassPoint+calcCentroid mpts + = MP  (sumP xs / mass) (sumP ys / mass) mass+ where	mass     = sumP [: m | MP _ _ m  <- mpts :]+	(xs, ys) = unzipP [: (m * x, m * y) | MP x y m <- mpts :]   +++-- | Calculate the accelleration of a point due to the points in the given tree.+calcAccel :: Double -> MassPoint -> BHTree -> (Double, Double)+calcAccel epsilon mpt (BHT s x y m subtrees)+	| lengthP subtrees I.== 0+	= accel epsilon mpt (MP x y m)++	| isFar mpt s x y +	= accel epsilon mpt (MP x y m)++	| otherwise+	= let	(xs, ys) = unzipP [: calcAccel epsilon mpt st | st <- subtrees :]+	  in	(sumP xs, sumP ys)+++-- | Calculate the acceleration on a point due to some other point.+accel 	:: Double 	-- ^ If the distance between the points is smaller than this+			--   then ignore the forces between them.+	-> MassPoint	-- ^ The point being acclerated.+	-> MassPoint	-- ^ Neibouring point.+	-> Accel++accel epsilon (MP x1 y1 _) (MP x2 y2 m)  + = (aabs * dx / r , aabs * dy / r)  + where	rsqr = (dx * dx) + (dy * dy) + epsilon+	r    = sqrt rsqr +	dx   = x1 - x2 +	dy   = y1 - y2 +	aabs = m / rsqr +++-- | If the point is far from a cell in the tree then we can use+--   it's centroid as an approximation of all the points in the region.+isFar 	:: MassPoint 	-- point being accelerated+	-> Double	-- size of region+	-> Double	-- position of center of mass of cell+	-> Double	-- position of center of mass of cell+	-> Bool++isFar (MP x1 y1 m) s x2 y2 + = let	dx	= x2 - x1+	dy	= y2 - y1+	dist	= sqrt (dx * dx + dy * dy)+   in	(s / dist) < 1+

@@ -0,0 +1,163 @@+{-# LANGUAGE BangPatterns, PatternGuards #-}++-- | The list version of the solver also builds the bounding box at every+--   node of the tree, which is good for visualisation.+module Solver.VectorBH.Solver+	( MassPoint	(..)+	, BoundingBox	(..)+	, BHTree	(..)+	, calcAccels+	, buildTree+	, findBounds)+where+import Common.Body+import Data.Vector.Unboxed			(Vector)+import qualified Data.Vector.Unboxed		as V+	+type BoundingBox+	= (Double, Double, Double, Double)+	+sizeOfBox :: BoundingBox -> Double+{-# INLINE sizeOfBox #-}+sizeOfBox (llx, lly, rux, ruy)+	= min (abs (rux - llx)) (abs (ruy - lly))+++-- | The Barnes-Hut tree we use to organise the points.+data BHTree+	= BHT+	{ bhTreeSize	:: {-# UNPACK #-} !Double	-- minimum of hight and width of cell+	, bhTreeCenterX	:: {-# UNPACK #-} !Double+	, bhTreeCenterY	:: {-# UNPACK #-} !Double+	, bhTreeMass	:: {-# UNPACK #-} !Double+	, bhTreeBranch	:: ![BHTree] }+	deriving Show+++-- | Compute the acclerations on all these points.+calcAccels :: Double -> Vector MassPoint -> Vector Accel+calcAccels epsilon mpts+ = V.map (calcAccel epsilon (buildTree mpts)) mpts+++-- | Build a Barnes-Hut tree from these points.+buildTree :: Vector MassPoint -> BHTree+buildTree mpts+	= buildTreeWithBox (findBounds mpts) mpts+++-- | Find the coordinates of the bounding box that contains these points.+findBounds :: Vector MassPoint -> (Double, Double, Double, Double)+{-# INLINE findBounds #-}+findBounds bounds+ = V.foldl' acc (x1, y1, x1, y1) bounds+ where+	(x1, y1, _)	= bounds V.! 0++	acc (!llx, !lly, !rux, !ruy) (x, y, _)+	 = let	!llx'	= min llx  x+		!lly'	= min lly  y+		!rux'	= max rux  x+		!ruy'	= max ruy  y+	   in	(llx', lly', rux', ruy')+++-- | Given a bounding box that contains all the points, +--   build the Barnes-Hut tree for them.+buildTreeWithBox+	:: BoundingBox		-- ^ bounding box containing all the points.+	-> Vector MassPoint	-- ^ points in the box.+	-> BHTree++buildTreeWithBox bb mpts+  | V.length mpts <= 1		= BHT s x y m []+  | otherwise			= BHT s x y m subTrees+  where	+	s			= sizeOfBox bb+	(x, y, m)		= calcCentroid   mpts+    	(boxes, splitPnts)	= splitPoints bb mpts+    	subTrees		= [buildTreeWithBox bb' ps+					| (bb', ps) <- zip boxes splitPnts]++  +-- | Split massPoints according to their locations in the quadrants.+splitPoints+	:: BoundingBox		-- ^ bounding box containing all the points.+	-> Vector MassPoint	-- ^ points in the box.+	-> ( [BoundingBox]	-- +	   , [Vector MassPoint])++splitPoints b@(llx, lly, rux, ruy) mpts+  | noOfPoints <= 1 = ([b], [mpts])+  | otherwise         +  = unzip [ (b,p) +		| (b,p) <- zip boxes splitPars+		, V.length p > 0]+  where+        noOfPoints	= V.length mpts++	-- The midpoint of the parent bounding box.+        (midx,  midy)	= ((llx + rux) / 2.0 , (lly + ruy) / 2.0) ++	-- Split the parent bounding box into four quadrants.+        b1		= (llx,  lly,  midx,  midy)+        b2		= (llx,  midy, midx,  ruy)+        b3		= (midx, midy, rux,   ruy)+        b4		= (midx, lly,  rux,   midy)+        boxes		= [b1,   b2,   b3,   b4]++	-- Sort the particles into the smaller boxes.+        lls		= V.filter (inBox b1) mpts+        lus		= V.filter (inBox b2) mpts+        rus		= V.filter (inBox b3) mpts+        rls		= V.filter (inBox b4) mpts+        splitPars	= [lls, lus, rus, rls]+++-- | Check if a particle is in box (excluding left and lower border)+inBox:: BoundingBox -> MassPoint -> Bool+{-# INLINE inBox #-}+inBox (llx, lly, rux, ruy) (px, py, _) +	= (px > llx) && (px <= rux) && (py > lly) && (py <= ruy)+++-- | Calculate the centroid of some points.+calcCentroid :: Vector MassPoint -> MassPoint+{-# INLINE calcCentroid #-}+calcCentroid mpts +  = (V.sum xs / mass, V.sum ys / mass, mass)+  where	mass     = V.sum   $ V.map (\(_, _, m) -> m) mpts+	(xs, ys) = V.unzip $ V.map (\(x, y, m) -> (m * x, m * y)) mpts+++-- | Calculate the accelleration of a point due to the points in the given tree.+calcAccel:: Double -> BHTree -> MassPoint -> (Double, Double)+calcAccel !epsilon (BHT s x y m subtrees) mpt+	| []	<- subtrees+	= accel epsilon mpt (x, y, m)+	+	| isFar mpt s x y+	= accel epsilon mpt (x, y, m)++	| otherwise+	= let	(xs, ys)  = unzip [ calcAccel epsilon st mpt | st <- subtrees]+	  in	(sum xs, sum ys) +++-- | If the point is far from a cell in the tree then we can use+--   it's centroid as an approximation of all the points in the region.+--+isFar 	:: MassPoint 	-- point being accelerated+	-> Double	-- size of region+	-> Double	-- position of center of mass of cell+	-> Double	-- position of center of mass of cell+	-> Bool++{-# INLINE isFar #-}+isFar (x1, y1, m) s x2 y2 + = let	!dx	= x2 - x1+	!dy	= y2 - y1+	!dist	= sqrt (dx * dx + dy * dy)+   in	(s / dist) < 1++

@@ -0,0 +1,18 @@++-- | Naive n^2 computation of accelerations.+module Solver.VectorNaive.Solver+	(calcAccels)+where+import Common.Body+import Data.Vector.Unboxed		(Vector)+import qualified Data.Vector.Unboxed	as V++-- | Calculate accelerations on these point in a naive O(n^2) way+calcAccels :: Double -> Vector MassPoint -> Vector Accel+calcAccels epsilon mps+	= V.map (calcAccel epsilon mps) mps++calcAccel :: Double -> Vector MassPoint -> MassPoint -> Accel+calcAccel epsilon mps mp+ = let	(xs, ys)	= V.unzip $ V.map (accel epsilon mp) mps+   in	(V.sum xs, V.sum ys)

@@ -0,0 +1,59 @@++import QuickHullVect+import SVG+import Timing+import Points2D.Generate+import Points2D.Types+import System.Environment+import Data.Array.Parallel.PArray	as P++main :: IO ()+main + = do	args	<- getArgs+	case args of+	  [pointCount]	+	    -> run (read pointCount) Nothing+	+	  [pointCount, fileSVG]+	    -> run (read pointCount) (Just fileSVG)++	  _ -> do+		putStr usage+		return ()+++-- | Command line usage information.+usage :: String+usage	= unlines+	[ "Usage: quickhull <points> [out.svg]"	]+++-- | Run the benchmark.+run 	:: Int 			-- ^ How many points to use.+	-> Maybe String 	-- ^ File name to dump an SVG of the output to.+	-> IO ()+	+run pointCount mFileSVG+ = do+	vPoints	<- pointsPArrayOfUArray+		$ genPointsDisc pointCount (400, 400) 350 ++	-- Force points to create the input vector.+	vPoints `seq` return ()++	-- Compute the convex hull.+	(vHull, tElapsed)+		<- time +		$  let 	vHull	= quickhullPA vPoints+		   in	vHull `seq` return vHull+					+	-- Print how long it took.+	putStr $ prettyTime tElapsed+	+	-- If we were asked for an SVG then write it out to file.+	maybe 	(return ())+	 	(\fileSVG -> +			writeFile fileSVG+			 $ makeSVG 	(roundPoints $ P.toList vPoints) +					(roundPoints $ P.toList vHull))+		mFileSVG

@@ -0,0 +1,42 @@+{-# LANGUAGE ParallelArrays #-}+{-# OPTIONS -fvectorise #-}++module QuickHullVect (quickhullPA) where+import Points2D.Types+import Data.Array.Parallel+import Data.Array.Parallel.Prelude.Double+import qualified Data.Array.Parallel.Prelude.Int as Int+import qualified Prelude as P+++distance :: Point -> Line -> Double+distance (xo, yo) ((x1, y1), (x2, y2))+  = (x1-xo) * (y2 - yo) - (y1 - yo) * (x2 - xo)+++hsplit :: [:Point:] -> Line -> [:Point:]+hsplit points line@(p1, p2)+  | lengthP packed Int.== 0 = [:p1:]+  | otherwise+  = concatP [: hsplit packed ends | ends <- [:(p1, pm), (pm, p2):] :]+  where+    cross  = [: distance p line | p <- points :]+    packed = [: p | (p,c) <- zipP points cross, c > 0.0 :]+    pm     = points !: maxIndexP cross+++quickHull :: [:Point:] -> [:Point:]+quickHull points+  | lengthP points Int.== 0 = points+  | otherwise+  = concatP [: hsplit points ends | ends <- [: (minx, maxx), (maxx, minx) :] :]+  where+    xs   = [: x | (x, y) <- points :]+    minx = points !: minIndexP xs+    maxx = points !: maxIndexP xs+++quickhullPA :: PArray Point -> PArray Point+{-# NOINLINE quickhullPA #-}+quickhullPA ps = toPArrayP (quickHull (fromPArrayP ps))+

@@ -0,0 +1,34 @@++module SVG where++-- Making a SVG diagram of the points and hull+makeSVG :: [(Int, Int)] -> [(Int, Int)] -> String+makeSVG points hull+	= unlines+	$  [ "<svg width=\"100%\" height=\"100%\" version=\"1.1\" xmlns=\"http://www.w3.org/2000/svg\">" ]+	++ [svgPolygon hull]+	++ map svgPoint points+	++ map svgPointHull hull+	++ ["</svg>"]++svgPolygon  :: [(Int, Int)] -> String+svgPolygon points+	=  "<polygon"+	++ " points=\"" ++ (concat [show x ++ "," ++ show y ++ " " | (x, y) <- points]) ++ "\"" +	++ " style=\"fill:#d0d0ff;stroke:#000000;stroke-width:1\""+	++ "/>"++svgPoint :: (Int, Int) -> String+svgPoint (x, y)+	= "<circle cx=\"" ++ show x ++ "\" cy=\"" ++ show y ++ "\" r=\"0.5\""+	++ " style=\"stroke:#000000\""+	++ "/>"++svgPointHull :: (Int, Int) -> String+svgPointHull (x, y)+	= "<circle cx=\"" ++ show x ++ "\" cy=\"" ++ show y ++ "\" r=\"1\""+	++ " style=\"stroke:#ff0000\""+	++ "/>"+	+roundPoints :: [(Double, Double)] -> [(Int, Int)]+roundPoints ps = [(round x, round y) | (x, y) <- ps]

@@ -0,0 +1,71 @@+{-# LANGUAGE PatternGuards #-}++import qualified QuickHullVector+import qualified QuickHullIO+import qualified QuickHullSplit+import Points2D.Types+import Points2D.Generate+import Timing+import SVG++import System.Environment+import Data.Function+import qualified Data.Vector.Unboxed	as V+import Data.Vector.Unboxed		(Vector)+import Control.Monad++algs = 	[ ("vector",	(\v -> return $ QuickHullVector.quickHull v))+	, ("io",	QuickHullIO.quickHull)+	, ("split",	(\v -> return $ QuickHullSplit.quickHull v)) ]++parseArgs args+	| [alg, strCount]	<- args+	, Just fun 		<- lookup alg algs+	= Just (fun, read strCount, Nothing)++	| [alg, strCount, file]	<- args+	, Just fun 		<- lookup alg algs+	= Just (fun, read strCount, Just file)++	| otherwise+	= Nothing+++main :: IO ()+main+ = do	argStrs		<- getArgs+	case parseArgs argStrs of+	 Just args	-> run args+	 _ 		-> putStr $ unlines+				[ "usage: quickhull <alg> <points> [out.svg]"+				, "   algs: " ++ (show $ map fst algs) ++ "\n" ]++run ::	( Vector Point -> IO (Vector Point)+	, Int+	, Maybe FilePath)+    ->	IO ()+	+run (fun, pointCount, mFileSVG) + = do+	let vPoints	= genPointsDisc pointCount (400, 400) 350 ++	-- Force points to create the input vector.+	V.force vPoints `seq` return ()++	-- Compute the convex hull.+	(vHull, tElapsed)+	 	<- time+		$  do	vHull	<- fun vPoints+		     	vHull `seq` return vHull+		+	-- Print how long it took	+	putStr $ prettyTime tElapsed++	-- If we were asked for an SVG then write it out to file.+	maybe 	(return ())+	 	(\fileSVG -> +			writeFile fileSVG+			 $ makeSVG 	(roundPoints $ V.toList vPoints) +					(roundPoints $ V.toList vHull))+		mFileSVG+

@@ -0,0 +1,269 @@+{-# LANGUAGE BangPatterns, PatternGuards, RankNTypes #-}++module QuickHullIO+	(quickHull)+where+import Data.Function+import Control.Monad+import Control.Exception+import Control.Concurrent+import Control.Concurrent.MVar+import Control.Monad.ST+import GHC.Conc+import Data.IORef+import Data.List+import Data.Ord+import Data.Vector.Unboxed			(Vector)+import qualified Data.Vector.Unboxed		as V+import qualified Data.Vector.Unboxed.Mutable	as MV+import qualified Data.Vector.Generic		as G+import Debug.Trace++type Point	= (Double, Double)+type Line	= (Point, Point)+++-- | Compute the convex hull of a vector of points.+quickHull :: Vector Point -> IO (Vector Point)+quickHull !points+  | V.length points == 0	+  = return points++  | otherwise+  = do	-- Find the left and right-most points.+	let (minx, maxx) 	= minmax points++	-- Hull points get written to the vector in this IORef.+	hullRef	<- newIORef V.empty++	-- Fork off computations to handle half of the points each.+	-- For uniformly distributed points this first iteration takes most of the time.+	parIO 	[ hsplit hullRef points minx maxx+		, hsplit hullRef points maxx minx]++	-- Grab the finished hull points.+	hull	<- readIORef hullRef++	-- We've got the hull points, but they can appear in arbitrary order.+	-- Do a rubbish via-lists merge phase so that they appear clockwise around the edge.+	-- This isn't too expensive if there aren't many points on the hull.+	let (above, below) +		= V.unstablePartition +			(\p -> distance minx maxx p > 0)+			hull+	+	let aboveSorted	= V.fromList $ sortBy (comparing fst) $ V.toList above+	let belowSorted	= V.fromList $ sortBy (comparing fst) $ V.toList below+	let hull' = aboveSorted V.++ V.reverse belowSorted++	return hull'+	++hsplit :: IORef (Vector Point) -> Vector Point -> Point -> Point -> IO ()+{-# INLINE hsplit #-}+hsplit hullRef !points !p1@(!p1X, !p1Y) !p2@(!p2X, !p2Y)+	-- we've found one.+	| V.length packed == 0+	= addHullPoint hullRef p1+	+	-- do the two new segments in parallel.+	| V.length packed > 1000+	= parIO+		[ hsplit hullRef packed p1 pm+		, hsplit hullRef packed pm p2 ]+		+	| otherwise+	= do	hsplit hullRef packed p1 pm+		hsplit hullRef packed pm p2++	where	(packed, pm)	= parPackPoints points p1X p1Y p2X p2Y+	++-- | Copy points from the input vector that are on the left of the line into a+--	new buffer. While we're doing this, determine the point that is furthest+--	from the line.+--+--	If we have a big enough vector then split it in two and do both halves+--	in parallel. Doing this requires a copy afterwards to join the two+--	results back together. It's a trade off between decreased FP load and +--	increased memory traffic. +--+parPackPoints +	:: Vector Point +	-> Double -> Double+	-> Double -> Double+	-> ( Vector Point+	   , Point)+	+{-# INLINE parPackPoints #-}+parPackPoints !points !p1X !p1Y !p2X !p2Y+ |   numCapabilities == 1+  || V.length points < 1000+ = packPoints p1X p1Y p2X p2Y points++ | otherwise+ = let	+	numSegments	= numCapabilities++	-- Total number of points to process.+	lenPoints	= V.length points++	-- How many points to process in each segment.+	lenSeg		= lenPoints `div` numSegments++ 	-- If the total number of points doesn't divide evenly into segments+	-- then there may be an odd number. Make sure to get the rest into the last segment.+	splitPacked count ixStart +	    | count == 0	= []++	    | count == 1	+	    = let points'		= V.unsafeSlice ixStart (lenPoints - ixStart) points+	      	  result@(packed', _)	= packPoints p1X p1Y p2X p2Y points'+	      in  packed' `pseq` (result : [])++	    | otherwise	+	    = let points'		= V.unsafeSlice ixStart lenSeg points+	          result@(packed', _)	= packPoints p1X p1Y p2X p2Y points'+		  rest			= splitPacked (count - 1) (ixStart + lenSeg)+	      in  packed' `par` rest `par` (result : rest)++	results	= splitPacked numSegments 0+	vResult	= concatVectors $ map fst results+	pMax	= selectFurthest p1X p1Y p2X p2Y results+	+   in	(vResult, pMax)+++selectFurthest + 	:: Double -> Double +	-> Double -> Double+	-> [(Vector Point, Point)] +	-> Point+	+selectFurthest !p1X !p1Y !p2X !p2Y ps+ = go (0, 0) 0 ps++ where	go pMax !distMax []	+	 = pMax++	go pMax !distMax ((packed, pm):rest)+	 | V.length packed == 0+	 = go pMax distMax rest+	+	 | otherwise+	 , dist		<-  distance (p1X, p1Y) (p2X, p2Y) pm +  	 = if dist > distMax+		then go pm   dist    rest+		else go pMax distMax rest+++packPoints +	:: Double -> Double 		-- First point on dividing line.+	-> Double -> Double 		-- Second point on dividing line.+	-> Vector Point 		-- Source points.+	-> ( Vector Point		-- Packed vector containing only points on the left of the line.+	   , Point)			-- The point on the left that was furthest from the line.++{-# INLINE packPoints #-}+packPoints !p1X !p1Y !p2X !p2Y !points + = let+	result	+	 = G.create + 	 $ do	packed		 <- MV.new (V.length points + 1)+		(pMax, ixPacked) <- fill points packed p1X p1Y p2X p2Y 0 0++		-- We stash the maximum point on the end of the vector to get+		-- it through the create call.+		MV.unsafeWrite packed ixPacked pMax+		return $ MV.unsafeSlice 0 (ixPacked + 1) packed+	+   in	( V.unsafeSlice 0 (V.length result - 1) result+	, result `V.unsafeIndex` (V.length result - 1))+			++fill 	:: forall s+	.  Vector Point 		-- Source points.+	-> MV.MVector s Point 		-- Vector to write packed points into.+	-> Double -> Double 		-- First point on dividing line.+	-> Double -> Double		-- Second poitn on dividing line.+	-> Int 				-- Index into source points to start reading from.+	-> Int				-- Index into packed points to start writing to.+	-> ST s +		( Point			-- Furthest point from the line that was found.+		, Int)			-- The number of packed points written.++{-# INLINE fill #-}+fill !points !packed !p1X !p1Y !p2X !p2Y !ixPoints' !ixPacked'+ = go (0, 0) 0 ixPoints' ixPacked'+ where go pMax !distMax !ixPoints !ixPacked+	| ixPoints >= V.length points	+	= do	return (pMax, ixPacked)+		+	| p	<- points `V.unsafeIndex` ixPoints+	, d	<- distance (p1X, p1Y) (p2X, p2Y) p+	, d > 0+	= do	MV.unsafeWrite packed ixPacked p+		if d > distMax+		 then	go p    d       (ixPoints + 1) (ixPacked + 1)+		 else	go pMax distMax (ixPoints + 1) (ixPacked + 1)+			+	| otherwise+	= go pMax distMax (ixPoints + 1) ixPacked+++minmax :: Vector Point -> (Point, Point)+{-# INLINE minmax #-}+minmax !vec+ = go first first 0+ where	first	= vec V.! 0++	go pMin@(!minX, !minY) pMax@(!maxX, !maxY) !ix+	  | ix >= V.length vec	= (pMin, pMax)++	  | (x, y)	<- vec `V.unsafeIndex` ix+	  = if       x < minX then go (x, y) pMax   (ix + 1)+	    else if  x > maxX then go pMin   (x, y) (ix + 1)+	    else go pMin pMax (ix + 1)+	++distance :: Point -> Point -> Point -> Double+{-# INLINE distance #-}+distance (x1, y1) (x2, y2) (xo, yo)+  = (x1-xo) * (y2 - yo) - (y1 - yo) * (x2 - xo)+++addHullPoint :: IORef (Vector Point) -> Point -> IO ()+addHullPoint hullRef p+ = atomicModifyIORef hullRef+ $ \hull -> (V.singleton p V.++ hull, ())++++-- Can't find an equivalent for this in Control.Concurrent.+parIO :: [IO ()] -> IO ()+parIO stuff+ = do	mVars	<- replicateM (length stuff) newEmptyMVar+	zipWithM_ (\c v -> forkIO $ c `finally` putMVar v ()) stuff mVars+	mapM_ readMVar mVars+	++-- We really want a function in the vector library for this.+concatVectors :: [Vector Point] -> Vector Point+{-# NOINLINE concatVectors #-}+concatVectors vectors+ = G.create+ $ do	let len	= sum $ map V.length vectors+	vOut	<- MV.new len+	go vectors vOut 0+	return vOut++ where	{-# INLINE go #-}+	go [] _ _+	 = return ()++	go (vSrc:vsSrc) vDest !ixStart	+	 = do	let lenSrc	= V.length vSrc+		let vDestSlice	= MV.unsafeSlice ixStart lenSrc vDest+		V.copy vDestSlice vSrc+		go vsSrc vDest (ixStart + lenSrc) +

@@ -0,0 +1,38 @@++module QuickHullSplit+	(quickHull)+where+import qualified Data.Vector.Unboxed	as V+import Data.Vector.Unboxed		(Vector)++quickHull :: Vector (Double, Double) -> Vector (Double, Double)+quickHull vv+	= uncurry V.zip $ quickhull $ V.unzip vv++quickhull :: (Vector Double, Vector Double) -> (Vector Double, Vector Double)+{-# NOINLINE quickhull #-}+quickhull (xs, ys) = xs' `seq` ys' `seq` (xs',ys')+    where+      (xs',ys') = V.unzip+                $ hsplit points pmin pmax V.++ hsplit points pmax pmin++      imin = V.minIndex xs+      imax = V.maxIndex xs++      points = V.zip xs ys+      pmin   = points V.! imin+      pmax   = points V.! imax+++      hsplit points p1 p2+        | V.length packed < 2 = p1 `V.cons` packed+        | otherwise = hsplit packed p1 pm V.++ hsplit packed pm p2+        where+          cs     = V.map (\p -> cross p p1 p2) points+          packed = V.map fst+                 $ V.filter (\t -> snd t > 0)+                 $ V.zip points cs++          pm     = points V.! V.maxIndex cs++      cross (x,y) (x1,y1) (x2,y2) = (x1-x)*(y2-y) - (y1-y)*(x2-x)

@@ -0,0 +1,80 @@+{-# LANGUAGE BangPatterns, PatternGuards #-}++module QuickHullVector+	(quickHull)+where+import Points2D.Types+import Data.Function+import Data.Vector.Unboxed		as V+import Data.Vector.Unboxed.Mutable	as MV+import Data.Vector.Unboxed		(Vector)+import qualified Data.Vector.Generic	as G+import System.IO.Unsafe++++distance :: Point -> Point -> Point -> Double+{-# INLINE distance #-}+distance (x1, y1) (x2, y2) (xo, yo)+  = (x1-xo) * (y2 - yo) - (y1 - yo) * (x2 - xo)+++hsplit :: Vector Point -> Point -> Point -> Vector Point+{-# INLINE hsplit #-}+hsplit !points !p1@(!p1X, !p1Y) !p2@(!p2X, !p2Y)+ = let !packed	= packPoints points p1X p1Y p2X p2Y+   in  if V.length packed == 1+	 then V.singleton p1+	 else 	let pm		= packed `V.unsafeIndex` (V.length packed - 1)+		    packed'	= V.unsafeSlice 0 (V.length packed - 1) packed+		in  hsplit packed' p1 pm V.++ hsplit packed' pm p2+++packPoints :: Vector Point -> Double -> Double -> Double -> Double -> Vector Point+{-# INLINE packPoints #-}+packPoints !points !p1X !p1Y !p2X !p2Y+ = G.create + $ do	packed	<- MV.new (V.length points + 1)+	+	-- stash the furthest point on the end of the returned vector.	+	let fill !pMax !distMax !ixPoints !ixPacked+		| ixPoints >= V.length points	+		= do	MV.unsafeWrite packed ixPacked pMax+			return $ MV.unsafeSlice 0 (ixPacked + 1) packed++		| p	<- points `V.unsafeIndex` ixPoints+		, d	<- distance (p1X, p1Y) (p2X, p2Y) p+		, d > 0+		= do	MV.unsafeWrite packed ixPacked p+			if d > distMax+			 then	fill p    d       (ixPoints + 1) (ixPacked + 1)+			 else	fill pMax distMax (ixPoints + 1) (ixPacked + 1)+			+		| otherwise+		= fill pMax distMax (ixPoints + 1) ixPacked+			+	fill (0, 0) 0 0 0+++quickHull :: Vector Point -> Vector Point+quickHull !points+  	| V.length points == 0	= points++	| (minx, maxx) 		<- minmax points+	= hsplit points minx maxx V.++ hsplit points maxx minx+++minmax :: Vector Point -> (Point, Point)+{-# INLINE minmax #-}+minmax !vec+ = go first first 0+ where	first	= vec V.! 0++	go pMin@(!minX, !minY) pMax@(!maxX, !maxY) !ix+	  | ix >= V.length vec	= (pMin, pMax)++	  | (x, y)	<- vec `V.unsafeIndex` ix+	  = if       x < minX then go (x, y) pMax   (ix + 1)+	    else if  x > maxX then go pMin   (x, y) (ix + 1)+	    else go pMin pMax (ix + 1)+	

@@ -0,0 +1,50 @@++import QuickSortVect+import Timing+import Randomish+import System.Environment+import Data.Vector.Unboxed		(Vector)+import Data.Array.Parallel	as P+import Data.Array.Parallel.PArray	as P+import qualified Data.Vector.Unboxed	as V+import Data.Maybe++main :: IO ()+main + = do	args	<- getArgs+	case args of+	  [len]	-> run (read len) +	  _	-> putStr usage+++-- | Command line usage information.+usage :: String+usage	= unlines+	[ "Usage: quicksort <length>" ]+++-- | Run the benchmark.+run :: Int -> IO ()+run len+ = do	-- Create the input vector+	let vInts 	= fromUArrPA'+			$ randomishDoubles len 0 1 1234++	vInts `seq` return ()++	-- Compute the convex hull.+	(vSorted, tElapsed)+		<- time +		$  let 	vSorted	= quicksortPA vInts+		   in	vSorted `seq` return vSorted+					+	-- Print how long it took.+	putStr $ prettyTime tElapsed+	+	-- Check they're really sorted.+	print	$ isSorted $ P.toUArrPA vSorted+	++-- | Check if a vector is sorted (monotonically increasing)+isSorted :: Vector Double -> Bool+isSorted = isJust . V.fold1M (\x y -> if y >= x then Just y else Nothing)

@@ -0,0 +1,25 @@+{-# LANGUAGE ParallelArrays #-}+{-# OPTIONS -fvectorise #-}+{-# OPTIONS -fno-spec-constr-count #-}+module QuickSortVect (quicksortPA) where++import Data.Array.Parallel+import Data.Array.Parallel.Prelude.Double+import qualified Data.Array.Parallel.Prelude.Int as I++import qualified Prelude++quicksortPA:: PArray Double -> PArray Double +quicksortPA xs = toPArrayP  (qsortVect' (fromPArrayP xs))++qsortVect':: [: Double :] -> [: Double :]+{-# NOINLINE qsortVect' #-}+qsortVect' xs | lengthP xs I.<=  1 = xs+              | otherwise =+  let p  = xs !: (lengthP xs `I.div` 2)+      ls = [:x | x <- xs, x < p:]+      gs = [:x | x <- xs, x > p:]++      ss = mapP qsortVect' [:ls, gs:]+ in+ (ss !: 0) +:+ [:x | x <- xs, x == p:] +:+ (ss !: 1)

@@ -0,0 +1,102 @@+{-# LANGUAGE ScopedTypeVariables #-}++import Timing+import SMVMVectorised+import System.IO+import Foreign.Storable+import Foreign.Marshal.Alloc+import Data.Array.Parallel.PArray	as P+import Data.Array.Parallel+import qualified Data.Array.Parallel.Unlifted as U+import System.Environment+import Control.Exception (evaluate)++main :: IO ()+main + = do	args	<- getArgs+	case args of+	  [fileName] -> run fileName+	  _	     -> usage+	+usage	+ = putStr $ unlines+	[ "usage: smvm <file>" ]++run fileName+ = do	(matrix, vector) <- loadPA fileName++	matrix `seq` return ()+	vector `seq` return ()++	-- Multiply sparse matrix by the dense vector.+	(vResult, tElapsed)+	 <- time $ let result	= smvmPA matrix vector+		   in  P.nf result `seq` return result+					+	-- Print how long it took.+	putStr $ prettyTime tElapsed++	-- Print some info about the test setup.+	putStrLn $ "vector length   = " ++ show (U.length (toUArrPA vector))+--	putStrLn $ "matrix height   = " ++ show (U.length (toUArrPA matrix))+	+	+	-- Print checksum of resulting vector.+	putStrLn $ "result sum      = " ++ show (U.sum (toUArrPA vResult))++++-- | Load a test file containing a sparse matrix and dense vector.+loadPA 	:: String 				-- ^ filename.+	-> IO  ( PArray (PArray (Int, Double))	-- sparse matrix+	       , PArray Double)			-- dense vector++loadPA fileName+ = do 	(segd, arrMatrixElems, arrVector) <- loadUArr fileName++    	let paMatrix	= nestUSegdPA' segd (fromUArrPA_2' arrMatrixElems)+	let paVector	= fromUArrPA' arrVector+	return (paMatrix, paVector)+++-- | Load a test file containing a sparse matrix and dense vector.+loadUArr :: String				-- ^ filename+	 -> IO ( U.Segd				-- segment descriptor saying what array elements+						--    belong to each row of the matrix.+	       , U.Array (Int, Double)		-- column indices and matrix elements+	       , U.Array Double)		-- the dense vector++loadUArr fname + = do	h <- openBinaryFile fname ReadMode++	-- check magic numbers at start of file to guard against word-size screwups.+	alloca $ \ptr -> do+		hGetBuf h ptr (sizeOf (undefined :: Int))+		magic1 :: Int	<- peek ptr+		hGetBuf h ptr (sizeOf (undefined :: Int))+		magic2	:: Int <- peek ptr+		if magic1 == 0xc0ffee00 Prelude.&& magic2 == 0x12345678 +			then return ()+			else error $ "bad magic in " ++ fname++	-- number of elements in each row of the matrix.+	lengths <- U.hGet h++	-- indices of all the elements.+	indices <- U.hGet h++	-- values of the matrix elements.+	values  <- U.hGet h++	-- the dense vector.+	vector  <- U.hGet h++	evaluate lengths+	evaluate indices+	evaluate values+	evaluate vector++	let segd    = U.lengthsToSegd lengths+	    matrix  = U.zip indices values++	return (segd, matrix, vector)

@@ -0,0 +1,18 @@++{-# LANGUAGE ParallelArrays #-}+{-# OPTIONS -fvectorise #-}+module SMVMVectorised (smvmPA) where++import Data.Array.Parallel+import Data.Array.Parallel.Prelude.Double as D+import Data.Array.Parallel.Prelude.Int    as I++import qualified Prelude as P++smvmPA :: PArray (PArray (Int, Double)) -> PArray Double -> PArray Double+{-# NOINLINE smvmPA #-}+smvmPA m v = toPArrayP (smvm (fromNestedPArrayP m) (fromPArrayP v))++smvm :: [:[: (Int, Double) :]:] -> [:Double:] -> [:Double:]+smvm m v = [: D.sumP [: x D.* (v !: i) | (i,x) <- row :] | row <- m :]+