[Haskell-cafe] Squashing space leaks

[Greg Buchholz]

    I've been dabbling with implementing the "n-body" benchmark over at
the "Great Computer Language Shootout"...

http://shootout.alioth.debian.org/benchmark.php?test=nbody

...and I've stumbled on to a problem with space leaks.  The code below
works great for low numbers of iterations, but crashes and burns with
large numbers.  You can temporarily fix it up by specifying larger
stack and heap sizes (with the +RTS -K & -H options), but this only gets
you so far.  I haven't been able to glean much information from the
profiling reports.  You can see the reports I generated (like the
biography report) over at <http://sleepingsquirrel.org/nbody>.  I've
tried several different variations on the "advance" function, as well as
trying "seq" in numerous spots, but I haven't stumbled on to the right
combination yet.  Any pointers?

Thanks,

Greg Buchholz

>-- n-body simulation
>-- compile:  ghc -O2 -o nbody nbody.hs
>-- run:  nbody 100000
>
>import System(getArgs)
>import Numeric
>
>data Vec = V !Double !Double !Double deriving Show
>instance Num Vec where
>    (V a b c) + (V x y z) = (V (a+x) (b+y) (c+z))
>    (V a b c) - (V x y z) = (V (a-x) (b-y) (c-z))
>    fromInteger 0 = V 0.0 0.0 0.0 -- for sum function
>instance Eq Vec where (V a b c) == (V x y z) = (a==x) && (b==y) && (c==z)
>
>dot (V a b c) (V x y z) = a*x + b*y + c*z
>scale (V a b c) n = V (n*a) (n*b) (n*c)
>mag (V x y z) =  sqrt (x*x + y*y + z*z)
>
>data Planet = Planet Vec Vec Double deriving Show --Position Velocity Mass
>dist (Planet p1 _ _) (Planet p2 _ _) = mag $ p1 - p2
>mass (Planet _ _ m) = m
>vel  (Planet _ v _) = v
>pos  (Planet p _ _) = p
>
>main = do 
>        [arg] <- getArgs
>        let iter = read arg
>        let n = 5::Int
>        let bodies = offset_momentum n [sun, jupiter, saturn, neptune, uranus]
>        let begin = energy n bodies
>        putStrLn $ showFFloat (Just 9) begin ""
>        let final = (iterate (advance n 0.01) bodies)
>        let end = energy n (final !! iter)
>        putStrLn $ showFFloat (Just 9) end ""
>
>days_per_year = 365.24 
>solar_mass = 4.0 * pi * pi
>
>advance :: Int -> Double -> [Planet] -> [Planet]
>advance n dt (p:ps) = take n ((Planet (pos p + delta_x)
>                              (new_v) (mass p) ) : advance n dt (ps++[p]))
> where
>   delta_v = sum (map (\q -> 
>                  (pos p - pos q) `scale` ((mass q)*dt/(dist p q)^3)) ps)
>   new_v   = (vel p) - delta_v
>   delta_x = new_v `scale` dt 
>
>energy:: Int -> [Planet] -> Double
>energy n ps = kinetic - potential
>  where
>    kinetic   = 0.5 * (sum (map (\q->(mass q)*((vel q) `dot` (vel q))) ps))
>    potential = sum [(mass (ps!!i))*(mass (ps!!j))/(dist (ps!!i) (ps!!j))
>                      | i<-[0..n-1], j<-[i+1..n-1]] 
>
>offset_momentum n ((Planet p v m):ps) = (Planet p new_v m):ps
>  where new_v = (sum (map (\n->(vel n) `scale` (mass n)) ps)) 
>                `scale` ((-1.0)/solar_mass)
>
>{-
>rotations 0 _ = []
>rotations n xss@(x:xs) = flipped:rotations (n-1) flipped
>    where flipped = xs++[x]
>
>advance :: Int -> Double -> [Planet] -> [Planet]
>advance n dt pss = map (adv dt) (rotations n pss)
>
>adv dt (p:ps) = Planet (pos p + delta_x) (new_v) (mass p) 
> where
>   delta_v = sum (map (\q -> 
>                  (pos p - pos q) `scale` ((mass q)*dt/(dist p q)^3)) ps)
>   new_v   = (vel p) - delta_v
>   delta_x = new_v `scale` dt 
>--}
>jupiter = (Planet 
> (V 4.84143144246472090e+00 (-1.16032004402742839e+00)
> (-1.03622044471123109e-01))
> (V ( 1.66007664274403694e-03 * days_per_year)
>    ( 7.69901118419740425e-03 * days_per_year)
>    ((-6.90460016972063023e-05) * days_per_year)) 
> (9.54791938424326609e-04 * solar_mass))
>
>saturn = (Planet
> (V 8.34336671824457987e+00 4.12479856412430479e+00
> (-4.03523417114321381e-01))
> (V (-2.76742510726862411e-03 * days_per_year) 
>    (4.99852801234917238e-03 * days_per_year)
>    (2.30417297573763929e-05 * days_per_year))
> (2.85885980666130812e-04 * solar_mass))
>
>uranus = (Planet
> (V 1.28943695621391310e+01 (-1.51111514016986312e+01)
> (-2.23307578892655734e-01))
> (V (2.96460137564761618e-03 * days_per_year)
>    (2.37847173959480950e-03 * days_per_year)
>    (-2.96589568540237556e-05 * days_per_year))
> (4.36624404335156298e-05 * solar_mass))
>
>neptune = (Planet
> (V 1.53796971148509165e+01 (-2.59193146099879641e+01)
> 1.79258772950371181e-01)
> (V (2.68067772490389322e-03 * days_per_year)
>    (1.62824170038242295e-03 * days_per_year)
>    (-9.51592254519715870e-05 * days_per_year))
> (5.15138902046611451e-05 * solar_mass))
>
>sun = (Planet (V 0.0 0.0 0.0) (V 0.0 0.0 0.0) solar_mass)