diff --git a/LICENSE b/LICENSE --- a/LICENSE +++ b/LICENSE @@ -1,5 +1,4 @@ -Copyright (c) 2001-2011, The DPH Team -All rights reserved. +Copyright (c) 2001-2012, The DPH Team The DPH Team is: Manuel M T Chakravarty @@ -10,22 +9,28 @@ 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. +- 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 name of the University 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 UNIVERSITY COURT OF THE UNIVERSITY OF +GLASGOW AND THE CONTRIBUTORS "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 THE +UNIVERSITY COURT OF THE UNIVERSITY OF GLASGOW OR THE CONTRIBUTORS 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/dph-examples.cabal b/dph-examples.cabal --- a/dph-examples.cabal +++ b/dph-examples.cabal @@ -1,8 +1,7 @@ - Name: dph-examples -Version: 0.5.1.2 +Version: 0.6.0.1 License: BSD3 License-file: LICENSE Author: The DPH Team @@ -12,215 +11,154 @@ Stability: experimental Category: Data Structures Homepage: http://www.haskell.org/haskellwiki/GHC/Data_Parallel_Haskell -Description: Examples using Data Parallel Haskell -Synopsis: Examples using Data Parallel Haskell - --- 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 +Description: Data Parallel Haskell example programs. +Synopsis: Data Parallel Haskell example programs. --- 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.* +-- Smoke ---------------------------------------------------------------------- +-- examples/smoke/data +Executable dph-smoke-bool + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* 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 + other-modules: Vectorised + hs-source-dirs: examples/smoke/data/Bool + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 -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 +-- examples/smoke/prims +Executable dph-smoke-concat + build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* + main-is: Main.hs + other-modules: Vectorised + hs-source-dirs: examples/smoke/prims/Concat + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 --- 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.* + +Executable dph-smoke-sumsq + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* Main-is: Main.hs - other-modules: DotProductVector - DotProductVectorised + other-modules: Vector + Vectorised 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.* + hs-source-dirs: examples/smoke/prims/SumSquares/dph lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 + + +Executable dph-smoke-evens + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* Main-is: Main.hs - other-modules: DotProductVector - DotProductVectorised + other-modules: Vector + Vectorised 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 + hs-source-dirs: examples/smoke/prims/Evens/dph lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 --- 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.* +-- examples/smoke/sharing +Executable dph-smoke-indices + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* 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.* + other-modules: Vectorised + hs-source-dirs: examples/smoke/sharing/Indices lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 + + +Executable dph-smoke-rank + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* 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 + other-modules: Vectorised Util Timing Randomish + hs-source-dirs: examples/smoke/sharing/Rank lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 --- 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-smoke-reverse + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* + Main-is: Main.hs + other-modules: Vectorised Randomish + hs-source-dirs: examples/smoke/sharing/Reverse lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 --- 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 +-- Imaginary ------------------------------------------------------------------ +-- Executable dph-imaginary-primes +-- Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* +-- Main-is: Main.hs +-- other-modules: Vectorised +-- hs-source-dirs: examples/imaginary/Primes lib +-- ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 --- 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.* +Executable dph-imaginary-words + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* 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 + other-modules: Vectorised + hs-source-dirs: examples/imaginary/Words lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 -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.* +-- Spectral ------------------------------------------------------------------- +Executable dph-spectral-dotp + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* 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 + other-modules: Vector + Vectorised + Timing Randomish + hs-source-dirs: examples/spectral/DotProduct/dph lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 --- 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.* +Executable dph-spectral-smvm + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* Main-is: Main.hs - other-modules: SMVMVectorised + other-modules: Vectorised Timing - hs-source-dirs: spectral/SMVM/dph lib - ghc-options: -rtsopts -threaded -fllvm -Odph -fdph-seq -fcpr-off -fsimplifier-phases=4 -fstrictness-before=3 - + hs-source-dirs: examples/spectral/SMVM/dph lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 --- 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.* +Executable dph-spectral-quickhull + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* Main-is: Main.hs - other-modules: QuickHullVect + other-modules: Vectorised 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 + hs-source-dirs: examples/spectral/QuickHull/dph examples/spectral/QuickHull/lib lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 -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.* +Executable dph-spectral-quickhull-vector + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* 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 + hs-source-dirs: examples/spectral/QuickHull/vector examples/spectral/QuickHull/lib lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 --- 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.* +Executable dph-spectral-quicksort + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* Main-is: Main.hs - other-modules: QuickSortVect + other-modules: Vectorised 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 + hs-source-dirs: examples/spectral/QuickSort/dph lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 --- NBody ---------------------------------------------------------------------- -Executable dph-nbody-gloss +-- Real ----------------------------------------------------------------------- +Executable dph-real-nbody-gloss Main-is: MainGloss.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 - 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 - + Solver Solver.ListBH.Solver + Solver.NestedBH.Solver + Solver.VectorBH.Solver + Solver.VectorNaive.Solver + Timing Points2D.Types Points2D.Generate + System.Console.ParseArgs + Gloss.MainArgs Gloss.Draw Gloss.Config + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.*, gloss == 1.6.* + hs-source-dirs: examples/real/NBody examples/real/NBody/Gloss lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 -Executable dph-nbody-batch +Executable dph-real-nbody Main-is: MainBatch.hs other-modules: Common.Dump Common.World Common.Body Common.Util Solver Solver.ListBH.Solver @@ -229,18 +167,7 @@ Solver.VectorNaive.Solver Timing Points2D.Types Points2D.Generate 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 + Build-depends: base == 4.5.*, vector == 0.9.*, random == 1.0.*, old-time == 1.1.*, containers == 0.4.*, dph-base == 0.6.*, dph-prim-par == 0.6.*, dph-lifted-vseg == 0.6.*, HUnit == 1.2.* + hs-source-dirs: examples/real/NBody examples/real/NBody/Batch lib + ghc-options: -rtsopts -threaded -fllvm -Odph -package dph-lifted-vseg -fcpr-off -fno-liberate-case -fsimpl-tick-factor=1000 --- 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/examples/imaginary/Words/Main.hs b/examples/imaginary/Words/Main.hs new file mode 100644 --- /dev/null +++ b/examples/imaginary/Words/Main.hs @@ -0,0 +1,37 @@ + +import Vectorised +import Data.Array.Parallel +import qualified Data.Array.Parallel.Prelude.Word8 as W +import qualified Data.Array.Parallel.PArray as P +import qualified Data.Array.Parallel.Unlifted as U +import Data.Char + +main + = do -- take the filename containing the words as the first arg + let str = "When I look into the looking glass I'm always sure to see" + ++ " no matter how I dodge about, me looking back at me." + + -- convert string to a PArray + let paStr :: PArray W.Word8 + paStr = P.fromUArray $ U.map W.fromInt $ U.fromList $ map ord str + + + -- break the string into words then flatten it back + let str' :: String + str' = map chr + $ map fromIntegral + $ P.toList + $ wordsOfPArray paStr + + + -- count the number of words in the string, using the vectorised program + let wordCountVect = fromIntegral $ wordCountOfPArray paStr + + -- count the number of words with the ye'olde list way + let wordCountList = length $ words str + + -- + putStr $ show str' ++ "\n" + ++ "word count vect = " ++ show wordCountVect ++ "\n" + ++ "word count lists = " ++ show wordCountList ++ "\n" + diff --git a/examples/imaginary/Words/Vectorised.hs b/examples/imaginary/Words/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/imaginary/Words/Vectorised.hs @@ -0,0 +1,113 @@ +{-# LANGUAGE ParallelArrays, ParallelListComp #-} +{-# OPTIONS -fvectorise #-} + +module Vectorised + ( wordsOfPArray + , wordCountOfPArray ) +where +import qualified Data.Array.Parallel.Prelude.Word8 as W +import Data.Array.Parallel.Prelude.Word8 (Word8) +import Data.Array.Parallel.Prelude.Int as I +import Data.Array.Parallel + +import qualified Prelude as Prel + + +-- We can't use the Prelude Char and String types in vectorised code yet.. +type Char = Word8 +char_space = W.fromInt 32 + +type String = [: Char :] + + +-- | Word state +data State + = Chunk String + | Seg String -- initial word chunk + [:String:] -- complete words in the middle of the segment + String -- final word chunk + + +-- | Compose two wordstates. +plusState :: State -> State -> State +plusState str1 str2 + = case (str1, str2) of + (Chunk as, Chunk bs) -> Chunk (as +:+ bs) + (Chunk as, Seg bl bss br) -> Seg (as +:+ bl) bss br + (Seg al ass ar, Chunk bs) -> Seg al ass (ar +:+ bs) + (Seg al ass ar, Seg bl bss br) -> Seg al (ass +:+ joinEmpty [:ar +:+ bl:] +:+ bss) br + +joinEmpty :: [:[:Word8:]:] -> [:[:Word8:]:] +joinEmpty ws + | lengthP ws I.== 1 && lengthP (ws !: 0) I.== 0 = [::] + | otherwise = ws + + +-- | Convert a single char to a wordstate. +stateOfChar :: Char -> State +stateOfChar c + | c W.== char_space = Seg [::] [::] [::] + | otherwise = Chunk [:c:] + + +-- | Break this string into words. +stateOfString :: String -> State +stateOfString str + = let len = lengthP str + result + | len I.== 0 = Chunk [::] + | len I.== 1 = stateOfChar (str !: 0) + | otherwise + = let half = len `div` 2 + s1 = sliceP 0 half str + s2 = sliceP half (len I.- half) str + in plusState (stateOfString s1) (stateOfString s2) + in result + + +-- | Count the number of words in a string. +countWordsOfState :: State -> Int +countWordsOfState state + = case state of + Chunk c -> wordsInChunkArr c + Seg c1 ws c2 -> wordsInChunkArr c1 I.+ lengthP ws I.+ wordsInChunkArr c2 + +wordsInChunkArr :: [:Word8:] -> Int +wordsInChunkArr arr + | lengthP arr I.== 0 = 0 + | otherwise = 1 + + +-- | Flatten a state back to an array of Word8s, +-- inserting spaces between the words. +flattenState :: State -> [:Word8:] +flattenState ss + = case ss of + Chunk s -> s + + Seg w1 ws w2 + -> w1 + +:+ [:char_space:] + +:+ concatP [: w +:+ [:char_space:] | w <- ws :] + +:+ w2 + +-- Interface ------------------------------------------------------------------ + +-- | Break up an array of chars into words then flatten it back. +{-# NOINLINE wordsOfPArray #-} +wordsOfPArray :: PArray Word8 -> PArray Word8 +wordsOfPArray arr + = let str = fromPArrayP arr + state = stateOfString str + strOut = flattenState state + in toPArrayP strOut + + +-- | Count the number of words in an array +{-# NOINLINE wordCountOfPArray #-} +wordCountOfPArray :: PArray Word8 -> Int +wordCountOfPArray arr + = let str = fromPArrayP arr + state = stateOfString str + in countWordsOfState state + diff --git a/examples/real/NBody/Batch/Config.hs b/examples/real/NBody/Batch/Config.hs new file mode 100644 --- /dev/null +++ b/examples/real/NBody/Batch/Config.hs @@ -0,0 +1,65 @@ + +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 + + -- Print points to stdout + , configPrintFinal :: Bool + + -- Print timings to stdout + , configPrintTimings :: Bool } + + +-- | 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 + , configPrintFinal = gotArg args ArgPrintFinal + , configPrintTimings = gotArg args ArgPrintTimings } diff --git a/examples/real/NBody/Batch/MainArgs.hs b/examples/real/NBody/Batch/MainArgs.hs new file mode 100644 --- /dev/null +++ b/examples/real/NBody/Batch/MainArgs.hs @@ -0,0 +1,99 @@ + +module Batch.MainArgs + ( MainArg(..) + , mainArgs) +where +import System.Console.ParseArgs + +data MainArg + = ArgHelp + | ArgSolver + | ArgMaxSteps + | ArgTimeStep + | ArgBodyCount + | ArgBodyMass + | ArgEpsilon + | ArgDiscSize + | ArgStartSpeed + | ArgDumpFinal + | ArgPrintFinal + | ArgPrintTimings + 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-bh" + , argDesc = "One of: list-bh, vector-naive, vector-bh, nested-bh. (default vector-bh)" } + + -- 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" } + + , Arg { argIndex = ArgPrintFinal + , argAbbr = Nothing + , argName = Just "print-final" + , argData = Nothing + , argDesc = "Print final body positions and masses to stdout" } + + , Arg { argIndex = ArgPrintTimings + , argAbbr = Nothing + , argName = Just "print-timings" + , argData = Nothing + , argDesc = "Print timings of the solver" } + + ] diff --git a/examples/real/NBody/Common/Body.hs b/examples/real/NBody/Common/Body.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/real/NBody/Common/Dump.hs b/examples/real/NBody/Common/Dump.hs new file mode 100644 --- /dev/null +++ b/examples/real/NBody/Common/Dump.hs @@ -0,0 +1,46 @@ +{-# LANGUAGE ScopedTypeVariables #-} +module Common.Dump + ( dumpWorld + , printWorld) +where +import Common.Body +import Common.World +import System.IO +import Data.List +import Text.Printf +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 + hClose h + + +-- | Print the bodies in a world to stdout +printWorld :: World -> IO () +printWorld world + = do mapM_ (hWriteBody stdout) + $ 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 (\f -> printf "%15.8f" f) [ 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/examples/real/NBody/Common/Util.hs b/examples/real/NBody/Common/Util.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/real/NBody/Common/World.hs b/examples/real/NBody/Common/World.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/real/NBody/Gloss/Config.hs b/examples/real/NBody/Gloss/Config.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/real/NBody/Gloss/Draw.hs b/examples/real/NBody/Gloss/Draw.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/real/NBody/Gloss/MainArgs.hs b/examples/real/NBody/Gloss/MainArgs.hs new file mode 100644 --- /dev/null +++ b/examples/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-bh" + , argDesc = "One of: list-bh, vector-naive, vector-bh, nested-bh. (default vector-bh)" } + + -- 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/examples/real/NBody/MainBatch.hs b/examples/real/NBody/MainBatch.hs new file mode 100644 --- /dev/null +++ b/examples/real/NBody/MainBatch.hs @@ -0,0 +1,102 @@ +{-# LANGUAGE ParallelListComp, BangPatterns #-} + +import Batch.MainArgs +import Batch.Config +import Common.Dump +import Common.World +import Common.Body +import Common.Util +import Solver +import Timing +import Points2D.Generate +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 mainBatch config calcAccels worldStart + + +-- | Run the simulation in batch mode. +mainBatch :: Config -> Solver -> World -> IO () +mainBatch config calcAccels worldStart + = do + worldStart `seq` return () + + (world', tElapsed) + <- time + $ let world = mainBatchRun config calcAccels worldStart + in world `seq` return world + + when (configPrintTimings config) + $ putStr $ prettyTime tElapsed + + mainEnd (configDumpFinal config) + (configPrintFinal config) + world' + + +mainBatchRun config calcAccels worldStart + = go worldStart + where go !world + = let world' = advanceWorld + (calcAccels $ configEpsilon config) + (configTimeStep config) + world + + in if worldSteps world' < configMaxSteps config + then go world' + else world' + + +-- | Called at end of run to dump final world state. +mainEnd :: Maybe FilePath -- ^ Write final bodies to this file. + -> Bool -- ^ Print final bodies to stdout + -> World -- ^ Final world state. + -> IO () + +mainEnd mDumpFinal printFinal world + = do -- Dump the final world state to file if requested. + maybe (return ()) (dumpWorld world) mDumpFinal + when printFinal (printWorld world) + + diff --git a/examples/real/NBody/MainGloss.hs b/examples/real/NBody/MainGloss.hs new file mode 100644 --- /dev/null +++ b/examples/real/NBody/MainGloss.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 simulate + (InWindow "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/examples/real/NBody/Solver.hs b/examples/real/NBody/Solver.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/real/NBody/Solver/ListBH/Solver.hs b/examples/real/NBody/Solver/ListBH/Solver.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/real/NBody/Solver/NestedBH/Solver.hs b/examples/real/NBody/Solver/NestedBH/Solver.hs new file mode 100644 --- /dev/null +++ b/examples/real/NBody/Solver/NestedBH/Solver.hs @@ -0,0 +1,156 @@ +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +module Solver.NestedBH.Solver + (calcAccelsWithBoxPA) +where +import Data.Array.Parallel +import Data.Array.Parallel.Prelude.Bool +import Data.Array.Parallel.Prelude.Double as D +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' +{-# NOINLINE calcAccelsWithBoxPA #-} + +-- | 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 D.- llx + sy = ruy D.- lly + s = if sx D.< 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 D.+ rux) D./ 2.0, (lly D.+ ruy) D./ 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 D.> llx) && (px D.<= rux) && (py D.> lly) && (py D.<= ruy) + + +-- | Calculate the centroid of some points. +calcCentroid:: [:MassPoint:] -> MassPoint +calcCentroid mpts + = MP (sumP xs D./ mass) (sumP ys D./ mass) mass + where mass = sumP [: m | MP _ _ m <- mpts :] + (xs, ys) = unzipP [: (m D.* x, m D.* 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 D.* dx D./ r , aabs D.* dy D./ r) + where rsqr = (dx D.* dx) D.+ (dy D.* dy) D.+ epsilon D.* epsilon + r = sqrt rsqr + dx = x1 D.- x2 + dy = y1 D.- y2 + aabs = m D./ 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 D.- x1 + dy = y2 D.- y1 + dist = sqrt (dx D.* dx D.+ dy D.* dy) + in (s D./ dist) D.< 1 + diff --git a/examples/real/NBody/Solver/VectorBH/Solver.hs b/examples/real/NBody/Solver/VectorBH/Solver.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/real/NBody/Solver/VectorNaive/Solver.hs b/examples/real/NBody/Solver/VectorNaive/Solver.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/smoke/data/Bool/Main.hs b/examples/smoke/data/Bool/Main.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/data/Bool/Main.hs @@ -0,0 +1,7 @@ + +-- | Test vectorisation of enumerations. +-- This tests the conversion to and from our generic representation. +import Vectorised +import qualified Data.Array.Parallel.PArray as P + +main = print $ P.toList $ test $ P.fromList [0, 1, 1, 0, 1, 0, 1, 1, 0, 0, 1] diff --git a/examples/smoke/data/Bool/Vectorised.hs b/examples/smoke/data/Bool/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/data/Bool/Vectorised.hs @@ -0,0 +1,28 @@ + +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +module Vectorised (test) where +import Data.Array.Parallel hiding (Bool, True, False, not, fromBool, toBool) +import Data.Array.Parallel.Prelude hiding (Bool, True, False, not, fromBool, toBool) +import Data.Array.Parallel.Prelude.Int as I +import qualified Prelude as P + +data Bool = True | False + +toBool :: Int -> Bool +toBool n + | n I.== 0 = False + | otherwise = True + +fromBool :: Bool -> Int +fromBool False = 0 +fromBool True = 1 + +not :: Bool -> Bool +not False = True +not True = False + +{-# NOINLINE test #-} +test :: PArray Int -> PArray Int +test xs = toPArrayP (mapP test' (fromPArrayP xs)) +test' x = fromBool (not (toBool x)) diff --git a/examples/smoke/prims/Concat/Main.hs b/examples/smoke/prims/Concat/Main.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/prims/Concat/Main.hs @@ -0,0 +1,23 @@ + +-- | Test out contatentation. +-- Concatenation is a core operation that exercises the representation of nested +-- arrays, as well as the extractsPR function from the dph-lifted-vseg library. +import Vectorised +import Test.HUnit +import qualified Data.Array.Parallel.PArray as PA + +arr2 :: [[Int]] +arr2 = [ [0, 1, 2, 3], [4, 5], [], [6, 7, 8], [9]] + +arr3 :: [[[Int]]] +arr3 = [[[0]], [[1], [2, 3]], [[4, 5]], [], [[6, 7], [8]], [[]], [[9]] ] + + +main = runTestTT $ test $ + [ "test0" ~: PA.toList test0 ~?= concat arr2 + , "test1" ~: PA.toList test1 ~?= concat (concat arr3) + , "test2" ~: PA.toList test2 ~?= concat (map concat arr3) + , "test3" ~: PA.toList test3 ~?= concat [ map (+1) x | x <- [ [1, 2], [3, 4] ] ] + , "test4" ~: PA.toList test4 ~?= concat (map (map (+1)) arr2) + ] + diff --git a/examples/smoke/prims/Concat/Vectorised.hs b/examples/smoke/prims/Concat/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/prims/Concat/Vectorised.hs @@ -0,0 +1,38 @@ + +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +module Vectorised + (test0, test1, test2, test3, test4) +where +import Data.Array.Parallel +import Data.Array.Parallel.Prelude.Int as I +import qualified Prelude as P + +arr2 :: [:[:Int:]:] +arr2 = [: [:0, 1, 2, 3:], [:4, 5:], [::], [:6, 7, 8:], [:9:]:] + +arr3 :: [:[:[:Int:]:]:] +arr3 = [: [:[:0:]:], [:[:1:], [:2, 3:]:], [:[:4, 5:]:], [::], [:[:6, 7:], [:8:]:], [:[::]:], [:[:9:]:] :] + + +test0' :: [:Int:] +test0' = concatP arr2 +test0 = toPArrayP test0' + +test1' :: [:Int:] +test1' = concatP (concatP arr3) +test1 = toPArrayP test1' + +test2' :: [:Int:] +test2' = concatP (mapP concatP arr3) +test2 = toPArrayP test2' + + +test3' :: [:Int:] +test3' = concatP [: mapP (I.+ 1) x | x <- [: [:1, 2:], [:3, 4:] :] :] +test3 = toPArrayP test3' + + +test4' :: [:Int:] +test4' = concatP (mapP (mapP (I.+ 1)) arr2) +test4 = toPArrayP test4' diff --git a/examples/smoke/prims/Evens/dph/Main.hs b/examples/smoke/prims/Evens/dph/Main.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/prims/Evens/dph/Main.hs @@ -0,0 +1,48 @@ + +-- | Test filter operations. +-- These are fairly involved because the implementation uses selectors and +-- the packByTag array primitive. +import Vectorised as Z +import Vector 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.fromUArray 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/examples/smoke/prims/Evens/dph/Vector.hs b/examples/smoke/prims/Evens/dph/Vector.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/prims/Evens/dph/Vector.hs @@ -0,0 +1,8 @@ + +module Vector (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/examples/smoke/prims/Evens/dph/Vectorised.hs b/examples/smoke/prims/Evens/dph/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/prims/Evens/dph/Vectorised.hs @@ -0,0 +1,16 @@ + +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +module Vectorised (evensPA) where +import Data.Array.Parallel +import Data.Array.Parallel.Prelude.Int as I +import Data.Array.Parallel.Prelude.Bool +import qualified Prelude as P + +evens :: [:Int:] -> [:Int:] +evens ints = filterP (\x -> x `mod` 2 I.== 0) ints + +evensPA :: PArray Int -> PArray Int +{-# NOINLINE evensPA #-} +evensPA arr = toPArrayP (evens (fromPArrayP arr)) + diff --git a/examples/smoke/prims/SumSquares/dph/Main.hs b/examples/smoke/prims/SumSquares/dph/Main.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/prims/SumSquares/dph/Main.hs @@ -0,0 +1,35 @@ + +-- | Test maps, sums and enumerations. +import Timing +import Randomish +import System.Environment +import qualified Vector as V +import qualified Vectorised as Z +import qualified Data.Vector.Unboxed as V + +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/examples/smoke/prims/SumSquares/dph/Vector.hs b/examples/smoke/prims/SumSquares/dph/Vector.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/prims/SumSquares/dph/Vector.hs @@ -0,0 +1,11 @@ + +module Vector 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/examples/smoke/prims/SumSquares/dph/Vectorised.hs b/examples/smoke/prims/SumSquares/dph/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/prims/SumSquares/dph/Vectorised.hs @@ -0,0 +1,16 @@ + +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +module Vectorised (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/examples/smoke/sharing/Indices/Main.hs b/examples/smoke/sharing/Indices/Main.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/sharing/Indices/Main.hs @@ -0,0 +1,42 @@ + +-- | This divide and conquer program accesseses top-level array from the +-- computations at the bottom of the tree. In lifted backends that +-- don't manage sharing properly, this program will blow up when it tries +-- to replicate the top-level a array at every step in the division phase. +import Util +import Timing +import Randomish +import System.Environment +import Control.Exception +import qualified Vectorised as ID +import qualified Data.Array.Parallel.PArray as P +import qualified Data.Vector.Unboxed as V + +main + = do args <- getArgs + + case args of + [alg, count] -> run alg (read count) + _ -> usage + + +run "vectorised" count + = do let arr = P.fromList [0 .. count - 1] + arr `seq` return () + + (arrResult, tElapsed) + <- time + $ let arr' = ID.indicesPA arr arr + in P.nf arr' `seq` return arr' + + print $ P.length arrResult + putStr $ prettyTime tElapsed + +run _ _ + = usage + + +usage = putStr $ unlines + [ "usage: indices \n" + , " algorithm one of " ++ show ["vectorised"] + , ""] diff --git a/examples/smoke/sharing/Indices/Vectorised.hs b/examples/smoke/sharing/Indices/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/sharing/Indices/Vectorised.hs @@ -0,0 +1,31 @@ + +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +module Vectorised (indicesPA, indices) where +import Data.Array.Parallel +import Data.Array.Parallel.Prelude.Int as I +import Data.Array.Parallel.Prelude.Bool +import qualified Prelude as P + + +{-# NOINLINE indicesPA #-} +indicesPA :: PArray Int -> PArray Int -> PArray Int +indicesPA arr ixs + = toPArrayP (indices (fromPArrayP arr) (fromPArrayP ixs)) + +indices :: [:Int:] -> [:Int:] -> [:Int:] +indices arr ixs + = treeLookup arr ixs + +{-# NOINLINE treeLookup #-} +treeLookup :: [:Int:] -> [:Int:] -> [:Int:] +treeLookup table xx + | lengthP xx I.== 1 + = [: table !: (xx !: 0) :] + + | otherwise + = let len = lengthP xx + half = len `div` 2 + s1 = sliceP 0 half xx + s2 = sliceP half half xx + in concatP (mapP (treeLookup table) [: s1, s2 :]) diff --git a/examples/smoke/sharing/Rank/Main.hs b/examples/smoke/sharing/Rank/Main.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/sharing/Rank/Main.hs @@ -0,0 +1,40 @@ + +-- | In with lifted backends that don't manage sharing properly, the rank +-- example suffers similar work complexity problems as the indices example. +import Util +import Timing +import Randomish +import System.Environment +import Control.Exception +import qualified Vectorised as RD +import qualified Data.Array.Parallel.PArray as PA +import qualified Data.Vector.Unboxed as V + +main + = do args <- getArgs + + case args of + [alg, count] -> run alg (read count) + _ -> usage + + +run "vectorised" count + = do let arr = PA.fromList [0 .. count - 1] + arr `seq` return () + + (arrRanks, tElapsed) + <- time + $ let arr' = RD.ranksPA arr + in PA.nf arr' `seq` return arr' + + print $ PA.length arrRanks + putStr $ prettyTime tElapsed + +run _ _ + = usage + + +usage = putStr $ unlines + [ "usage: rank \n" + , " algorithm one of " ++ show ["vectorised"] + , ""] diff --git a/examples/smoke/sharing/Rank/Vectorised.hs b/examples/smoke/sharing/Rank/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/sharing/Rank/Vectorised.hs @@ -0,0 +1,18 @@ +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +module Vectorised (ranksPA) where +import Data.Array.Parallel +import Data.Array.Parallel.Prelude.Int as I +import Data.Array.Parallel.Prelude.Bool +import qualified Prelude as P + + +ranksPA :: PArray Int -> PArray Int +ranksPA ps = toPArrayP (ranks (fromPArrayP ps)) +{-# NOINLINE ranksPA #-} + + +ranks :: [:Int:] -> [:Int:] +ranks arr = [: lengthP [: a | a <- arr, a I.< b :] | b <- arr :] +{-# NOINLINE ranks #-} + diff --git a/examples/smoke/sharing/Reverse/Main.hs b/examples/smoke/sharing/Reverse/Main.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/sharing/Reverse/Main.hs @@ -0,0 +1,58 @@ + +import Util +import Timing +import Randomish +import System.Environment +import System.Random +import Control.Exception +import qualified Vector as RV +import qualified Vectorised 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/examples/smoke/sharing/Reverse/Vectorised.hs b/examples/smoke/sharing/Reverse/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/smoke/sharing/Reverse/Vectorised.hs @@ -0,0 +1,28 @@ +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +module Vectorised (treeReversePA) where +import Data.Array.Parallel +import Data.Array.Parallel.Prelude.Int as I +import Data.Array.Parallel.Prelude.Bool +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 I.== 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/examples/spectral/DotProduct/dph/Main.hs b/examples/spectral/DotProduct/dph/Main.hs new file mode 100644 --- /dev/null +++ b/examples/spectral/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 Vector as V +import qualified Vectorised 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 = P.fromUArray vec1 + let arr2 = P.fromUArray 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/examples/spectral/DotProduct/dph/Vector.hs b/examples/spectral/DotProduct/dph/Vector.hs new file mode 100644 --- /dev/null +++ b/examples/spectral/DotProduct/dph/Vector.hs @@ -0,0 +1,7 @@ + +module Vector (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/examples/spectral/DotProduct/dph/Vectorised.hs b/examples/spectral/DotProduct/dph/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/spectral/DotProduct/dph/Vectorised.hs @@ -0,0 +1,15 @@ +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +module Vectorised ( dotPA ) where + +import Data.Array.Parallel +import Data.Array.Parallel.Prelude.Double as D + +import qualified Prelude + +dotPA :: PArray Double -> PArray Double -> Double +dotPA v w = dotp' (fromPArrayP v) (fromPArrayP w) +{-# NOINLINE dotPA #-} + +dotp' :: [:Double:] -> [:Double:] -> Double +dotp' v w = D.sumP (zipWithP (D.*) v w) diff --git a/examples/spectral/QuickHull/dph/Main.hs b/examples/spectral/QuickHull/dph/Main.hs new file mode 100644 --- /dev/null +++ b/examples/spectral/QuickHull/dph/Main.hs @@ -0,0 +1,59 @@ + +import Vectorised +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/examples/spectral/QuickHull/dph/Vectorised.hs b/examples/spectral/QuickHull/dph/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/spectral/QuickHull/dph/Vectorised.hs @@ -0,0 +1,43 @@ +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} + +module Vectorised (quickhullPA) where +import Points2D.Types +import Data.Array.Parallel +import Data.Array.Parallel.Prelude.Bool +import Data.Array.Parallel.Prelude.Double as D +import qualified Data.Array.Parallel.Prelude.Int as I +import qualified Prelude as P + + +distance :: Point -> Line -> Double +distance (xo, yo) ((x1, y1), (x2, y2)) + = (x1 D.- xo) D.* (y2 D.- yo) D.- (y1 D.- yo) D.* (x2 D.- xo) + + +hsplit :: [:Point:] -> Line -> [:Point:] +hsplit points line@(p1, p2) + | lengthP packed I.== 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 D.> 0.0 :] + pm = points !: maxIndexP cross + + +quickHull :: [:Point:] -> [:Point:] +quickHull points + | lengthP points I.== 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 +quickhullPA ps = toPArrayP (quickHull (fromPArrayP ps)) +{-# NOINLINE quickhullPA #-} + diff --git a/examples/spectral/QuickHull/lib/SVG.hs b/examples/spectral/QuickHull/lib/SVG.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/spectral/QuickHull/vector/Main.hs b/examples/spectral/QuickHull/vector/Main.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/spectral/QuickHull/vector/QuickHullIO.hs b/examples/spectral/QuickHull/vector/QuickHullIO.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/spectral/QuickHull/vector/QuickHullSplit.hs b/examples/spectral/QuickHull/vector/QuickHullSplit.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/spectral/QuickHull/vector/QuickHullVector.hs b/examples/spectral/QuickHull/vector/QuickHullVector.hs new file mode 100644 --- /dev/null +++ b/examples/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/examples/spectral/QuickSort/dph/Main.hs b/examples/spectral/QuickSort/dph/Main.hs new file mode 100644 --- /dev/null +++ b/examples/spectral/QuickSort/dph/Main.hs @@ -0,0 +1,50 @@ + +import Vectorised +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 = fromUArray + $ 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.toUArray 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/examples/spectral/QuickSort/dph/Vectorised.hs b/examples/spectral/QuickSort/dph/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/spectral/QuickSort/dph/Vectorised.hs @@ -0,0 +1,25 @@ +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +{-# OPTIONS -fno-spec-constr-count #-} +module Vectorised (quicksortPA) where +import Data.Array.Parallel +import Data.Array.Parallel.Prelude.Double as D +import qualified Data.Array.Parallel.Prelude.Int as I +import qualified Prelude + + +{-# NOINLINE quicksortPA #-} +quicksortPA:: PArray Double -> PArray Double +quicksortPA xs = toPArrayP (qsortVect' (fromPArrayP xs)) + + +qsortVect':: [: Double :] -> [: Double :] +qsortVect' xs | lengthP xs I.<= 1 = xs + | otherwise = + let p = xs !: (lengthP xs `I.div` 2) + ls = [:x | x <- xs, x D.< p:] + gs = [:x | x <- xs, x D.> p:] + + ss = mapP qsortVect' [:ls, gs:] + in + (ss !: 0) +:+ [:x | x <- xs, x D.== p:] +:+ (ss !: 1) diff --git a/examples/spectral/SMVM/dph/Main.hs b/examples/spectral/SMVM/dph/Main.hs new file mode 100644 --- /dev/null +++ b/examples/spectral/SMVM/dph/Main.hs @@ -0,0 +1,102 @@ +{-# LANGUAGE ScopedTypeVariables #-} + +import Timing +import Vectorised +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 (P.toUArray vector)) +-- putStrLn $ "matrix height = " ++ show (U.length (toUArrPA matrix)) + + + -- Print checksum of resulting vector. + putStrLn $ "result sum = " ++ show (U.sum (P.toUArray 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 = P.nestUSegd segd (P.fromUArray2 arrMatrixElems) + let paVector = P.fromUArray 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/examples/spectral/SMVM/dph/Vectorised.hs b/examples/spectral/SMVM/dph/Vectorised.hs new file mode 100644 --- /dev/null +++ b/examples/spectral/SMVM/dph/Vectorised.hs @@ -0,0 +1,17 @@ + +{-# LANGUAGE ParallelArrays #-} +{-# OPTIONS -fvectorise #-} +module Vectorised (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 +smvmPA m v = toPArrayP (smvm (fromNestedPArrayP m) (fromPArrayP v)) +{-# NOINLINE smvmPA #-} + + +smvm :: [:[: (Int, Double) :]:] -> [:Double:] -> [:Double:] +smvm m v = [: D.sumP [: x D.* (v !: i) | (i,x) <- row :] | row <- m :] + diff --git a/imaginary/DotProduct/dph/DotProductVector.hs b/imaginary/DotProduct/dph/DotProductVector.hs deleted file mode 100644 --- a/imaginary/DotProduct/dph/DotProductVector.hs +++ /dev/null @@ -1,7 +0,0 @@ - -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 deleted file mode 100644 --- a/imaginary/DotProduct/dph/DotProductVectorised.hs +++ /dev/null @@ -1,15 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/imaginary/DotProduct/dph/Main.hs +++ /dev/null @@ -1,45 +0,0 @@ - -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 deleted file mode 100644 --- a/imaginary/Evens/dph/EvensVector.hs +++ /dev/null @@ -1,8 +0,0 @@ - -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 deleted file mode 100644 --- a/imaginary/Evens/dph/EvensVectorised.hs +++ /dev/null @@ -1,16 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/imaginary/Evens/dph/Main.hs +++ /dev/null @@ -1,45 +0,0 @@ - -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 deleted file mode 100644 --- a/imaginary/Reverse/dph/Main.hs +++ /dev/null @@ -1,58 +0,0 @@ - -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 deleted file mode 100644 --- a/imaginary/Reverse/dph/ReverseVector.hs +++ /dev/null @@ -1,21 +0,0 @@ - -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 deleted file mode 100644 --- a/imaginary/Reverse/dph/ReverseVectorised.hs +++ /dev/null @@ -1,29 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/imaginary/SumSquares/dph/Main.hs +++ /dev/null @@ -1,34 +0,0 @@ - -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 deleted file mode 100644 --- a/imaginary/SumSquares/dph/SumSquaresVector.hs +++ /dev/null @@ -1,11 +0,0 @@ - -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 deleted file mode 100644 --- a/imaginary/SumSquares/dph/SumSquaresVectorised.hs +++ /dev/null @@ -1,15 +0,0 @@ -{-# 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 --- a/lib/Points2D/Generate.hs +++ b/lib/Points2D/Generate.hs @@ -91,8 +91,8 @@ pointsPArrayOfUArray ps = do let pts = makePointsPA - (P.fromUArrPA' (U.fsts ps)) - (P.fromUArrPA' (U.snds ps)) + (P.fromUArray (U.fsts ps)) + (P.fromUArray (U.snds ps)) evaluate $ P.nf pts return pts diff --git a/lib/Points2D/Types.hs b/lib/Points2D/Types.hs --- a/lib/Points2D/Types.hs +++ b/lib/Points2D/Types.hs @@ -9,6 +9,9 @@ , ysOf, ysOfPA) where import Data.Array.Parallel +import Data.Array.Parallel.Prelude.Double +import qualified Data.Array.Parallel.Prelude.Double as D +import qualified Prelude as P type Point = (Double, Double) type Line = (Point, Point) diff --git a/lib/System/Console/ParseArgs.hs b/lib/System/Console/ParseArgs.hs new file mode 100644 --- /dev/null +++ b/lib/System/Console/ParseArgs.hs @@ -0,0 +1,667 @@ +{-# LANGUAGE FlexibleInstances, DeriveDataTypeable #-} +-- Full-featured argument parsing library for Haskell programs +-- Bart Massey + +-- Copyright © 2007-2010 Bart Massey +-- ALL RIGHTS RESERVED + +-- You can redistribute and/or modify this library under the +-- terms of the "3-clause BSD LICENSE", as stated in the file +-- COPYING in the top-level directory of this distribution. +-- +-- This library is distributed in the hope that it will be +-- useful, but WITHOUT ANY WARRANTY; without even the +-- implied warranty of MERCHANTABILITY or FITNESS FOR A +-- PARTICULAR PURPOSE. + +-- |This module supplies an argument parser. +-- Given a description of type [`Arg`] of the legal +-- arguments to the program, a list of argument strings, +-- and a bit of extra information, the `parseArgs` function +-- in this module returns an +-- `Args` data structure suitable for querying using the +-- provided functions `gotArg`, `getArg`, etc. +module System.Console.ParseArgs ( + -- * Describing allowed arguments + -- |The argument parser requires a description of + -- the arguments that will be parsed. This is + -- supplied as a list of `Arg` records, built up + -- using the functions described here. + Arg(..), + Argtype(..), + ArgsComplete(..), + -- ** DataArg and its pseudo-constructors + DataArg, + argDataRequired, argDataOptional, argDataDefaulted, + -- * Argument processing + -- |The argument descriptions are used to parse + -- the command line arguments, and the results + -- of the parse can later be (efficiently) queried + -- to determine program behavior. + + -- ** Getting parse results + -- |The argument parser returns an opaque map + -- from argument index to parsed argument data + -- (plus some convenience information). + ArgRecord, Args(..), + parseArgs, parseArgsIO, + -- ** Using parse results + -- |Query functions permit checking for the existence + -- and values of command-line arguments. + gotArg, ArgType(..), + getArgString, getArgFile, getArgStdio, + getArgInteger, getArgInt, + getArgDouble, getArgFloat, + ArgFileOpener(..), + -- * Misc + ParseArgsException(..), + baseName, parseError, usageError, + System.IO.IOMode(ReadMode, WriteMode, AppendMode)) +where + +import Control.Exception +import Control.Monad +import Control.Monad.ST +import Data.List +import qualified Data.Map as Map +import Data.Maybe +import Data.Typeable +import System.Environment +import System.IO + +-- The main job of this module is to provide parseArgs. +-- See below for its contract. + +-- +-- Provided datatypes. +-- + +-- |The description of an argument, suitable for +-- messages and for parsing. The `argData` field +-- is used both for flags with a data argument, and +-- for positional data arguments. +-- +-- There are two cases: +-- +-- (1) The argument is a flag, in which case at least +-- one of `argAbbr` and `argName` is provided; +-- +-- (2) The argument is positional, in which case neither +-- `argAbbr` nor `argName` are provided, but `argData` is. +-- +-- If none of `argAbbr`, `argName`, or `argData` are +-- provided, this is an error. See also the +-- `argDataRequired`, `argDataOptional`, and +-- `argDataDefaulted` functions below, which are used to +-- generate `argData`. +data Arg a = + Arg { argIndex :: a -- ^Connects the input description + -- to the output argument. + , argAbbr :: Maybe Char -- ^One-character flag name. + , argName :: Maybe String -- ^\"Long name\" of flag. + , argData :: Maybe DataArg -- ^Datum description. + , argDesc :: String -- ^Documentation for the argument. + } + + +-- |The types of an argument carrying data. The constructor +-- argument is used to carry a default value. +-- +-- The constructor argument should really be hidden. +-- Values of this type are normally constructed within +-- the pseudo-constructors pseudo-constructors +-- `argDataRequired`, `argDataOptional`, and +-- `argDataDefaulted`, to which only the constructor +-- function itself is passed. +data Argtype = ArgtypeString (Maybe String) + | ArgtypeInteger (Maybe Integer) + | ArgtypeInt (Maybe Int) + | ArgtypeDouble (Maybe Double) + | ArgtypeFloat (Maybe Float) + + +-- |Information specific to an argument carrying a datum. This +-- is an opaque type, whose instances are constructed using the +-- pseudo-constructors `argDataRequired`, `argDataOptional`, +-- and `argDataDefaulted`. +data DataArg = DataArg { dataArgName :: String -- ^Print name of datum. + , dataArgArgtype :: Argtype -- ^Type of datum. + , dataArgOptional :: Bool -- ^Datum is not required. + } + +-- |Generate the `argData` for the given non-optional argument. +argDataRequired :: String -- ^Datum print name. + -> (Maybe a -> Argtype) -- ^Type constructor for datum. + -> Maybe DataArg -- ^Result is `argData`-ready. +argDataRequired s c = Just (DataArg { dataArgName = s, + dataArgArgtype = c Nothing, + dataArgOptional = False }) + +-- |Generate the `argData` for the given optional argument with no default. +argDataOptional :: String -- ^Datum print name. + -> (Maybe a -> Argtype) -- ^Type constructor for datum. + -> Maybe DataArg -- ^Result is `argData`-ready. +argDataOptional s c = Just (DataArg { dataArgName = s, + dataArgArgtype = c Nothing, + dataArgOptional = True }) + +-- |Generate the `argData` for the given optional argument with the +-- given default. +argDataDefaulted :: String -- ^Datum print name. + -> (Maybe a -> Argtype) -- ^Type constructor for datum. + -> a -- ^Datum default value. + -> Maybe DataArg -- ^Result is `argData`-ready. +argDataDefaulted s c d = Just (DataArg { dataArgName = s, + dataArgArgtype = c (Just d), + dataArgOptional = True }) +-- +-- Returned datatypes. +-- + +-- |The \"kinds of values\" an argument can have. +data Argval = ArgvalFlag -- ^For simple present vs not-present flags. + | ArgvalString String + | ArgvalInteger Integer + | ArgvalInt Int + | ArgvalDouble Double + | ArgvalFloat Float + +-- |The type of the mapping from argument index to value. +newtype ArgRecord a = ArgRecord (Map.Map a Argval) + +-- |The data structure `parseArgs` produces. The key +-- element is the `ArgRecord` `args`. +data Args a = + Args { args :: ArgRecord a -- ^The argument map. + , argsProgName :: String -- ^Basename of 0th argument. + , argsUsage :: String -- ^Full usage string. + , argsRest :: [ String ] -- ^Remaining unprocessed arguments. + } + +-- +-- Exception type. +-- + +-- |This exception is raised with an appropriate error message +-- when argument parsing fails. The first argument is the usage +-- message, the second the actual error message from the parser. +data ParseArgsException = ParseArgsException String String + deriving (Eq, Typeable) + +instance Exception ParseArgsException + +instance Show ParseArgsException where + show (ParseArgsException usage msg) = msg ++ "\n" ++ usage + +-- +-- Implementation. +-- + +-- |True if the described argument is positional. +arg_posn :: (Ord a) => + Arg a -- ^Argument. + -> Bool -- ^True if argument is positional. +arg_posn (Arg { argAbbr = Nothing, + argName = Nothing }) = True +arg_posn _ = False + +-- |True if the described argument is a flag. +arg_flag :: (Ord a) => + Arg a -- ^Argument. + -> Bool -- ^True if argument is a flag. +arg_flag a = not (arg_posn a) + +-- |True if the described argument is optional. +arg_optional :: (Ord a) => + Arg a -- ^Argument. + -> Bool -- ^False if argument is required to be present. +arg_optional (Arg { argData = Just (DataArg { dataArgOptional = b }) }) = b +arg_optional _ = True + +-- |Return the value of a defaulted argument. +arg_default_value :: (Ord a) + => Arg a -- ^Argument. + -> Maybe Argval -- ^Optional default value. +arg_default_value arg@(Arg { argData = Just + (DataArg { dataArgArgtype = da }) }) | + arg_optional arg = + defval da + where + defval (ArgtypeString (Just v)) = Just (ArgvalString v) + defval (ArgtypeInteger (Just v)) = Just (ArgvalInteger v) + defval (ArgtypeInt (Just v)) = Just (ArgvalInt v) + defval (ArgtypeDouble (Just v)) = Just (ArgvalDouble v) + defval (ArgtypeFloat (Just v)) = Just (ArgvalFloat v) + defval _ = Nothing +arg_default_value _ = Nothing + +-- |There's probably a better way to do this. +perhaps b s = if b then s else "" + +-- |Format the described argument as a string. +arg_string :: (Ord a) => + Arg a -- ^Argument to be described. + -> String -- ^String describing argument. +arg_string a@(Arg { argAbbr = abbr, + argName = name, + argData = arg }) = + (optionally "[") ++ + (sometimes flag_abbr abbr) ++ + (perhaps ((isJust abbr) && (isJust name)) ",") ++ + (sometimes flag_name name) ++ + (perhaps ((arg_flag a) && (isJust arg)) " ") ++ + (sometimes data_arg arg) ++ + (optionally "]") + where + sometimes = maybe "" + optionally s = perhaps (arg_optional a) s + flag_name s = "--" ++ s + flag_abbr c = [ '-', c ] + data_arg (DataArg {dataArgName = s}) = "<" ++ s ++ ">" + +-- |Filter out the empty keys for a hash. +filter_keys :: [ (Maybe a, b) ] -- ^List of (optional key, value) pairs. + -> [ (a, b) ] -- ^Pairs with actual keys. +filter_keys l = + foldr check_key [] l + where + check_key (Nothing, _) rest = rest + check_key (Just k, v) rest = (k, v) : rest + +-- |Fail with an error if the argument description is bad +-- for some reason. +argdesc_error :: String -- ^Error message. + -> a -- ^Bogus polymorphic result. +argdesc_error msg = + error ("internal error: argument description: " ++ msg) + +-- |Make a keymap. +keymap_from_list :: (Ord k, Show k) => + [ (k, a) ] -- ^List of key-value pairs. + -- Will be checked for duplicate keys. + -> Map.Map k a -- ^Key-value map. +keymap_from_list l = + foldl add_entry Map.empty l + where + add_entry m (k, a) = + case Map.member k m of + False -> Map.insert k a m + True -> argdesc_error ("duplicate argument description name " ++ + (show k)) + +-- |Make a keymap for looking up a flag argument. +make_keymap :: (Ord a, Ord k, Show k) => + ((Arg a) -> Maybe k) -- ^Mapping from argdesc to flag key. + -> [ Arg a ] -- ^List of argdesc. + -> (Map.Map k (Arg a)) -- ^Map from key to argdesc. +make_keymap f_field args = + (keymap_from_list . + filter_keys . + map (\arg -> (f_field arg, arg))) args + +-- |How \"sloppy\" the parse is. +data ArgsComplete = ArgsComplete -- ^Any extraneous arguments + -- (unparseable from description) + -- will cause the parser to fail. + | ArgsTrailing String -- ^Trailing extraneous arguments are + -- permitted, and will be skipped, + -- saved, and returned. The + -- constructor argument is the + -- name of the args. + | ArgsInterspersed -- ^All extraneous arguments are + -- permitted, and will be skipped, + -- saved, and returned. + +-- |The iteration function is given a state and a list, and +-- expected to produce a new state and list. The function +-- is again invoked with the resulting state and list. +-- When the function returns the empty list, `exhaust` returns +-- the final state produced. +exhaust :: (s -> [e] -> ([e], s)) -- ^Function to iterate. + -> s -- ^Initial state. + -> [e] -- ^Initial list. + -> s -- ^Final state. +exhaust f s [] = s +exhaust f s l = + let (l', s') = f s l + in exhaust f s' l' + +-- |Generate a usage error with the given supplementary message string. +parseError :: String -- ^Usage message. + -> String -- ^Specific error message. + -> a -- ^Bogus polymorphic result. +parseError usage msg = + throw (ParseArgsException usage msg) + +-- |Given a description of the arguments, `parseArgs` produces +-- a map from the arguments to their \"values\" and some other +-- useful byproducts. `parseArgs` requires that the argument +-- descriptions occur in the order 1) flag arguments, 2) required +-- positional arguments, 3) optional positional arguments; otherwise +-- a runtime error will be thrown. +parseArgs :: (Show a, Ord a) => + ArgsComplete -- ^Degree of completeness of parse. + -> [ Arg a ] -- ^Argument descriptions. + -> String -- ^Full program pathname. + -> [ String ] -- ^Incoming program argument list. + -> Args a -- ^Outgoing argument parse results. +parseArgs acomplete argd pathname argv = + runST (do + check_argd + let flag_args = takeWhile arg_flag argd + let posn_args = dropWhile arg_flag argd + let name_hash = make_keymap argName flag_args + let abbr_hash = make_keymap argAbbr flag_args + let prog_name = baseName pathname + let usage = make_usage_string prog_name + let (am, posn, rest) = exhaust (parse usage name_hash abbr_hash) + (Map.empty, posn_args, []) + argv + let required_args = filter (not . arg_optional) argd + unless (and (map (check_present usage am) required_args)) + (error "internal error") + let am' = foldl supply_defaults am argd + return (Args { args = ArgRecord am', + argsProgName = prog_name, + argsUsage = usage, + argsRest = rest })) + where + supply_defaults am ad@(Arg { argIndex = k }) = + case Map.lookup k am of + Just _ -> am + Nothing -> case arg_default_value ad of + Just v -> Map.insert k v am + Nothing -> am + check_present usage am ad@(Arg { argIndex = k }) = + case Map.lookup k am of + Just _ -> True + Nothing -> parseError usage ("missing required argument " ++ + (arg_string ad)) + --- Check for various possible misuses. + check_argd :: ST s () + check_argd = do + --- Order must be flags, posn args, optional posn args + let residue = dropWhile arg_flag argd + let residue' = dropWhile arg_fixed_posn residue + let residue'' = dropWhile arg_opt_posn residue' + unless (null residue'') + (argdesc_error "argument description in wrong order") + --- No argument may be "nullary". + when (or (map arg_nullary argd)) + (argdesc_error "bogus 'nothing' argument") + return () + where + arg_fixed_posn a = (arg_posn a) && (not (arg_optional a)) + arg_opt_posn a = (arg_posn a) && (arg_optional a) + arg_nullary (Arg { argName = Nothing, + argAbbr = Nothing, + argData = Nothing }) = True + arg_nullary _ = False + --- Generate a usage message string + make_usage_string prog_name = + summary_line ++ arg_lines + where + flag_args = filter arg_flag argd + posn_args = filter arg_posn argd + n = maximum (map (length . arg_string) argd) + --- top (summary) line + summary_line = + "usage: " ++ prog_name ++ + perhaps + (not (null flag_args)) + " [options]" ++ + perhaps + (not (null posn_args)) + (" " ++ unwords (map arg_string posn_args)) ++ + (case acomplete of + ArgsComplete -> "" + ArgsTrailing s -> " [--] [" ++ s ++ " ...]" + ArgsInterspersed -> " ... [--] ...") ++ "\n" + --- argument lines + arg_lines = concatMap (arg_line n) argd where + arg_line n a = + let s = arg_string a in + " " ++ s ++ + replicate (n - (length s)) ' ' ++ + " " ++ argDesc a ++ "\n" + --- simple recursive-descent parser + parse _ _ _ av@(_, _, []) [] = ([], av) + parse usage _ _ av [] = + case acomplete of + ArgsComplete -> parseError usage "unexpected extra arguments" + _ -> ([], av) + parse usage name_hash abbr_hash (am, posn, rest) av@(aa : aas) = + case aa of + "--" -> case acomplete of + ArgsComplete -> parseError usage + ("unexpected -- " ++ + "(extra arguments not allowed)") + _ -> ([], (am, posn, (rest ++ aas))) + s@('-' : '-' : name) -> + case Map.lookup name name_hash of + Just ad -> peel s ad aas + Nothing -> + case acomplete of + ArgsInterspersed -> + (aas, (am, posn, rest ++ ["--" ++ name])) + _ -> parseError usage + ("unknown argument --" ++ name) + ('-' : abbr : abbrs) -> + case Map.lookup abbr abbr_hash of + Just ad -> + let p@(args', state') = peel ['-', abbr] ad aas + in case abbrs of + [] -> p + ('-' : _) -> parseError usage + ("bad internal '-' in argument " ++ aa) + _ -> (['-' : abbrs] ++ args', state') + Nothing -> + case acomplete of + ArgsInterspersed -> + (['-' : abbrs] ++ aas, + (am, posn, rest ++ [['-', abbr]])) + _ -> parseError usage + ("unknown argument -" ++ [abbr]) + aa -> case posn of + (ad@(Arg { argData = Just adata }) : ps) -> + let (argl', (am', _, rest')) = + peel_process (dataArgName adata) ad av + in (argl', (am', ps, rest')) + [] -> case acomplete of + ArgsComplete -> parseError usage + ("unexpected argument " ++ aa) + _ -> (aas, (am, [], rest ++ [aa])) + where + add_entry s m (k, a) = + case Map.member k m of + False -> Map.insert k a m + True -> parseError usage ("duplicate argument " ++ s) + peel name ad@(Arg { argData = Nothing, argIndex = index }) argl = + let am' = add_entry name am (index, ArgvalFlag) + in (argl, (am', posn, rest)) + peel name (Arg { argData = Just (DataArg {}) }) [] = + parseError usage (name ++ " is missing its argument") + peel name ad argl = peel_process name ad argl + peel_process name + ad@(Arg { argData = Just (DataArg { + dataArgArgtype = atype }), + argIndex = index }) + (a : argl) = + let read_arg constructor kind = + case reads a of + [(v, "")] -> constructor v + _ -> parseError usage ("argument " ++ + a ++ " to " ++ name ++ + " is not " ++ kind) + v = case atype of + ArgtypeString _ -> ArgvalString a + ArgtypeInteger _ -> read_arg ArgvalInteger + "an integer" + ArgtypeInt _ -> read_arg ArgvalInt "an int" + ArgtypeDouble _ -> read_arg ArgvalDouble "a double" + ArgtypeFloat _ -> read_arg ArgvalFloat "a float" + am' = add_entry name am (index, v) + in (argl, (am', posn, rest)) + + +-- |Most of the time, you just want the environment's +-- arguments and are willing to live in the IO monad. +-- This version of `parseArgs` digs the pathname and arguments +-- out of the system directly. +parseArgsIO :: (Show a, Ord a) => + ArgsComplete -- ^Degree of completeness of parse. + -> [ Arg a ] -- ^Argument descriptions. + -> IO (Args a) -- ^Argument parse results. +parseArgsIO acomplete argd = do + argv <- getArgs + pathname <- getProgName + return (parseArgs acomplete argd pathname argv) + + +-- |Check whether a given optional argument was supplied. Works on all types. +gotArg :: (Ord a) => + Args a -- ^Parsed arguments. + -> a -- ^Index of argument to be checked for. + -> Bool -- ^True if the arg was present. +gotArg (Args { args = ArgRecord am }) k = + case Map.lookup k am of + Just _ -> True + Nothing -> False + +-- |Type of values that can be parsed by the argument parser. +class ArgType b where + -- |Fetch an argument's value if it is present. + getArg :: (Show a, Ord a) + => Args a -- ^Parsed arguments. + -> a -- ^Index of argument to be retrieved. + -> Maybe b -- ^Argument value if present. + -- |Fetch the value of a required argument. + getRequiredArg :: (Show a, Ord a) + => Args a -- ^Parsed arguments. + -> a -- ^Index of argument to be retrieved. + -> b -- ^Argument value. + getRequiredArg args index = + case getArg args index of + Just v -> v + Nothing -> error ("internal error: required argument " + ++ show index ++ "not supplied") + +getArgPrimitive decons (Args { args = ArgRecord am }) k = + case Map.lookup k am of + Just v -> Just (decons v) + Nothing -> Nothing + +instance ArgType ([] Char) where + getArg = getArgPrimitive (\(ArgvalString s) -> s) + +-- |[Deprecated] Return the `String` value, if any, of the given argument. +getArgString :: (Show a, Ord a) => + Args a -- ^Parsed arguments. + -> a -- ^Index of argument to be retrieved. + -> Maybe String -- ^Argument value if present. +getArgString = getArg + +instance ArgType Integer where + getArg = getArgPrimitive (\(ArgvalInteger i) -> i) + +-- |[Deprecated] Return the `Integer` value, if any, of the given argument. +getArgInteger :: (Show a, Ord a) => + Args a -- ^Parsed arguments. + -> a -- ^Index of argument to be retrieved. + -> Maybe Integer -- ^Argument value if present. +getArgInteger = getArg + +instance ArgType Int where + getArg = getArgPrimitive (\(ArgvalInt i) -> i) + +-- |[Deprecated] Return the `Int` value, if any, of the given argument. +getArgInt :: (Show a, Ord a) => + Args a -- ^Parsed arguments. + -> a -- ^Index of argument to be retrieved. + -> Maybe Int -- ^Argument value if present. +getArgInt = getArg + +instance ArgType Double where + getArg = getArgPrimitive (\(ArgvalDouble i) -> i) + +-- |[Deprecated] Return the `Double` value, if any, of the given argument. +getArgDouble :: (Show a, Ord a) => + Args a -- ^Parsed arguments. + -> a -- ^Index of argument to be retrieved. + -> Maybe Double -- ^Argument value if present. +getArgDouble = getArg + +instance ArgType Float where + getArg = getArgPrimitive (\(ArgvalFloat i) -> i) + +-- |[Deprecated] Return the `Float` value, if any, of the given argument. +getArgFloat :: (Show a, Ord a) => + Args a -- ^Parsed arguments. + -> a -- ^Index of argument to be retrieved. + -> Maybe Float -- ^Argument value if present. +getArgFloat = getArg + +-- |`ArgType` instance for opening a file from its string name. +newtype ArgFileOpener = ArgFileOpener { + argFileOpener :: IOMode -> IO Handle -- ^Function to open the file + } + +instance ArgType ArgFileOpener where + getArg args index = + case getArg args index of + Nothing -> Nothing + Just s -> Just (ArgFileOpener { argFileOpener = openFile s }) + +-- |[Deprecated] Treat the `String` value, if any, of the given argument as +-- a file handle and try to open it as requested. +getArgFile :: (Show a, Ord a) => + Args a -- ^Parsed arguments. + -> a -- ^Index of argument to be retrieved. + -> IOMode -- ^IO mode the file should be opened in. + -> IO (Maybe Handle) -- ^Handle of opened file, if the argument + -- was present. +getArgFile args k m = + case getArg args k of + Just fo -> (do h <- argFileOpener fo m; return (Just h)) + Nothing -> return Nothing + + +-- |Treat the `String` value, if any, of the given argument as a +-- file handle and try to open it as requested. If not +-- present, substitute the appropriate one of stdin or +-- stdout as indicated by `IOMode`. +getArgStdio :: (Show a, Ord a) => + Args a -- ^Parsed arguments. + -> a -- ^Index of argument to be retrieved. + -> IOMode -- ^IO mode the file should be opened in. + -- Must not be `ReadWriteMode`. + -> IO Handle -- ^Appropriate file handle. +getArgStdio args k m = + case getArg args k of + Just s -> openFile s m + Nothing -> + case m of + ReadMode -> return stdin + WriteMode -> return stdout + AppendMode -> return stdout + ReadWriteMode -> + error ("internal error: tried to open stdio " + ++ "in ReadWriteMode") + +--- +--- Misc +--- + +-- |Return the filename part of a pathname. +-- Unnecessarily efficient implementation does a single +-- tail-call traversal with no construction. +baseName :: String -- ^Pathname. + -> String -- ^Rightmost component of pathname. +baseName s = + let s' = dropWhile (/= '/') s in + if null s' then s else baseName (tail s') + + +-- |Generate a usage error with the given supplementary message string. +usageError :: (Ord a) => Args a -> String -> b +usageError args msg = error (argsUsage args ++ "\n" ++ msg) diff --git a/lib/Timing.hs b/lib/Timing.hs --- a/lib/Timing.hs +++ b/lib/Timing.hs @@ -1,10 +1,9 @@ module Timing (time, showTime, prettyTime) where -import GHC.Exts (traceEvent) import System.CPUTime import System.Time - +import Debug.Trace -- Time ----------------------------------------------------------------------- data Time @@ -63,8 +62,8 @@ {-# NOINLINE time #-} time p = do start <- getTime - traceEvent "dph-examples: start timing" + traceEventIO "dph-examples: start timing" x <- p - traceEvent "dph-examples: finished timing" + traceEventIO "dph-examples: finished timing" end <- getTime return (x, end `minus` start) diff --git a/real/NBody/Batch/Config.hs b/real/NBody/Batch/Config.hs deleted file mode 100644 --- a/real/NBody/Batch/Config.hs +++ /dev/null @@ -1,57 +0,0 @@ - -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 deleted file mode 100644 --- a/real/NBody/Batch/MainArgs.hs +++ /dev/null @@ -1,87 +0,0 @@ - -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 deleted file mode 100644 --- a/real/NBody/Common/Body.hs +++ /dev/null @@ -1,103 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/real/NBody/Common/Dump.hs +++ /dev/null @@ -1,32 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/real/NBody/Common/Util.hs +++ /dev/null @@ -1,18 +0,0 @@ - -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 deleted file mode 100644 --- a/real/NBody/Common/World.hs +++ /dev/null @@ -1,48 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/real/NBody/Gloss/Config.hs +++ /dev/null @@ -1,74 +0,0 @@ - -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 deleted file mode 100644 --- a/real/NBody/Gloss/Draw.hs +++ /dev/null @@ -1,104 +0,0 @@ -{-# 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/MainArgs.hs b/real/NBody/Gloss/MainArgs.hs deleted file mode 100644 --- a/real/NBody/Gloss/MainArgs.hs +++ /dev/null @@ -1,113 +0,0 @@ - -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/MainBatch.hs b/real/NBody/MainBatch.hs deleted file mode 100644 --- a/real/NBody/MainBatch.hs +++ /dev/null @@ -1,95 +0,0 @@ -{-# LANGUAGE ParallelListComp, BangPatterns #-} - -import Batch.MainArgs -import Batch.Config -import Common.Dump -import Common.World -import Common.Body -import Common.Util -import Solver -import Timing -import Points2D.Generate -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 mainBatch config calcAccels worldStart - - --- | Run the simulation in batch mode. -mainBatch :: Config -> Solver -> 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 if worldSteps world' < configMaxSteps config - then go world' - else 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/MainGloss.hs b/real/NBody/MainGloss.hs deleted file mode 100644 --- a/real/NBody/MainGloss.hs +++ /dev/null @@ -1,151 +0,0 @@ -{-# 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/Solver.hs b/real/NBody/Solver.hs deleted file mode 100644 --- a/real/NBody/Solver.hs +++ /dev/null @@ -1,59 +0,0 @@ - --- | 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 deleted file mode 100644 --- a/real/NBody/Solver/ListBH/Solver.hs +++ /dev/null @@ -1,161 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/real/NBody/Solver/NestedBH/Solver.hs +++ /dev/null @@ -1,154 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/real/NBody/Solver/VectorBH/Solver.hs +++ /dev/null @@ -1,163 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/real/NBody/Solver/VectorNaive/Solver.hs +++ /dev/null @@ -1,18 +0,0 @@ - --- | 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 deleted file mode 100644 --- a/spectral/QuickHull/dph/Main.hs +++ /dev/null @@ -1,59 +0,0 @@ - -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 deleted file mode 100644 --- a/spectral/QuickHull/dph/QuickHullVect.hs +++ /dev/null @@ -1,42 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/spectral/QuickHull/lib/SVG.hs +++ /dev/null @@ -1,34 +0,0 @@ - -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 deleted file mode 100644 --- a/spectral/QuickHull/vector/Main.hs +++ /dev/null @@ -1,71 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/spectral/QuickHull/vector/QuickHullIO.hs +++ /dev/null @@ -1,269 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/spectral/QuickHull/vector/QuickHullSplit.hs +++ /dev/null @@ -1,38 +0,0 @@ - -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 deleted file mode 100644 --- a/spectral/QuickHull/vector/QuickHullVector.hs +++ /dev/null @@ -1,80 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/spectral/QuickSort/dph/Main.hs +++ /dev/null @@ -1,50 +0,0 @@ - -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 deleted file mode 100644 --- a/spectral/QuickSort/dph/QuickSortVect.hs +++ /dev/null @@ -1,25 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/spectral/SMVM/dph/Main.hs +++ /dev/null @@ -1,102 +0,0 @@ -{-# 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 deleted file mode 100644 --- a/spectral/SMVM/dph/SMVMVectorised.hs +++ /dev/null @@ -1,18 +0,0 @@ - -{-# 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 :] -