cl3-3.0.0.0: benchmarks/NbodyGameCl3.hs
{-# LANGUAGE ViewPatterns #-}
{-# LANGUAGE BangPatterns #-}
{-# LANGUAGE CPP #-}
{-# LANGUAGE ScopedTypeVariables #-}
{-# LANGUAGE InstanceSigs #-}
#if __GLASGOW_HASKELL__ == 810
-- Work around to fix GHC Issue #15304, issue popped up again in GHC 8.10, it should be fixed in GHC 8.12
-- This code is meant to reproduce MR 2608 for GHC 8.10
{-# OPTIONS_GHC -funfolding-keeness-factor=1 -funfolding-use-threshold=80 #-}
#endif
------------------------------------------------------------------
-- |
-- Description: Benchmark for Algebra.Geometric.Cl3
-- Uses standard "The Computer language Benchmarks Game" implementation
-- with Criterion to compare the performance difference between
-- using the Cl3 library vs. the benchmarks implementation using doubles.
--
--
--
-- This derivative work was derived from the excellent work of:
-- Branimir Maksimovic.
--
-------------------------------------------------------------------
module Main (main) where
import Criterion.Main (defaultMain, bgroup, bench, nfIO) -- To add Criterion to the benchmark
import Algebra.Geometric.Cl3 (Cl3(..), toR, toV3, toCl3_V3, toCl3_R, fromCl3_V3, fromCl3_R, Cl3_V3(), Cl3_R()) -- To add Cl3
import Foreign.Ptr (Ptr, castPtr, plusPtr)
import Foreign.Storable (Storable(..))
import Foreign.Marshal.Alloc (mallocBytes)
import Control.Monad (when)
import Text.Printf (printf)
-------------------------------------------------------------------
-- |
-- An n-body solver for the Sun and gas giants using the semi-implicit Euler method.
-- The benchmark uses the Criterion tool to measure timing.
-- The benchmark advances 50 million steps.
--
-- The @nbodyBaseline@ is derived from Branimir Maksimovic's work with some changes suggested by @hlint@.
-- @nbodyCl3@ aggressively uses the "to" functions like 'toR' and 'toV3' to take advantage of the case of
-- known constructor optimization. @nbodyAPS@ does not make aggressive use "to" functions. The
-- program would not finish compiling (memory would grow to 54GB while compiling) before INLINE pragmas
-- were added to the Cl3 library.
--
-- * Baseline 20171007:
--
-- Initial Commit of the Code, GHC 8.0.2
--
-- >>> ./benchmarks/NbodyGameCl3
-- benchmarking nbodyBaseline/50000000
-- time 37.99 s (37.92 s .. 38.08 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 37.93 s (37.89 s .. 37.95 s)
-- std dev 36.25 ms (0.0 s .. 38.29 ms)
-- variance introduced by outliers: 19% (moderately inflated)
-- benchmarking nbodyCl3/50000000
-- time 38.95 s (38.91 s .. 39.03 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 38.95 s (38.94 s .. 38.96 s)
-- std dev 14.70 ms (0.0 s .. 16.66 ms)
-- variance introduced by outliers: 19% (moderately inflated)
-- benchmarking nbodyAPS/50000000
-- time 179.4 s (177.4 s .. 181.4 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 179.3 s (178.8 s .. 179.5 s)
-- std dev 395.0 ms (0.0 s .. 413.9 ms)
-- variance introduced by outliers: 19% (moderately inflated)
--
-- * 20171016:
--
-- Change of INLINE phase control to fix ghc simplifier ticks issue, GHC 8.0.2
--
-- >>> ./benchmarks/NbodyGameCl3
-- benchmarking nbodyBaseline/50000000
-- time 37.21 s (36.80 s .. 37.46 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 37.02 s (36.89 s .. 37.10 s)
-- std dev 125.4 ms (0.0 s .. 144.2 ms)
-- variance introduced by outliers: 19% (moderately inflated)
-- benchmarking nbodyCl3/50000000
-- time 38.57 s (38.07 s .. 39.20 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 38.55 s (38.44 s .. 38.64 s)
-- std dev 128.8 ms (0.0 s .. 147.5 ms)
-- variance introduced by outliers: 19% (moderately inflated)
-- benchmarking nbodyAPS/50000000
-- time 176.9 s (169.3 s .. 181.7 s)
-- 1.000 R² (0.999 R² .. 1.000 R²)
-- mean 177.1 s (176.1 s .. 177.9 s)
-- std dev 1.180 s (0.0 s .. 1.316 s)
-- variance introduced by outliers: 19% (moderately inflated)
--
-- * 20171130:
--
-- Though, thru, some unknown reason, an extreme increase of performance has occurred.
-- Not sure what changed but a huge improvement, dev-lang/ghc-8.0.2,
-- upgraded sys-devel/gcc-6.4.0, sys-devel/llvm-3.9.1-r1, linux-4.12.12-gentoo?
-- Spooky thing is that the parent compiled executable from before that ran ~35.6 Sec
-- is now ~10.7 Sec.
--
-- >>> ./benchmarks/NbodyGameCl3
-- benchmarking nbodyBaseline/50000000
-- time 14.04 s (13.67 s .. 14.60 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 13.74 s (13.59 s .. 13.85 s)
-- std dev 168.6 ms (0.0 s .. 190.4 ms)
-- variance introduced by outliers: 19% (moderately inflated)
-- benchmarking nbodyCl3/50000000
-- time 15.68 s (14.80 s .. 16.12 s)
-- 1.000 R² (0.999 R² .. 1.000 R²)
-- mean 15.79 s (15.64 s .. 15.88 s)
-- std dev 141.7 ms (0.0 s .. 161.5 ms)
-- variance introduced by outliers: 19% (moderately inflated)
-- benchmarking nbodyAPS/50000000
-- time 182.9 s (182.4 s .. 183.5 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 182.7 s (182.5 s .. 182.8 s)
-- std dev 162.1 ms (0.0 s .. 183.2 ms)
-- variance introduced by outliers: 19% (moderately inflated)
--
-- 20180309 :
--
-- Well, while correcting spelling I found several issues re: Ord;
-- Ord is now more what you'd expect (now -2 is smaller than -1).
-- It appears to have sped up the performance a bit.
--
-- benchmarking nbodyBaseline/50000000
-- time 13.47 s (12.95 s .. 13.89 s)
-- 1.000 R² (1.000 R² .. NaN R²)
-- mean 13.54 s (13.42 s .. 13.59 s)
-- std dev 98.11 ms (0.0 s .. 103.4 ms)
-- variance introduced by outliers: 19% (moderately inflated)
-- benchmarking nbodyCl3/50000000
-- time 15.50 s (15.45 s .. 15.59 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 15.54 s (15.52 s .. 15.56 s)
-- std dev 27.35 ms (0.0 s .. 30.00 ms)
-- variance introduced by outliers: 19% (moderately inflated)
-- benchmarking nbodyAPS/50000000
-- time 154.5 s (149.7 s .. 159.5 s)
-- 1.000 R² (0.999 R² .. 1.000 R²)
-- mean 153.9 s (152.8 s .. 154.6 s)
-- std dev 1.119 s (0.0 s .. 1.283 s)
-- variance introduced by outliers: 19% (moderately inflated)
--
-- 20180518:
--
-- Submitted the library to review on Haskell :: Reddit, implemented
-- several suggestions, the main one that may affect performance
-- would be removing the orphan instances and moving the instances
-- to Cl3 and using CPP gating of the optional instances.
-- Performance seems to stay the same.
--
-- benchmarking nbodyBaseline/50000000
-- time 13.38 s (12.37 s .. 14.02 s)
-- 0.999 R² (0.998 R² .. 1.000 R²)
-- mean 13.35 s (13.23 s .. 13.52 s)
-- std dev 167.1 ms (45.82 ms .. 225.3 ms)
-- variance introduced by outliers: 19% (moderately inflated)
--
-- benchmarking nbodyCl3/50000000
-- time 15.67 s (15.58 s .. 15.72 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 15.73 s (15.70 s .. 15.79 s)
-- std dev 60.02 ms (583.2 μs .. 70.85 ms)
-- variance introduced by outliers: 19% (moderately inflated)
--
-- benchmarking nbodyAPS/50000000
-- time 156.2 s (154.5 s .. 159.3 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 155.1 s (154.8 s .. 155.7 s)
-- std dev 560.8 ms (88.83 ms .. 714.2 ms)
-- variance introduced by outliers: 19% (moderately inflatemarking nbodyBaseline/50000000
--
-- 20180816: cl3-1.0.0.3
--
-- Modified code to use ghc-8.4.2, ghc greater than 8.0 require -fno-worker-wrapper
-- to compile within reasonable time and space constraints.
--
-- time 9.271 s (9.231 s .. 9.344 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 9.495 s (9.407 s .. 9.636 s)
-- std dev 133.9 ms (53.72 ms .. 174.3 ms)
-- variance introduced by outliers: 19% (moderately inflated)
--
-- benchmarking nbodyCl3/50000000
-- time 10.10 s (9.797 s .. 10.30 s)
-- 1.000 R² (NaN R² .. 1.000 R²)
-- mean 10.10 s (10.05 s .. 10.16 s)
-- std dev 58.98 ms (28.18 ms .. 72.96 ms)
-- variance introduced by outliers: 19% (moderately inflated)
--
-- benchmarking nbodyAPS/50000000
-- time 178.6 s (NaN s .. 181.1 s)
-- 1.000 R² (1.000 R² .. 1.000 R²)
-- mean 177.2 s (175.5 s .. 178.2 s)
-- std dev 1.670 s (618.8 ms .. 2.283 s)
-- variance introduced by outliers: 19% (moderately inflated)
--
-------------------------------------------------------------------
main :: IO ()
main = do n <- return (50000000 :: Int)
defaultMain [bgroup "nbodyBaseline" [ bench (show n) $ nfIO (nbodyBaseline n)]
,bgroup "nbodyCl3" [ bench (show n) $ nfIO (nbodyCl3 n)]
,bgroup "nbodyAPS" [ bench (show n) $ nfIO (nbodyAPS n)]
,bgroup "nbodyCl3_V3" [ bench (show n) $ nfIO (nbodyCl3_V3 n)]
]
nbodyBaseline :: (Num t, Eq t) => t -> IO ()
nbodyBaseline n = do
pPlanets <- fromList planets -- load planets into a stack
nbodyInit pPlanets -- initialize the system momentum
energy pPlanets >>= print -- printf "%.9f\n" -- calculate and print out the initial energy of the system
run n pPlanets -- solve the initial value problem, increment the system 'n' times
energy pPlanets >>= print -- printf "%.9f\n" -- calculate and print out the final energy of the system
nbodyCl3 :: (Eq t, Num t) => t -> IO ()
nbodyCl3 n = do
pPlanetsCl3 <- fromListCl3 planetsCl3
nbodyInitCl3 pPlanetsCl3
energyCl3 pPlanetsCl3 >>= print.toR -- printf "%0.9f\n"
runCl3 n pPlanetsCl3
energyCl3 pPlanetsCl3 >>= print.toR -- printf "%0.9f\n"
nbodyAPS :: (Eq t, Num t) => t -> IO ()
nbodyAPS n = do
pPlanetsAPS <- fromListCl3 planetsCl3
nbodyInitAPS pPlanetsAPS
energyAPS pPlanetsAPS >>= print.toR -- printf "%0.9f\n"
runAPS n pPlanetsAPS
energyAPS pPlanetsAPS >>= print.toR -- printf "%0.9f\n"
nbodyCl3_V3 :: (Eq t, Num t) => t -> IO ()
nbodyCl3_V3 n = do
pPlanetsCl3_V3 <- fromListCl3_V3 planetsCl3_V3
nbodyInitCl3_V3 pPlanetsCl3_V3
energyCl3_V3 pPlanetsCl3_V3 >>= print
runCl3_V3 n pPlanetsCl3_V3
energyCl3_V3 pPlanetsCl3_V3 >>= print
-- | 'run' advances the system 'n' times by recursively calling 'run'
run :: (Eq t, Num t) => t -> Ptr Planet -> IO ()
run 0 _ = return ()
run i p = do
advance p
run (i-1) p
runCl3 :: (Eq t, Num t) => t -> Ptr PlanetCl3 -> IO ()
runCl3 0 _ = return ()
runCl3 i p = do
advanceCl3 p
runCl3 (i-1) p
runAPS :: (Eq t, Num t) => t -> Ptr PlanetCl3 -> IO ()
runAPS 0 _ = return ()
runAPS i p = do
advanceAPS p
runAPS (i-1) p
-- runCl3_V3 is a tail recursive action to advance the system
runCl3_V3 :: (Eq t, Num t) => t -> Ptr PlanetCl3_V3 -> IO ()
runCl3_V3 0 _ = return ()
runCl3_V3 i p = do
advanceCl3_V3 p
runCl3_V3 (i-1) p
-- | 'offsetMomentum' is used to update the sun's momentum
-- so that the total momentum will be zero for the whole
-- system at initialization.
offsetMomentum :: Planet -> (Double, Double, Double) -> Planet
offsetMomentum p (px,py,pz) = p {vx = -px / solarMass,
vy = -py / solarMass,
vz = -pz / solarMass
}
offsetMomentumCl3 :: PlanetCl3 -> Cl3 -> PlanetCl3
offsetMomentumCl3 (PlanetCl3 (pos') _ (mass')) (momentum) = PlanetCl3 pos' (toV3.negate $ toV3 momentum * (toR.recip.toR $ mass')) mass'
offsetMomentumAPS :: PlanetCl3 -> Cl3 -> PlanetCl3
offsetMomentumAPS p momCl3 = p{velCl3 = negate momCl3 / massCl3 p}
offsetMomentumCl3_V3 :: PlanetCl3_V3 -> Cl3 -> PlanetCl3_V3
offsetMomentumCl3_V3 p momCl3 = p{velCl3_V3 = toCl3_V3 (negate momCl3 / (fromCl3_R.massCl3_R $ p))}
-- | 'nbodyInit' initializes the velocity of the sun so that
-- the total system momentum is zero.
nbodyInit :: Ptr Planet -> IO ()
nbodyInit pPlanets = do
let initBl (px,py,pz) i = -- do "Is this do do?" -- init cycles through the stack and builds up the total momentum
if i < length planets
then do
p <- peekElemOff pPlanets i
initBl (px + vx p * mass p, py + vy p * mass p, pz + vz p * mass p) (i+1) -- requires velocity and mass
else return (px,py,pz)
s <- initBl (0,0,0) 0 -- total momentum of the system
p <- peek pPlanets -- grab the sun, don't get burned
poke pPlanets $ offsetMomentum p s -- update the sun's momentum
-- end of nbodyInit
nbodyInitCl3 :: Ptr PlanetCl3 -> IO ()
nbodyInitCl3 pPlanetsCl3 = do
let initCl3 (toV3 -> pCl3) i =
if i < length planetsCl3
then do
(PlanetCl3 _ (toV3 -> vel') (toR -> mass')) <- peekElemOff pPlanetsCl3 i
initCl3 (toV3 $! pCl3 + vel' * mass') (i+1)
else return pCl3
totalMomentum <- initCl3 (V3 0 0 0) 0
sun0 <- peek pPlanetsCl3
poke pPlanetsCl3 $ offsetMomentumCl3 sun0 totalMomentum
-- end of nbodyInitCl3
nbodyInitAPS :: Ptr PlanetCl3 -> IO ()
nbodyInitAPS pPlanetsCl3 = do
let initCl3 !pCl3 i =
if i < length planetsCl3
then do
p <- peekElemOff pPlanetsCl3 i
initCl3 (pCl3 + velCl3 p * massCl3 p) (i+1)
else return pCl3
totalMomentum <- initCl3 (APS 0 0 0 0 0 0 0 0) 0
sun0 <- peek pPlanetsCl3
poke pPlanetsCl3 $ offsetMomentumAPS sun0 totalMomentum
-- end of nbodyInitAPS
nbodyInitCl3_V3 :: Ptr PlanetCl3_V3 -> IO ()
nbodyInitCl3_V3 pPlanetsCl3_V3 = do
let initCl3_V3 !pCl3_V3 i = -- initCl3_V3 takes a momentum accumulator and an iteration count
if i < length planetsCl3_V3 -- if the iteration count is less than the length of the list of bodies; add more momentum
then do
p <- peekElemOff pPlanetsCl3_V3 i -- pull off a body at the current iteration count
initCl3_V3 (pCl3_V3 + (fromCl3_V3.velCl3_V3 $ p) * (fromCl3_R.massCl3_R $ p)) (i+1) -- accumulate the momentum in the accumulator, increase the iteration count and recursively call the function again
else return pCl3_V3 -- no more body's momentum left to accumulate, so return the total momentum of all of the bodies
totalMomentum <- initCl3_V3 (V3 0 0 0) 0 -- calculate the total momentum, start with a zero momentum vector in the accumulator, and an iteration count of zero, returns a Cl3
sun0 <- peek pPlanetsCl3_V3 -- grab the Sun, don't get burned
poke pPlanetsCl3_V3 $ offsetMomentumCl3_V3 sun0 totalMomentum -- update the Sun's velocity with the total momentum, such that the system has zero total momentum
-- | 'squared' scalar product of a vector with its self
squared :: Num a => a -> a -> a -> a
squared x' y' z' = x' * x' + y' * y' + z' * z'
-- | 'energy' calculate the total energy of the system
-- this is the gravitational potential energy added to
-- the kinetic energy of all of the planets.
energy :: Ptr Planet -> IO Double
energy pPlanets = do
let energy' e i = if i < length planets -- ''
then do
p <- peekElemOff pPlanets i
e1 <- energy'' p (i+1) e
e2 <- energy' e (i+1)
return $ e + 0.5 * mass p * squared (vx p) (vy p) (vz p) + e1 + e2 -- requires mass and velocity
else return e
energy'' p j e = if j < length planets
then do
pj <- peekElemOff pPlanets j
let distance = sqrt $ squared dx dy dz
dx = x pj - x p
dy = y pj - y p
dz = z pj - z p
e1 <- energy'' p (j+1) e
return $ e - (mass p * mass pj) / distance + e1 -- requires mass and position
else return e
energy' 0.0 0 -- starts off the recursive calculation of energy
energyCl3 :: Ptr PlanetCl3 -> IO Cl3
energyCl3 pPlanetsCl3 = do
let energy' e i =
if i < length planetsCl3
then do
p@(PlanetCl3 _ (toV3 -> vel1) (toR -> mass1)) <- peekElemOff pPlanetsCl3 i
let !ke = toR $ 0.5 * mass1 * toR (vel1^(2 :: Int))
e1 <- energy'' p (i+1) e
e2 <- energy' e (i+1)
return $! e + ke + e1 + e2
else return e
energy'' p@(PlanetCl3 (toV3 -> pos1) _ (toR -> mass1)) j e =
if j < length planetsCl3
then do
(PlanetCl3 (toV3 -> pos2) _ (toR -> mass2)) <- peekElemOff pPlanetsCl3 j
let !dpos = toV3 $ pos2 - pos1
!distance = abs dpos
!pe = toR $ (mass1 * mass2) / distance
e1 <- energy'' p (j+1) e
return $! e - pe + e1 -- requires mass and position
else return e
energy' (R 0) 0 -- starts off the recursive calculation of energy
-- end of energyCl3
energyAPS :: Ptr PlanetCl3 -> IO Cl3
energyAPS pPlanetsCl3 = do
let energy' e i =
if i < length planetsCl3
then do
p <- peekElemOff pPlanetsCl3 i
e1 <- energy'' p (i+1) e
e2 <- energy' e (i+1)
return $! e + 0.5 * massCl3 p * (velCl3 p)^(2 :: Int) + e1 + e2
else return e
energy'' p j e =
if j < length planetsCl3
then do
pj <- peekElemOff pPlanetsCl3 j
let !distance = abs dpos
!dpos = posCl3 pj - posCl3 p
e1 <- energy'' p (j+1) e
return $! e - (massCl3 p * massCl3 pj) / distance + e1
else return e
energy' (APS 0 0 0 0 0 0 0 0) 0 -- starts off the recursive calculation of energy
-- end of energyAPS
energyCl3_V3 :: Ptr PlanetCl3_V3 -> IO Cl3
energyCl3_V3 pPlanetsCl3_V3 = do
let energy' e i =
if i < length planetsCl3_V3
then do
p <- peekElemOff pPlanetsCl3_V3 i
e1 <- energy'' p (i+1) e
e2 <- energy' e (i+1)
return $! e + 0.5 * (fromCl3_R.massCl3_R $ p) * (toR $ (fromCl3_V3.velCl3_V3 $ p)^(2 :: Int)) + e1 + e2 -- the Kenitic Energy and other accumulated energies
else return e -- all of the bodies energies have been calculated, return the total energy
energy'' p j e = -- a Potential Energy calculation
if j < length planetsCl3_V3
then do
pj <- peekElemOff pPlanetsCl3_V3 j
let !distance = abs dpos
!dpos = (fromCl3_V3.posCl3_V3 $ pj) - (fromCl3_V3.posCl3_V3 $ p)
e1 <- energy'' p (j+1) e
return $! e - (((fromCl3_R.massCl3_R $ p) * (fromCl3_R.massCl3_R $ pj)) / distance) + e1
else return e
energy' (R 0) 0 -- starts off the recursive calculation of energy, returns the total energy as a Cl3
-- | 'advance' integrates the system of differential equations using a
-- semi-implicit Euler method, also called symplectic Euler method.
-- The first order method updates the velocities based on the acceleration
-- calculation, then it updates the positions. The velocity verlet
-- would be an improvement in accuracy because it is a second order
-- symplectic integrator.
advance :: Ptr Planet -> IO ()
advance pPlanets = do
let advance' i = when (i < length planets) $ do -- loops through all the planets and updates the velocity
let loop j = when (j < length planets) $ do
ii <- peekElemOff pPlanets i
jj <- peekElemOff pPlanets j
let mag = dt / (dSquared * sqrt dSquared)
dSquared = squared dx dy dz
dx = x ii - x jj
dy = y ii - y jj
dz = z ii - z jj
pokeV pPlanets i ii{ -- requires position and mass to update the velocity
vx = vx ii - dx * mass jj * mag,
vy = vy ii - dy * mass jj * mag,
vz = vz ii - dz * mass jj * mag
}
pokeV pPlanets j jj{ -- requires position and mass to update the velocity
vx = vx jj + dx * mass ii * mag,
vy = vy jj + dy * mass ii * mag,
vz = vz jj + dz * mass ii * mag
}
loop (j+1)
loop (i+1)
advance' (i+1)
advance'' i = when (i < length planets) $ do -- loops through all of the planets and updates the position
p <- peekElemOff pPlanets i
pokeC pPlanets i p{ -- requires position and velocity to update the position
x = x p + dt * vx p,
y = y p + dt * vy p,
z = z p + dt * vz p
}
advance'' (i+1)
advance' 0 -- update all of the planets velocities
advance'' 0 -- update all of the planets positions
-- end of advance
advanceCl3 :: Ptr PlanetCl3 -> IO ()
advanceCl3 pPlanetsCl3 = do
let advance' i = when (i < length planetsCl3) $ do -- loops through all the planets and updates the velocity
let loop j = when (j < length planetsCl3) $ do
(PlanetCl3 (toV3 -> posi) (toV3 -> veli) (toR -> massi)) <- peekElemOff pPlanetsCl3 i
(PlanetCl3 (toV3 -> posj) (toV3 -> velj) (toR -> massj)) <- peekElemOff pPlanetsCl3 j
let !dpos = posi - posj
!(R dMag) = abs dpos
!mag = R (dt / (dMag^(3 :: Int)))
!veli' = veli - dpos * massj * mag
!velj' = velj + dpos * massi * mag
pokeVCl3 pPlanetsCl3 i (PlanetCl3 posi veli' massi)
pokeVCl3 pPlanetsCl3 j (PlanetCl3 posj velj' massj)
loop (j+1)
loop (i+1)
advance' (i+1)
advance'' i = when (i < length planetsCl3) $ do -- loops through all of the planets and updates the position
(PlanetCl3 (toV3 -> posi) (toV3 -> veli) massi) <- peekElemOff pPlanetsCl3 i
let !posi' = posi + R dt * veli
pokeCCl3 pPlanetsCl3 i (PlanetCl3 posi' veli massi)
advance'' (i+1)
advance' 0 -- update all of the planets velocities
advance'' 0 -- update all of the planets positions
-- end of advanceCl3
advanceAPS :: Ptr PlanetCl3 -> IO ()
advanceAPS pPlanetsCl3 = do
let advance' i = when (i < length planetsCl3) $ do -- loops through all the planets and updates the velocity
let loop j = when (j < length planetsCl3) $ do
ii <- peekElemOff pPlanetsCl3 i
jj <- peekElemOff pPlanetsCl3 j
let !mag = R dt / (dSquared * abs dpos)
!dSquared = dpos * dpos
!dpos = posCl3 ii - posCl3 jj
pokeVCl3 pPlanetsCl3 i ii{velCl3 = velCl3 ii - dpos * massCl3 jj * mag}
pokeVCl3 pPlanetsCl3 j jj{velCl3 = velCl3 jj + dpos * massCl3 ii * mag}
loop (j+1)
loop (i+1)
advance' (i+1)
advance'' i = when (i < length planetsCl3) $ do -- loops through all of the planets and updates the position
p <- peekElemOff pPlanetsCl3 i
pokeCCl3 pPlanetsCl3 i p{posCl3 = posCl3 p + R dt * velCl3 p}
advance'' (i+1)
advance' 0 -- update all of the planets velocities
advance'' 0 -- update all of the planets positions
-- end of advanceAPS
advanceCl3_V3 :: Ptr PlanetCl3_V3 -> IO ()
advanceCl3_V3 pPlanetsCl3_V3 = do
let advance' i = when (i < length planetsCl3_V3) $ do -- loops through all the planets and updates the velocity
let loop j = when (j < length planetsCl3_V3) $ do
ii <- peekElemOff pPlanetsCl3_V3 i
jj <- peekElemOff pPlanetsCl3_V3 j
let -- !mag = R dt / (toR dSquared * abs dpos)
-- !dSquared = dpos * dpos
!dpos = (fromCl3_V3.posCl3_V3 $ ii) - (fromCl3_V3.posCl3_V3 $ jj)
!(R !dMag) = abs dpos
!mag = R (dt / (dMag^(3 :: Int)))
pokeVCl3_V3 pPlanetsCl3_V3 i ii{velCl3_V3 = toCl3_V3 $ (fromCl3_V3.velCl3_V3 $ ii) - dpos * (fromCl3_R.massCl3_R $ jj) * mag}
pokeVCl3_V3 pPlanetsCl3_V3 j jj{velCl3_V3 = toCl3_V3 $ (fromCl3_V3.velCl3_V3 $ jj) + dpos * (fromCl3_R.massCl3_R $ ii) * mag}
loop (j+1)
loop (i+1)
advance' (i+1)
advance'' i = when (i < length planetsCl3_V3) $ do -- loops through all of the planets and updates the position
p <- peekElemOff pPlanetsCl3_V3 i
pokeCCl3_V3 pPlanetsCl3_V3 i p{posCl3_V3 = toCl3_V3 $ (fromCl3_V3.posCl3_V3 $ p) + R dt * (fromCl3_V3.velCl3_V3 $ p)}
advance'' (i+1)
advance' 0 -- update all of the planets velocities
advance'' 0 -- update all of the planets positions
data Planet = Planet {x :: !Double,
y :: !Double,
z :: !Double,
vx :: !Double,
vy :: !Double,
vz :: !Double,
mass :: !Double } deriving (Show)
planets :: [Planet]
planets = [sun, jupiter, saturn, uranus, neptune]
-- | A planet can be described by three cliffors.
data PlanetCl3 = PlanetCl3 {posCl3 :: !Cl3,
velCl3 :: !Cl3,
massCl3 :: !Cl3} deriving (Show)
planetsCl3 :: [PlanetCl3]
planetsCl3 = [PlanetCl3 sunPos sunVel sunMass,
PlanetCl3 jupiterPos jupiterVel jupiterMass,
PlanetCl3 saturnPos saturnVel saturnMass,
PlanetCl3 uranusPos uranusVel uranusMass,
PlanetCl3 neptunePos neptuneVel neptuneMass]
data PlanetCl3_V3 = PlanetCl3_V3 {posCl3_V3 :: !Cl3_V3
,velCl3_V3 :: !Cl3_V3
,massCl3_R :: !Cl3_R}
instance Show PlanetCl3_V3 where
show p = "PlanetCl3_V3 " ++ show (show.fromCl3_V3.posCl3_V3 $ p) ++ show (show.fromCl3_V3.velCl3_V3 $ p) ++ show (show.fromCl3_R.massCl3_R $ p)
planetsCl3_V3 :: [PlanetCl3_V3]
planetsCl3_V3 = [PlanetCl3_V3 (toCl3_V3 sunPos) (toCl3_V3 sunVel) (toCl3_R sunMass)
,PlanetCl3_V3 (toCl3_V3 jupiterPos) (toCl3_V3 jupiterVel) (toCl3_R jupiterMass)
,PlanetCl3_V3 (toCl3_V3 saturnPos) (toCl3_V3 saturnVel) (toCl3_R saturnMass)
,PlanetCl3_V3 (toCl3_V3 uranusPos) (toCl3_V3 uranusVel) (toCl3_R uranusMass)
,PlanetCl3_V3 (toCl3_V3 neptunePos) (toCl3_V3 neptuneVel) (toCl3_R neptuneMass)]
sun :: Planet
sun = Planet {x = 0, y = 0, z = 0,
vx = 0, vy = 0, vz = 0,
mass = solarMass
}
sunPos :: Cl3
sunPos = V3 0 0 0
sunVel :: Cl3
sunVel = V3 0 0 0
sunMass :: Cl3
sunMass = R solarMass
jupiter :: Planet
jupiter = Planet
{x = 4.84143144246472090e+00, y = -1.16032004402742839e+00, z= -1.03622044471123109e-01,
vx = 1.66007664274403694e-03*dp, vy = 7.69901118419740425e-03*dp, vz = -6.90460016972063023e-05*dp,
mass = 9.54791938424326609e-04 * solarMass
}
jupiterPos :: Cl3
jupiterPos = V3 (4.84143144246472090e+00) (-1.16032004402742839e+00) (-1.03622044471123109e-01)
jupiterVel :: Cl3
jupiterVel = V3 (1.66007664274403694e-03*dp) (7.69901118419740425e-03*dp) (-6.90460016972063023e-05*dp)
jupiterMass :: Cl3
jupiterMass = R (9.54791938424326609e-04 * solarMass)
saturn :: Planet
saturn = Planet
{x = 8.34336671824457987e+00, y = 4.12479856412430479e+00, z = -4.03523417114321381e-01,
vx = -2.76742510726862411e-03*dp, vy = 4.99852801234917238e-03*dp, vz = 2.30417297573763929e-05*dp,
mass = 2.85885980666130812e-04 * solarMass
}
saturnPos :: Cl3
saturnPos = V3 (8.34336671824457987e+00) (4.12479856412430479e+00) (-4.03523417114321381e-01)
saturnVel :: Cl3
saturnVel = V3 (-2.76742510726862411e-03*dp) (4.99852801234917238e-03*dp) (2.30417297573763929e-05*dp)
saturnMass :: Cl3
saturnMass = R (2.85885980666130812e-04 * solarMass)
uranus :: Planet
uranus = Planet
{x = 1.28943695621391310e+01,y = -1.51111514016986312e+01,z = -2.23307578892655734e-01,
vx = 2.96460137564761618e-03*dp,vy = 2.37847173959480950e-03*dp, vz = -2.96589568540237556e-05*dp,
mass = 4.36624404335156298e-05 * solarMass
}
uranusPos :: Cl3
uranusPos = V3 (1.28943695621391310e+01) (-1.51111514016986312e+01) (-2.23307578892655734e-01)
uranusVel :: Cl3
uranusVel = V3 (2.96460137564761618e-03*dp) (2.37847173959480950e-03*dp) (-2.96589568540237556e-05*dp)
uranusMass :: Cl3
uranusMass = R (4.36624404335156298e-05 * solarMass)
neptune :: Planet
neptune = Planet
{x = 1.53796971148509165e+01,y = -2.59193146099879641e+01,z = 1.79258772950371181e-01,
vx = 2.68067772490389322e-03*dp,vy = 1.62824170038242295e-03*dp, vz = -9.51592254519715870e-05*dp,
mass = 5.15138902046611451e-05 * solarMass
}
neptunePos :: Cl3
neptunePos = V3 (1.53796971148509165e+01) (-2.59193146099879641e+01) (1.79258772950371181e-01)
neptuneVel :: Cl3
neptuneVel = V3 (2.68067772490389322e-03*dp) (1.62824170038242295e-03*dp) (-9.51592254519715870e-05*dp)
neptuneMass :: Cl3
neptuneMass = R (5.15138902046611451e-05 * solarMass)
daysPerYear :: Double
daysPerYear = 365.24
solarMass :: Double
solarMass = 4 * pi^(2 :: Int)
dp :: Double
dp = daysPerYear
dt :: Double
dt = 0.01
-- | 'Storable' instance of Planet
instance Storable Planet where
sizeOf _ = 8 * dblSz
alignment _ = dblSz
peekElemOff p i = peek (plusPtr p (i * sizeOf (undefined :: Planet)))
pokeElemOff p i e = poke (plusPtr p (i * sizeOf e)) e
peek :: Ptr Planet -> IO Planet
peek p = do
x' <- peek (offset 0)
y' <- peek (offset 1)
z' <- peek (offset 2)
vx' <- peek (offset 3)
vy' <- peek (offset 4)
vz' <- peek (offset 5)
mass' <- peek (offset 6)
return Planet {x=x',y=y',z=z',vx=vx',vy=vy',vz=vz',mass=mass'}
where
offset i = plusPtr (castPtr p :: Ptr Double) (i*8)
poke :: Ptr Planet -> Planet -> IO ()
poke p e = do
poke (offset 0) $ x e
poke (offset 1) $ y e
poke (offset 2) $ z e
poke (offset 3) $ vx e
poke (offset 4) $ vy e
poke (offset 5) $ vz e
poke (offset 6) $ mass e
where
offset i = plusPtr (castPtr p :: Ptr Double) (i*8)
dblSz :: Int
dblSz = sizeOf (undefined :: Double)
-- | 'Storable' instance of PlanetCl3
instance Storable PlanetCl3 where
sizeOf _ = 4 * cl3Sz
alignment _ = cl3Sz
peekElemOff p i = peek ( p `plusPtr` (i * sizeOf (undefined :: PlanetCl3)))
pokeElemOff p i e = poke (p `plusPtr` (i * sizeOf e)) e
peek :: Ptr PlanetCl3 -> IO PlanetCl3
peek p = do
pos' <- peek (offset 0)
vel' <- peek (offset 1)
mass' <- peek (offset 2)
return $ PlanetCl3 pos' vel' mass'
where
offset i = (castPtr p :: Ptr Cl3) `plusPtr` (i * cl3Sz)
poke :: Ptr PlanetCl3 -> PlanetCl3 -> IO ()
poke p e = do
poke (offset 0) $! posCl3 e
poke (offset 1) $! velCl3 e
poke (offset 2) $! massCl3 e
where
offset i = (castPtr p :: Ptr Cl3) `plusPtr` (i * cl3Sz)
cl3Sz :: Int
cl3Sz = sizeOf (undefined :: Cl3)
-- | 'Storable' instance of PlanetCl3_V3
instance Storable PlanetCl3_V3 where
sizeOf _ = sizeOf (undefined :: Cl3_V3) + sizeOf (undefined :: Cl3_V3) + sizeOf (undefined :: Cl3_R)
alignment _ = undefined -- perhaps this should be some integer power of 2, sizeOf PlanetCl3_V3 is 54bytes, should I pick 64bytes? Does Haskell even support a mov64?, or should it be 8 because everything at the core is a Double?
peekElemOff ptr idx = peek (castPtr ptr `plusPtr` (idx * sizeOf (undefined :: PlanetCl3_V3)))
pokeElemOff ptr idx e = poke (castPtr ptr `plusPtr` (idx * sizeOf (undefined :: PlanetCl3_V3))) e
peek :: Ptr PlanetCl3_V3 -> IO PlanetCl3_V3
peek ptr = do
pos :: Cl3_V3 <- peek (castPtr ptr `plusPtr` 0)
vel :: Cl3_V3 <- peek (castPtr ptr `plusPtr` (sizeOf pos))
m :: Cl3_R <- peek (castPtr ptr `plusPtr` (sizeOf pos) `plusPtr` (sizeOf vel))
return (PlanetCl3_V3 pos vel m)
poke :: Ptr PlanetCl3_V3 -> PlanetCl3_V3 -> IO ()
poke ptr (PlanetCl3_V3 pos vel m) = do
poke (castPtr ptr `plusPtr` 0) pos
poke (castPtr ptr `plusPtr` (sizeOf pos)) vel
poke (castPtr ptr `plusPtr` (sizeOf pos) `plusPtr` (sizeOf vel)) m
-- | 'pokeC' update the position
pokeC :: Ptr Planet -> Int -> Planet -> IO ()
pokeC p i e = do
poke (offset 0) $ x e
poke (offset 1) $ y e
poke (offset 2) $ z e
where
offset o = (castPtr p :: Ptr Double) `plusPtr` (o*8+i*64)
pokeCCl3 :: Ptr PlanetCl3 -> Int -> PlanetCl3 -> IO ()
pokeCCl3 p !i !e =
poke ((castPtr p :: Ptr Cl3) `plusPtr` (i * sizeOf (undefined :: PlanetCl3))) $! posCl3 e
-- Update the position of a planet at an index in the array of planets
pokeCCl3_V3 :: Ptr PlanetCl3_V3 -> Int -> PlanetCl3_V3 -> IO ()
pokeCCl3_V3 p idx e =
poke ((castPtr p) `plusPtr` (idx * sizeOf (undefined :: PlanetCl3_V3))) (posCl3_V3 e)
-- | 'pokeV' update the velocity
pokeV :: Ptr Planet -> Int -> Planet -> IO ()
pokeV p i e = do
poke (offset 3) $ vx e
poke (offset 4) $ vy e
poke (offset 5) $ vz e
where
offset o = (castPtr p::Ptr Double) `plusPtr` (o*8+i*64)
pokeVCl3 :: Ptr PlanetCl3 -> Int -> PlanetCl3 -> IO ()
pokeVCl3 p !i !e =
poke ((castPtr p :: Ptr Cl3) `plusPtr` (cl3Sz + i * sizeOf (undefined :: PlanetCl3))) $! velCl3 e
-- Update the velocity of a planet at an index in the array of planets
pokeVCl3_V3 :: Ptr PlanetCl3_V3 -> Int -> PlanetCl3_V3 -> IO ()
pokeVCl3_V3 p idx e =
poke ((castPtr p) `plusPtr` (idx * sizeOf (undefined :: PlanetCl3_V3)) `plusPtr` (sizeOf (undefined :: Cl3_V3))) (velCl3_V3 e)
-- | 'fromList' initialize the planets in the stack
fromList :: [Planet] -> IO (Ptr Planet)
fromList l = do
let len = length l
pa <- mallocBytes (len * sizeOf (undefined :: Planet))
let loop [] _ = return ()
loop (q:qs) i = do
poke (pa `plusPtr` (i * sizeOf (undefined :: Planet))) q
loop qs (i+1)
loop l 0
return pa
fromListCl3 :: [PlanetCl3] -> IO (Ptr PlanetCl3)
fromListCl3 l = do
let len = length l
pa <- mallocBytes (len * sizeOf (undefined :: PlanetCl3))
let loop [] _ = return ()
loop (q:qs) i = do
poke (pa `plusPtr` (i * sizeOf (undefined :: PlanetCl3))) q
loop qs (i+1)
loop l 0
return pa
fromListCl3_V3 :: [PlanetCl3_V3] -> IO (Ptr PlanetCl3_V3)
fromListCl3_V3 l = do
let len = length l
pa <- mallocBytes (len * sizeOf (undefined :: PlanetCl3_V3))
let loop [] _ = return ()
loop (q:qs) idx = do
poke (pa `plusPtr` (idx * sizeOf (undefined :: PlanetCl3_V3))) q
loop qs (idx+1)
loop l 0
return pa
-------------------------------------------------------------------
--
-- The Computer Language
-- Benchmarks Game
-- Revised BSD license
--
-- This is a specific instance of the Open Source Initiative (OSI) BSD license template.
--
-- Copyright © 2004-2008 Brent Fulgham, 2005-2017 Isaac Gouy
--
-- All rights reserved.
--
-- 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 Computer Language Benchmarks Game" nor the name of "The Computer Language Shootout Benchmarks" 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 AND 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 COPYRIGHT OWNER OR 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.
--
--
--
-- The Computer Language Benchmarks Game
-- http://benchmarksgame.alioth.debian.org/
--
-- Contributed by Branimir Maksimovic
--
-------------------------------------------------------------------