dph-examples 0.5.1.2 → 0.6.0.1
raw patch · 93 files changed
+3738/−2734 lines, 93 filesdep +HUnitdep +containersdep +dph-lifted-vsegdep −dph-pardep −dph-prim-seqdep −dph-seqdep ~basedep ~dph-basedep ~dph-prim-parnew-component:exe:dph-imaginary-wordsnew-component:exe:dph-real-nbodynew-component:exe:dph-real-nbody-glossnew-component:exe:dph-smoke-boolnew-component:exe:dph-smoke-concatnew-component:exe:dph-smoke-evensnew-component:exe:dph-smoke-indicesnew-component:exe:dph-smoke-ranknew-component:exe:dph-smoke-reversenew-component:exe:dph-smoke-sumsqnew-component:exe:dph-spectral-dotpnew-component:exe:dph-spectral-quickhullnew-component:exe:dph-spectral-quickhull-vectornew-component:exe:dph-spectral-quicksortnew-component:exe:dph-spectral-smvm
Dependencies added: HUnit, containers, dph-lifted-vseg
Dependencies removed: dph-par, dph-prim-seq, dph-seq, parseargs
Dependency ranges changed: base, dph-base, dph-prim-par, gloss, old-time, vector
Files
- LICENSE +25/−20
- dph-examples.cabal +106/−179
- examples/imaginary/Words/Main.hs +37/−0
- examples/imaginary/Words/Vectorised.hs +113/−0
- examples/real/NBody/Batch/Config.hs +65/−0
- examples/real/NBody/Batch/MainArgs.hs +99/−0
- examples/real/NBody/Common/Body.hs +103/−0
- examples/real/NBody/Common/Dump.hs +46/−0
- examples/real/NBody/Common/Util.hs +18/−0
- examples/real/NBody/Common/World.hs +48/−0
- examples/real/NBody/Gloss/Config.hs +74/−0
- examples/real/NBody/Gloss/Draw.hs +104/−0
- examples/real/NBody/Gloss/MainArgs.hs +113/−0
- examples/real/NBody/MainBatch.hs +102/−0
- examples/real/NBody/MainGloss.hs +151/−0
- examples/real/NBody/Solver.hs +59/−0
- examples/real/NBody/Solver/ListBH/Solver.hs +161/−0
- examples/real/NBody/Solver/NestedBH/Solver.hs +156/−0
- examples/real/NBody/Solver/VectorBH/Solver.hs +163/−0
- examples/real/NBody/Solver/VectorNaive/Solver.hs +18/−0
- examples/smoke/data/Bool/Main.hs +7/−0
- examples/smoke/data/Bool/Vectorised.hs +28/−0
- examples/smoke/prims/Concat/Main.hs +23/−0
- examples/smoke/prims/Concat/Vectorised.hs +38/−0
- examples/smoke/prims/Evens/dph/Main.hs +48/−0
- examples/smoke/prims/Evens/dph/Vector.hs +8/−0
- examples/smoke/prims/Evens/dph/Vectorised.hs +16/−0
- examples/smoke/prims/SumSquares/dph/Main.hs +35/−0
- examples/smoke/prims/SumSquares/dph/Vector.hs +11/−0
- examples/smoke/prims/SumSquares/dph/Vectorised.hs +16/−0
- examples/smoke/sharing/Indices/Main.hs +42/−0
- examples/smoke/sharing/Indices/Vectorised.hs +31/−0
- examples/smoke/sharing/Rank/Main.hs +40/−0
- examples/smoke/sharing/Rank/Vectorised.hs +18/−0
- examples/smoke/sharing/Reverse/Main.hs +58/−0
- examples/smoke/sharing/Reverse/Vectorised.hs +28/−0
- examples/spectral/DotProduct/dph/Main.hs +45/−0
- examples/spectral/DotProduct/dph/Vector.hs +7/−0
- examples/spectral/DotProduct/dph/Vectorised.hs +15/−0
- examples/spectral/QuickHull/dph/Main.hs +59/−0
- examples/spectral/QuickHull/dph/Vectorised.hs +43/−0
- examples/spectral/QuickHull/lib/SVG.hs +34/−0
- examples/spectral/QuickHull/vector/Main.hs +71/−0
- examples/spectral/QuickHull/vector/QuickHullIO.hs +269/−0
- examples/spectral/QuickHull/vector/QuickHullSplit.hs +38/−0
- examples/spectral/QuickHull/vector/QuickHullVector.hs +80/−0
- examples/spectral/QuickSort/dph/Main.hs +50/−0
- examples/spectral/QuickSort/dph/Vectorised.hs +25/−0
- examples/spectral/SMVM/dph/Main.hs +102/−0
- examples/spectral/SMVM/dph/Vectorised.hs +17/−0
- imaginary/DotProduct/dph/DotProductVector.hs +0/−7
- imaginary/DotProduct/dph/DotProductVectorised.hs +0/−15
- imaginary/DotProduct/dph/Main.hs +0/−45
- imaginary/Evens/dph/EvensVector.hs +0/−8
- imaginary/Evens/dph/EvensVectorised.hs +0/−16
- imaginary/Evens/dph/Main.hs +0/−45
- imaginary/Reverse/dph/Main.hs +0/−58
- imaginary/Reverse/dph/ReverseVector.hs +0/−21
- imaginary/Reverse/dph/ReverseVectorised.hs +0/−29
- imaginary/SumSquares/dph/Main.hs +0/−34
- imaginary/SumSquares/dph/SumSquaresVector.hs +0/−11
- imaginary/SumSquares/dph/SumSquaresVectorised.hs +0/−15
- lib/Points2D/Generate.hs +2/−2
- lib/Points2D/Types.hs +3/−0
- lib/System/Console/ParseArgs.hs +667/−0
- lib/Timing.hs +3/−4
- real/NBody/Batch/Config.hs +0/−57
- real/NBody/Batch/MainArgs.hs +0/−87
- real/NBody/Common/Body.hs +0/−103
- real/NBody/Common/Dump.hs +0/−32
- real/NBody/Common/Util.hs +0/−18
- real/NBody/Common/World.hs +0/−48
- real/NBody/Gloss/Config.hs +0/−74
- real/NBody/Gloss/Draw.hs +0/−104
- real/NBody/Gloss/MainArgs.hs +0/−113
- real/NBody/MainBatch.hs +0/−95
- real/NBody/MainGloss.hs +0/−151
- real/NBody/Solver.hs +0/−59
- real/NBody/Solver/ListBH/Solver.hs +0/−161
- real/NBody/Solver/NestedBH/Solver.hs +0/−154
- real/NBody/Solver/VectorBH/Solver.hs +0/−163
- real/NBody/Solver/VectorNaive/Solver.hs +0/−18
- spectral/QuickHull/dph/Main.hs +0/−59
- spectral/QuickHull/dph/QuickHullVect.hs +0/−42
- spectral/QuickHull/lib/SVG.hs +0/−34
- spectral/QuickHull/vector/Main.hs +0/−71
- spectral/QuickHull/vector/QuickHullIO.hs +0/−269
- spectral/QuickHull/vector/QuickHullSplit.hs +0/−38
- spectral/QuickHull/vector/QuickHullVector.hs +0/−80
- spectral/QuickSort/dph/Main.hs +0/−50
- spectral/QuickSort/dph/QuickSortVect.hs +0/−25
- spectral/SMVM/dph/Main.hs +0/−102
- spectral/SMVM/dph/SMVMVectorised.hs +0/−18
LICENSE view
@@ -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.+
dph-examples.cabal view
@@ -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
+ examples/imaginary/Words/Main.hs view
@@ -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"+
+ examples/imaginary/Words/Vectorised.hs view
@@ -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+
+ examples/real/NBody/Batch/Config.hs view
@@ -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 }
+ examples/real/NBody/Batch/MainArgs.hs view
@@ -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" }++ ]
+ examples/real/NBody/Common/Body.hs view
@@ -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)
+ examples/real/NBody/Common/Dump.hs view
@@ -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
+ examples/real/NBody/Common/Util.hs view
@@ -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
+ examples/real/NBody/Common/World.hs view
@@ -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' }++
+ examples/real/NBody/Gloss/Config.hs view
@@ -0,0 +1,74 @@++module Gloss.Config+ ( Config (..)+ , loadConfig)+where+import Gloss.MainArgs+import System.Console.ParseArgs+import Data.Maybe++-- | Program config.+data Config+ = Config {+ + -- How to present the output.+ configWindowSize :: Maybe Int+ , configShouldDrawTree :: Bool+ , configRate :: Int++ -- What solver to use+ , configSolverName :: String++ -- System setup+ , configBodyCount :: Int+ , configBodyMass :: Double+ , configTimeStep :: Double+ , configEpsilon :: Double++ -- Initial conditions.+ , configStartDiscSize :: Double+ , configStartSpeed :: Double++ -- Terminating conditions.+ , configMaxSteps :: Maybe Int+ + -- dump points to file+ , configDumpFinal :: Maybe FilePath }+ ++-- | Load program config from its command line arguments. +loadConfig :: Args MainArg -> Config+loadConfig args+ = let mWindowSize = getArgInt args ArgGloss+ Just solverName = getArgString args ArgSolver+ shouldDrawTree = gotArg args ArgDrawTree+ Just timeStep = getArgDouble args ArgTimeStep+ Just rate = getArgInt args ArgRate+ Just bodyCount = getArgInt args ArgBodyCount+ Just bodyMass = getArgDouble args ArgBodyMass+ Just epsilon = getArgDouble args ArgEpsilon+ Just discSize = getArgDouble args ArgDiscSize+ Just startSpeed = getArgDouble args ArgStartSpeed++ mMaxSteps = getArgInt args ArgMaxSteps+ mFilePath = getArgString args ArgDumpFinal++ checkMode x+ = if not (isJust mWindowSize || isJust mMaxSteps)+ then error "you must specify either --max-steps or --gloss <window size :: Int>"+ else x ++ in checkMode $+ Config+ { configWindowSize = mWindowSize+ , configShouldDrawTree = shouldDrawTree+ , configRate = rate+ , configSolverName = solverName+ , configBodyCount = bodyCount+ , configBodyMass = bodyMass+ , configTimeStep = timeStep+ , configEpsilon = epsilon+ , configStartDiscSize = discSize+ , configStartSpeed = startSpeed+ , configMaxSteps = mMaxSteps + , configDumpFinal = mFilePath }
+ examples/real/NBody/Gloss/Draw.hs view
@@ -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
+ examples/real/NBody/Gloss/MainArgs.hs view
@@ -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" }++ ]
+ examples/real/NBody/MainBatch.hs view
@@ -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)++
+ examples/real/NBody/MainGloss.hs view
@@ -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++
+ examples/real/NBody/Solver.hs view
@@ -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'
+ examples/real/NBody/Solver/ListBH/Solver.hs view
@@ -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+
+ examples/real/NBody/Solver/NestedBH/Solver.hs view
@@ -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+
+ examples/real/NBody/Solver/VectorBH/Solver.hs view
@@ -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++
+ examples/real/NBody/Solver/VectorNaive/Solver.hs view
@@ -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)
+ examples/smoke/data/Bool/Main.hs view
@@ -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]
+ examples/smoke/data/Bool/Vectorised.hs view
@@ -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))
+ examples/smoke/prims/Concat/Main.hs view
@@ -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)+ ]+
+ examples/smoke/prims/Concat/Vectorised.hs view
@@ -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'
+ examples/smoke/prims/Evens/dph/Main.hs view
@@ -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> <length>"+ , " alg one of " ++ show ["vectorised", "vector"] ]++run alg len+ = do let vec = V.fromList [0..len - 1]+ runAlg alg vec+ +runAlg "vectorised" vec+ = do let arr = P.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'+
+ examples/smoke/prims/Evens/dph/Vector.hs view
@@ -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+
+ examples/smoke/prims/Evens/dph/Vectorised.hs view
@@ -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))+
+ examples/smoke/prims/SumSquares/dph/Main.hs view
@@ -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> <length>"+ , " alg one of " ++ show ["vectorised", "vector"] ]+ +run alg num+ = do (result, tElapsed) <- runAlg alg num++ putStr $ prettyTime tElapsed+ print result++runAlg "vectorised" num+ = time $ let result = Z.sumSq num+ in result `seq` return result++runAlg "vector" num+ = time $ let result = V.sumSq num+ in result `seq` return result+
+ examples/smoke/prims/SumSquares/dph/Vector.hs view
@@ -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+
+ examples/smoke/prims/SumSquares/dph/Vectorised.hs view
@@ -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)))
+ examples/smoke/sharing/Indices/Main.hs view
@@ -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 <algorithm> <count>\n"+ , " algorithm one of " ++ show ["vectorised"]+ , ""]
+ examples/smoke/sharing/Indices/Vectorised.hs view
@@ -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 :])
+ examples/smoke/sharing/Rank/Main.hs view
@@ -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 <algorithm> <count>\n"+ , " algorithm one of " ++ show ["vectorised"]+ , ""]
+ examples/smoke/sharing/Rank/Vectorised.hs view
@@ -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 #-}+
+ examples/smoke/sharing/Reverse/Main.hs view
@@ -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 <algorithm> <count>\n"+ , " algorithm one of " ++ show ["vectorised"]+ , ""]
+ examples/smoke/sharing/Reverse/Vectorised.hs view
@@ -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 :])
+ examples/spectral/DotProduct/dph/Main.hs view
@@ -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> <length>"+ , " alg one of " ++ show ["vectorised", "vector"] ]+ +run alg len+ = do let vec1 = randomishDoubles len 0 1 1234+ let vec2 = randomishDoubles len 0 1 12345++ (result, tElapsed) <- runAlg alg vec1 vec2++ putStr $ prettyTime tElapsed+ putStr $ (take 12 $ show result) ++ "\n"++runAlg "vectorised" vec1 vec2 + = do let arr1 = 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+
+ examples/spectral/DotProduct/dph/Vector.hs view
@@ -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
+ examples/spectral/DotProduct/dph/Vectorised.hs view
@@ -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)
+ examples/spectral/QuickHull/dph/Main.hs view
@@ -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 <points> [out.svg]" ]+++-- | Run the benchmark.+run :: Int -- ^ How many points to use.+ -> Maybe String -- ^ File name to dump an SVG of the output to.+ -> IO ()+ +run pointCount mFileSVG+ = do+ vPoints <- pointsPArrayOfUArray+ $ genPointsDisc pointCount (400, 400) 350 ++ -- Force points to create the input vector.+ vPoints `seq` return ()++ -- Compute the convex hull.+ (vHull, tElapsed)+ <- time + $ let vHull = quickhullPA vPoints+ in vHull `seq` return vHull+ + -- Print how long it took.+ putStr $ prettyTime tElapsed+ + -- If we were asked for an SVG then write it out to file.+ maybe (return ())+ (\fileSVG -> + writeFile fileSVG+ $ makeSVG (roundPoints $ P.toList vPoints) + (roundPoints $ P.toList vHull))+ mFileSVG
+ examples/spectral/QuickHull/dph/Vectorised.hs view
@@ -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 #-}+
+ examples/spectral/QuickHull/lib/SVG.hs view
@@ -0,0 +1,34 @@++module SVG where++-- Making a SVG diagram of the points and hull+makeSVG :: [(Int, Int)] -> [(Int, Int)] -> String+makeSVG points hull+ = unlines+ $ [ "<svg width=\"100%\" height=\"100%\" version=\"1.1\" xmlns=\"http://www.w3.org/2000/svg\">" ]+ ++ [svgPolygon hull]+ ++ map svgPoint points+ ++ map svgPointHull hull+ ++ ["</svg>"]++svgPolygon :: [(Int, Int)] -> String+svgPolygon points+ = "<polygon"+ ++ " points=\"" ++ (concat [show x ++ "," ++ show y ++ " " | (x, y) <- points]) ++ "\"" + ++ " style=\"fill:#d0d0ff;stroke:#000000;stroke-width:1\""+ ++ "/>"++svgPoint :: (Int, Int) -> String+svgPoint (x, y)+ = "<circle cx=\"" ++ show x ++ "\" cy=\"" ++ show y ++ "\" r=\"0.5\""+ ++ " style=\"stroke:#000000\""+ ++ "/>"++svgPointHull :: (Int, Int) -> String+svgPointHull (x, y)+ = "<circle cx=\"" ++ show x ++ "\" cy=\"" ++ show y ++ "\" r=\"1\""+ ++ " style=\"stroke:#ff0000\""+ ++ "/>"+ +roundPoints :: [(Double, Double)] -> [(Int, Int)]+roundPoints ps = [(round x, round y) | (x, y) <- ps]
+ examples/spectral/QuickHull/vector/Main.hs view
@@ -0,0 +1,71 @@+{-# LANGUAGE PatternGuards #-}++import qualified QuickHullVector+import qualified QuickHullIO+import qualified QuickHullSplit+import Points2D.Types+import Points2D.Generate+import Timing+import SVG++import System.Environment+import Data.Function+import qualified Data.Vector.Unboxed as V+import Data.Vector.Unboxed (Vector)+import Control.Monad++algs = [ ("vector", (\v -> return $ QuickHullVector.quickHull v))+ , ("io", QuickHullIO.quickHull)+ , ("split", (\v -> return $ QuickHullSplit.quickHull v)) ]++parseArgs args+ | [alg, strCount] <- args+ , Just fun <- lookup alg algs+ = Just (fun, read strCount, Nothing)++ | [alg, strCount, file] <- args+ , Just fun <- lookup alg algs+ = Just (fun, read strCount, Just file)++ | otherwise+ = Nothing+++main :: IO ()+main+ = do argStrs <- getArgs+ case parseArgs argStrs of+ Just args -> run args+ _ -> putStr $ unlines+ [ "usage: quickhull <alg> <points> [out.svg]"+ , " algs: " ++ (show $ map fst algs) ++ "\n" ]++run :: ( Vector Point -> IO (Vector Point)+ , Int+ , Maybe FilePath)+ -> IO ()+ +run (fun, pointCount, mFileSVG) + = do+ let vPoints = genPointsDisc pointCount (400, 400) 350 ++ -- Force points to create the input vector.+ V.force vPoints `seq` return ()++ -- Compute the convex hull.+ (vHull, tElapsed)+ <- time+ $ do vHull <- fun vPoints+ vHull `seq` return vHull+ + -- Print how long it took + putStr $ prettyTime tElapsed++ -- If we were asked for an SVG then write it out to file.+ maybe (return ())+ (\fileSVG -> + writeFile fileSVG+ $ makeSVG (roundPoints $ V.toList vPoints) + (roundPoints $ V.toList vHull))+ mFileSVG+
+ examples/spectral/QuickHull/vector/QuickHullIO.hs view
@@ -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) +
+ examples/spectral/QuickHull/vector/QuickHullSplit.hs view
@@ -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)
+ examples/spectral/QuickHull/vector/QuickHullVector.hs view
@@ -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)+
+ examples/spectral/QuickSort/dph/Main.hs view
@@ -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 <length>" ]+++-- | 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)
+ examples/spectral/QuickSort/dph/Vectorised.hs view
@@ -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)
+ examples/spectral/SMVM/dph/Main.hs view
@@ -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 <file>" ]++run fileName+ = do (matrix, vector) <- loadPA fileName++ matrix `seq` return ()+ vector `seq` return ()++ -- Multiply sparse matrix by the dense vector.+ (vResult, tElapsed)+ <- time $ let result = smvmPA matrix vector+ in P.nf result `seq` return result+ + -- Print how long it took.+ putStr $ prettyTime tElapsed++ -- Print some info about the test setup.+ putStrLn $ "vector length = " ++ show (U.length (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)
+ examples/spectral/SMVM/dph/Vectorised.hs view
@@ -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 :]+
− imaginary/DotProduct/dph/DotProductVector.hs
@@ -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
− imaginary/DotProduct/dph/DotProductVectorised.hs
@@ -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)
− imaginary/DotProduct/dph/Main.hs
@@ -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> <length>"- , " alg one of " ++ show ["vectorised", "vector"] ]- -run alg len- = do let vec1 = randomishDoubles len 0 1 1234- let vec2 = randomishDoubles len 0 1 12345-- (result, tElapsed) <- runAlg alg vec1 vec2-- putStr $ prettyTime tElapsed- putStr $ (take 12 $ show result) ++ "\n"--runAlg "vectorised" vec1 vec2 - = do let arr1 = fromUArrPA' vec1- let arr2 = fromUArrPA' vec2- arr1 `seq` arr2 `seq` return ()-- time $ let result = Z.dotPA arr1 arr2- in result `seq` return result--runAlg "vector" vec1 vec2- = do vec1 `seq` vec2 `seq` return ()-- time $ let result = V.dotV vec1 vec2 - in result `seq` return result-
− imaginary/Evens/dph/EvensVector.hs
@@ -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-
− imaginary/Evens/dph/EvensVectorised.hs
@@ -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))-
− imaginary/Evens/dph/Main.hs
@@ -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> <length>"- , " alg one of " ++ show ["vectorised", "vector"] ]--run alg len- = do let vec = V.fromList [0..len - 1]- runAlg alg vec- -runAlg "vectorised" vec- = do let arr = P.fromUArrPA' vec- arr `seq` return ()-- (arr', tElapsed)- <- time $ let arr' = evensPA arr- in arr' `seq` return arr'- - putStr $ prettyTime tElapsed- putStrLn $ show $ P.length arr'--runAlg "vector" vec- = do vec `seq` return ()- - (vec', tElapsed)- <- time $ let vec' = evensV vec- in vec' `seq` return vec'- - putStr $ prettyTime tElapsed- putStrLn $ show $ V.length vec'-
− imaginary/Reverse/dph/Main.hs
@@ -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 <algorithm> <count>\n"- , " algorithm one of " ++ show ["vectorised"]- , ""]
− imaginary/Reverse/dph/ReverseVector.hs
@@ -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-
− imaginary/Reverse/dph/ReverseVectorised.hs
@@ -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 :])
− imaginary/SumSquares/dph/Main.hs
@@ -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> <length>"- , " alg one of " ++ show ["vectorised", "vector"] ]- -run alg num- = do (result, tElapsed) <- runAlg alg num-- putStr $ prettyTime tElapsed- print result--runAlg "vectorised" num- = time $ let result = Z.sumSq num- in result `seq` return result--runAlg "vector" num- = time $ let result = V.sumSq num- in result `seq` return result-
− imaginary/SumSquares/dph/SumSquaresVector.hs
@@ -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-
− imaginary/SumSquares/dph/SumSquaresVectorised.hs
@@ -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)))
lib/Points2D/Generate.hs view
@@ -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
lib/Points2D/Types.hs view
@@ -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)
+ lib/System/Console/ParseArgs.hs view
@@ -0,0 +1,667 @@+{-# LANGUAGE FlexibleInstances, DeriveDataTypeable #-}+-- Full-featured argument parsing library for Haskell programs+-- Bart Massey <bart@cs.pdx.edu>++-- 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)
lib/Timing.hs view
@@ -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)
− real/NBody/Batch/Config.hs
@@ -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 }
− real/NBody/Batch/MainArgs.hs
@@ -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" }-- ]
− real/NBody/Common/Body.hs
@@ -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)
− real/NBody/Common/Dump.hs
@@ -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
− real/NBody/Common/Util.hs
@@ -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
− real/NBody/Common/World.hs
@@ -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' }--
− real/NBody/Gloss/Config.hs
@@ -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 <window size :: Int>"- else x -- in checkMode $- Config- { configWindowSize = mWindowSize- , configShouldDrawTree = shouldDrawTree- , configRate = rate- , configSolverName = solverName- , configBodyCount = bodyCount- , configBodyMass = bodyMass- , configTimeStep = timeStep- , configEpsilon = epsilon- , configStartDiscSize = discSize- , configStartSpeed = startSpeed- , configMaxSteps = mMaxSteps - , configDumpFinal = mFilePath }
− real/NBody/Gloss/Draw.hs
@@ -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
− real/NBody/Gloss/MainArgs.hs
@@ -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" }-- ]
− real/NBody/MainBatch.hs
@@ -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-
− real/NBody/MainGloss.hs
@@ -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--
− real/NBody/Solver.hs
@@ -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'
− real/NBody/Solver/ListBH/Solver.hs
@@ -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-
− real/NBody/Solver/NestedBH/Solver.hs
@@ -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-
− real/NBody/Solver/VectorBH/Solver.hs
@@ -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--
− real/NBody/Solver/VectorNaive/Solver.hs
@@ -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)
− spectral/QuickHull/dph/Main.hs
@@ -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 <points> [out.svg]" ]----- | Run the benchmark.-run :: Int -- ^ How many points to use.- -> Maybe String -- ^ File name to dump an SVG of the output to.- -> IO ()- -run pointCount mFileSVG- = do- vPoints <- pointsPArrayOfUArray- $ genPointsDisc pointCount (400, 400) 350 -- -- Force points to create the input vector.- vPoints `seq` return ()-- -- Compute the convex hull.- (vHull, tElapsed)- <- time - $ let vHull = quickhullPA vPoints- in vHull `seq` return vHull- - -- Print how long it took.- putStr $ prettyTime tElapsed- - -- If we were asked for an SVG then write it out to file.- maybe (return ())- (\fileSVG -> - writeFile fileSVG- $ makeSVG (roundPoints $ P.toList vPoints) - (roundPoints $ P.toList vHull))- mFileSVG
− spectral/QuickHull/dph/QuickHullVect.hs
@@ -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))-
− spectral/QuickHull/lib/SVG.hs
@@ -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- $ [ "<svg width=\"100%\" height=\"100%\" version=\"1.1\" xmlns=\"http://www.w3.org/2000/svg\">" ]- ++ [svgPolygon hull]- ++ map svgPoint points- ++ map svgPointHull hull- ++ ["</svg>"]--svgPolygon :: [(Int, Int)] -> String-svgPolygon points- = "<polygon"- ++ " points=\"" ++ (concat [show x ++ "," ++ show y ++ " " | (x, y) <- points]) ++ "\"" - ++ " style=\"fill:#d0d0ff;stroke:#000000;stroke-width:1\""- ++ "/>"--svgPoint :: (Int, Int) -> String-svgPoint (x, y)- = "<circle cx=\"" ++ show x ++ "\" cy=\"" ++ show y ++ "\" r=\"0.5\""- ++ " style=\"stroke:#000000\""- ++ "/>"--svgPointHull :: (Int, Int) -> String-svgPointHull (x, y)- = "<circle cx=\"" ++ show x ++ "\" cy=\"" ++ show y ++ "\" r=\"1\""- ++ " style=\"stroke:#ff0000\""- ++ "/>"- -roundPoints :: [(Double, Double)] -> [(Int, Int)]-roundPoints ps = [(round x, round y) | (x, y) <- ps]
− spectral/QuickHull/vector/Main.hs
@@ -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 <alg> <points> [out.svg]"- , " algs: " ++ (show $ map fst algs) ++ "\n" ]--run :: ( Vector Point -> IO (Vector Point)- , Int- , Maybe FilePath)- -> IO ()- -run (fun, pointCount, mFileSVG) - = do- let vPoints = genPointsDisc pointCount (400, 400) 350 -- -- Force points to create the input vector.- V.force vPoints `seq` return ()-- -- Compute the convex hull.- (vHull, tElapsed)- <- time- $ do vHull <- fun vPoints- vHull `seq` return vHull- - -- Print how long it took - putStr $ prettyTime tElapsed-- -- If we were asked for an SVG then write it out to file.- maybe (return ())- (\fileSVG -> - writeFile fileSVG- $ makeSVG (roundPoints $ V.toList vPoints) - (roundPoints $ V.toList vHull))- mFileSVG-
− spectral/QuickHull/vector/QuickHullIO.hs
@@ -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) -
− spectral/QuickHull/vector/QuickHullSplit.hs
@@ -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)
− spectral/QuickHull/vector/QuickHullVector.hs
@@ -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)-
− spectral/QuickSort/dph/Main.hs
@@ -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 <length>" ]----- | Run the benchmark.-run :: Int -> IO ()-run len- = do -- Create the input vector- let vInts = fromUArrPA'- $ randomishDoubles len 0 1 1234-- vInts `seq` return ()-- -- Compute the convex hull.- (vSorted, tElapsed)- <- time - $ let vSorted = quicksortPA vInts- in vSorted `seq` return vSorted- - -- Print how long it took.- putStr $ prettyTime tElapsed- - -- Check they're really sorted.- print $ isSorted $ P.toUArrPA vSorted- ---- | Check if a vector is sorted (monotonically increasing)-isSorted :: Vector Double -> Bool-isSorted = isJust . V.fold1M (\x y -> if y >= x then Just y else Nothing)
− spectral/QuickSort/dph/QuickSortVect.hs
@@ -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)
− spectral/SMVM/dph/Main.hs
@@ -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 <file>" ]--run fileName- = do (matrix, vector) <- loadPA fileName-- matrix `seq` return ()- vector `seq` return ()-- -- Multiply sparse matrix by the dense vector.- (vResult, tElapsed)- <- time $ let result = smvmPA matrix vector- in P.nf result `seq` return result- - -- Print how long it took.- putStr $ prettyTime tElapsed-- -- Print some info about the test setup.- putStrLn $ "vector length = " ++ show (U.length (toUArrPA vector))--- putStrLn $ "matrix height = " ++ show (U.length (toUArrPA matrix))- - - -- Print checksum of resulting vector.- putStrLn $ "result sum = " ++ show (U.sum (toUArrPA vResult))------ | Load a test file containing a sparse matrix and dense vector.-loadPA :: String -- ^ filename.- -> IO ( PArray (PArray (Int, Double)) -- sparse matrix- , PArray Double) -- dense vector--loadPA fileName- = do (segd, arrMatrixElems, arrVector) <- loadUArr fileName-- let paMatrix = nestUSegdPA' segd (fromUArrPA_2' arrMatrixElems)- let paVector = fromUArrPA' arrVector- return (paMatrix, paVector)----- | Load a test file containing a sparse matrix and dense vector.-loadUArr :: String -- ^ filename- -> IO ( U.Segd -- segment descriptor saying what array elements- -- belong to each row of the matrix.- , U.Array (Int, Double) -- column indices and matrix elements- , U.Array Double) -- the dense vector--loadUArr fname - = do h <- openBinaryFile fname ReadMode-- -- check magic numbers at start of file to guard against word-size screwups.- alloca $ \ptr -> do- hGetBuf h ptr (sizeOf (undefined :: Int))- magic1 :: Int <- peek ptr- hGetBuf h ptr (sizeOf (undefined :: Int))- magic2 :: Int <- peek ptr- if magic1 == 0xc0ffee00 Prelude.&& magic2 == 0x12345678 - then return ()- else error $ "bad magic in " ++ fname-- -- number of elements in each row of the matrix.- lengths <- U.hGet h-- -- indices of all the elements.- indices <- U.hGet h-- -- values of the matrix elements.- values <- U.hGet h-- -- the dense vector.- vector <- U.hGet h-- evaluate lengths- evaluate indices- evaluate values- evaluate vector-- let segd = U.lengthsToSegd lengths- matrix = U.zip indices values-- return (segd, matrix, vector)
− spectral/SMVM/dph/SMVMVectorised.hs
@@ -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 :]-