{-# LANGUAGE ViewPatterns #-} {-# LANGUAGE BangPatterns #-} ------------------------------------------------------------------ -- | -- 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) -- 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 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)] ] 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 >>= 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 >>= 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" -- | '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 -- | '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} -- | '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 -- | '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 -- | '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 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] 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 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 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 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 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) -- | '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 -- | '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 -- | '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 ------------------------------------------------------------------- -- -- 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 -- -------------------------------------------------------------------