diff --git a/LICENSE b/LICENSE new file mode 100644 --- /dev/null +++ b/LICENSE @@ -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. diff --git a/Setup.hs b/Setup.hs new file mode 100644 --- /dev/null +++ b/Setup.hs @@ -0,0 +1,2 @@ +import Distribution.Simple +main = defaultMain diff --git a/dph-examples.cabal b/dph-examples.cabal new file mode 100644 --- /dev/null +++ b/dph-examples.cabal @@ -0,0 +1,240 @@ + + + +Name: dph-examples +Version: 0.5.1.1 +License: BSD3 +License-file: LICENSE +Author: The DPH Team +Maintainer: Ben Lippmeier +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 diff --git a/imaginary/DotProduct/dph/DotProductVector.hs b/imaginary/DotProduct/dph/DotProductVector.hs new file mode 100644 --- /dev/null +++ b/imaginary/DotProduct/dph/DotProductVector.hs @@ -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 diff --git a/imaginary/DotProduct/dph/DotProductVectorised.hs b/imaginary/DotProduct/dph/DotProductVectorised.hs new file mode 100644 --- /dev/null +++ b/imaginary/DotProduct/dph/DotProductVectorised.hs @@ -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) diff --git a/imaginary/DotProduct/dph/Main.hs b/imaginary/DotProduct/dph/Main.hs new file mode 100644 --- /dev/null +++ b/imaginary/DotProduct/dph/Main.hs @@ -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 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 + diff --git a/imaginary/Evens/dph/EvensVector.hs b/imaginary/Evens/dph/EvensVector.hs new file mode 100644 --- /dev/null +++ b/imaginary/Evens/dph/EvensVector.hs @@ -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 + diff --git a/imaginary/Evens/dph/EvensVectorised.hs b/imaginary/Evens/dph/EvensVectorised.hs new file mode 100644 --- /dev/null +++ b/imaginary/Evens/dph/EvensVectorised.hs @@ -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)) + diff --git a/imaginary/Evens/dph/Main.hs b/imaginary/Evens/dph/Main.hs new file mode 100644 --- /dev/null +++ b/imaginary/Evens/dph/Main.hs @@ -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 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' + diff --git a/imaginary/Reverse/dph/Main.hs b/imaginary/Reverse/dph/Main.hs new file mode 100644 --- /dev/null +++ b/imaginary/Reverse/dph/Main.hs @@ -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 \n" + , " algorithm one of " ++ show ["vectorised"] + , ""] diff --git a/imaginary/Reverse/dph/ReverseVector.hs b/imaginary/Reverse/dph/ReverseVector.hs new file mode 100644 --- /dev/null +++ b/imaginary/Reverse/dph/ReverseVector.hs @@ -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 + diff --git a/imaginary/Reverse/dph/ReverseVectorised.hs b/imaginary/Reverse/dph/ReverseVectorised.hs new file mode 100644 --- /dev/null +++ b/imaginary/Reverse/dph/ReverseVectorised.hs @@ -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 :]) diff --git a/imaginary/SumSquares/dph/Main.hs b/imaginary/SumSquares/dph/Main.hs new file mode 100644 --- /dev/null +++ b/imaginary/SumSquares/dph/Main.hs @@ -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 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 + diff --git a/imaginary/SumSquares/dph/SumSquaresVector.hs b/imaginary/SumSquares/dph/SumSquaresVector.hs new file mode 100644 --- /dev/null +++ b/imaginary/SumSquares/dph/SumSquaresVector.hs @@ -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 + diff --git a/imaginary/SumSquares/dph/SumSquaresVectorised.hs b/imaginary/SumSquares/dph/SumSquaresVectorised.hs new file mode 100644 --- /dev/null +++ b/imaginary/SumSquares/dph/SumSquaresVectorised.hs @@ -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))) diff --git a/lib/Points2D/Generate.hs b/lib/Points2D/Generate.hs new file mode 100644 --- /dev/null +++ b/lib/Points2D/Generate.hs @@ -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 + diff --git a/lib/Points2D/Types.hs b/lib/Points2D/Types.hs new file mode 100644 --- /dev/null +++ b/lib/Points2D/Types.hs @@ -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)) + + diff --git a/lib/Randomish.hs b/lib/Randomish.hs new file mode 100644 --- /dev/null +++ b/lib/Randomish.hs @@ -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 diff --git a/lib/Timing.hs b/lib/Timing.hs new file mode 100644 --- /dev/null +++ b/lib/Timing.hs @@ -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) diff --git a/lib/Util.hs b/lib/Util.hs new file mode 100644 --- /dev/null +++ b/lib/Util.hs @@ -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 + diff --git a/real/NBody/Batch/Config.hs b/real/NBody/Batch/Config.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Batch/Config.hs @@ -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 } diff --git a/real/NBody/Batch/MainArgs.hs b/real/NBody/Batch/MainArgs.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Batch/MainArgs.hs @@ -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" } + + ] diff --git a/real/NBody/Common/Body.hs b/real/NBody/Common/Body.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Common/Body.hs @@ -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) diff --git a/real/NBody/Common/Dump.hs b/real/NBody/Common/Dump.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Common/Dump.hs @@ -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 diff --git a/real/NBody/Common/Util.hs b/real/NBody/Common/Util.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Common/Util.hs @@ -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 diff --git a/real/NBody/Common/World.hs b/real/NBody/Common/World.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Common/World.hs @@ -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' } + + diff --git a/real/NBody/Gloss/Config.hs b/real/NBody/Gloss/Config.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Gloss/Config.hs @@ -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 " + 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 } diff --git a/real/NBody/Gloss/Draw.hs b/real/NBody/Gloss/Draw.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Gloss/Draw.hs @@ -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 diff --git a/real/NBody/Gloss/Main.hs b/real/NBody/Gloss/Main.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Gloss/Main.hs @@ -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 + + diff --git a/real/NBody/Gloss/MainArgs.hs b/real/NBody/Gloss/MainArgs.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Gloss/MainArgs.hs @@ -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" } + + ] diff --git a/real/NBody/Solver.hs b/real/NBody/Solver.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Solver.hs @@ -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' diff --git a/real/NBody/Solver/ListBH/Solver.hs b/real/NBody/Solver/ListBH/Solver.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Solver/ListBH/Solver.hs @@ -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 + diff --git a/real/NBody/Solver/NestedBH/Solver.hs b/real/NBody/Solver/NestedBH/Solver.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Solver/NestedBH/Solver.hs @@ -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 + diff --git a/real/NBody/Solver/VectorBH/Solver.hs b/real/NBody/Solver/VectorBH/Solver.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Solver/VectorBH/Solver.hs @@ -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 + + diff --git a/real/NBody/Solver/VectorNaive/Solver.hs b/real/NBody/Solver/VectorNaive/Solver.hs new file mode 100644 --- /dev/null +++ b/real/NBody/Solver/VectorNaive/Solver.hs @@ -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) diff --git a/spectral/QuickHull/dph/Main.hs b/spectral/QuickHull/dph/Main.hs new file mode 100644 --- /dev/null +++ b/spectral/QuickHull/dph/Main.hs @@ -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 [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 diff --git a/spectral/QuickHull/dph/QuickHullVect.hs b/spectral/QuickHull/dph/QuickHullVect.hs new file mode 100644 --- /dev/null +++ b/spectral/QuickHull/dph/QuickHullVect.hs @@ -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)) + diff --git a/spectral/QuickHull/lib/SVG.hs b/spectral/QuickHull/lib/SVG.hs new file mode 100644 --- /dev/null +++ b/spectral/QuickHull/lib/SVG.hs @@ -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 + $ [ "" ] + ++ [svgPolygon hull] + ++ map svgPoint points + ++ map svgPointHull hull + ++ [""] + +svgPolygon :: [(Int, Int)] -> String +svgPolygon points + = "" + +svgPoint :: (Int, Int) -> String +svgPoint (x, y) + = "" + +svgPointHull :: (Int, Int) -> String +svgPointHull (x, y) + = "" + +roundPoints :: [(Double, Double)] -> [(Int, Int)] +roundPoints ps = [(round x, round y) | (x, y) <- ps] diff --git a/spectral/QuickHull/vector/Main.hs b/spectral/QuickHull/vector/Main.hs new file mode 100644 --- /dev/null +++ b/spectral/QuickHull/vector/Main.hs @@ -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 [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 + diff --git a/spectral/QuickHull/vector/QuickHullIO.hs b/spectral/QuickHull/vector/QuickHullIO.hs new file mode 100644 --- /dev/null +++ b/spectral/QuickHull/vector/QuickHullIO.hs @@ -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) + diff --git a/spectral/QuickHull/vector/QuickHullSplit.hs b/spectral/QuickHull/vector/QuickHullSplit.hs new file mode 100644 --- /dev/null +++ b/spectral/QuickHull/vector/QuickHullSplit.hs @@ -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) diff --git a/spectral/QuickHull/vector/QuickHullVector.hs b/spectral/QuickHull/vector/QuickHullVector.hs new file mode 100644 --- /dev/null +++ b/spectral/QuickHull/vector/QuickHullVector.hs @@ -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) + diff --git a/spectral/QuickSort/dph/Main.hs b/spectral/QuickSort/dph/Main.hs new file mode 100644 --- /dev/null +++ b/spectral/QuickSort/dph/Main.hs @@ -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 " ] + + +-- | 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) diff --git a/spectral/QuickSort/dph/QuickSortVect.hs b/spectral/QuickSort/dph/QuickSortVect.hs new file mode 100644 --- /dev/null +++ b/spectral/QuickSort/dph/QuickSortVect.hs @@ -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) diff --git a/spectral/SMVM/dph/Main.hs b/spectral/SMVM/dph/Main.hs new file mode 100644 --- /dev/null +++ b/spectral/SMVM/dph/Main.hs @@ -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 " ] + +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) diff --git a/spectral/SMVM/dph/SMVMVectorised.hs b/spectral/SMVM/dph/SMVMVectorised.hs new file mode 100644 --- /dev/null +++ b/spectral/SMVM/dph/SMVMVectorised.hs @@ -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 :] +