haskell:

Donald Bruce Stewart wrote:

...

haskell:

...

There is a new combined benchmark, "partial sums" that subsumes several earlier benchmarks and runs 9 different numerical calculations:

http://haskell.org/hawiki/PartialSumsEntry

Ah! I had an entry too. I've posted it on the wiki. I was careful to watch that all loops are compiled into nice unboxed ones in the Core. It seems to run a little bit faster than your more abstracted code.

Timings on the page.

Also, -fasm seems to only be a benefit on the Mac, as you've pointed out previously. Maybe you could check the times on the Mac too?

-- Don

Yeah. I had not tried all the compiler options. Using -fasm is slower on this for me as well. I suspect that since your code will beat the entries that have been posted so far, so I thin you should submit it.

ok, I'll submit it.

Also, could you explain how to check the Core (un)boxing in a note on the (new?) wiki? I would be interested in learning that trick.

Ah, i just do: ghc A.hs -O2 -ddump-simpl | less and then read the Core, keeping an eye on the functions I'm interested in, and checking they're compiling to the kind of loops I'd write by hand. This is particularly useful for the kinds of tight numeric loops used in some of the shootout entries. Cheers, Don

Joel Koerwer wrote:

On 1/26/06, *Donald Bruce Stewart* mailto:dons@cse.unsw.edu.au> wrote:

Ah, i just do: ghc A.hs -O2 -ddump-simpl | less and then read the Core, keeping an eye on the functions I'm interested in, and checking they're compiling to the kind of loops I'd write by hand. This is particularly useful for the kinds of tight numeric loops used in some of the shootout entries.

Cheers, Don

In that case could you describe the kind of loops you'd write by hand?

See below for the pseudo-code loop and the Haskell version.

Seriously. And perhaps typical problems/fixes when the compiler doesn't produce what you want.

We don't have any fixes.

Thanks, Joel

More discussion and code is at http://haskell.org/hawiki/NbodyEntry The compiler produces code that runs 4 times slower than OCaml in our current best attempt at programming against a 40 element (IOUArray Int Double). The final programs speed is very architecture dependent, but more frustrating is that small referentially transparent changes to the source code produce up to factor-of-two fluctuations in run time. The small numeric functions in the shootout, where there is a recursive function with 1 or 2 parameters (Double's), perform quite well. But manipulating this medium number of Double's to model the solar system has been too slow. The main loop for the 5 planets looks quite simple in pseudo-c: deltaTime = 0.01 for (i=0 ; i<5; ++i) { "get mass m, position (x,y,z), velocity (vx,vy,vz) of particle number i" for (j=(i+1); j<5; ++j) { "get mass, position, velocity of particle j" dxyx = "position of i" - "position of j" mag = deltaTime /(length of dxyz)^3 "velocity of j" += "mass of i" * mag * dxyz "velocity of i" -= "mass of j" * mag * dxyz } "position of i" += deltaTime * "velocity of i" } Note that the inner loop "for j" starts a "j=(i+1)". The best performing Haskell code, for this loop, so far is: -- Offsets for each field x = 0; y = 1; z = 2; vx= 3; vy= 4; vz= 5; m = 6 -- This is the main code. Essentially all the time is spent here advance n = when (n > 0) $ updateVel 0 >> advance (pred n) where updateVel i = when (i <= nbodies) $ do let i' = (.|. shift i 3) im <- unsafeRead b (i' m) ix <- unsafeRead b (i' x) iy <- unsafeRead b (i' y) iz <- unsafeRead b (i' z) ivx <- unsafeRead b (i' vx) ivy <- unsafeRead b (i' vy) ivz <- unsafeRead b (i' vz) let updateVel' ivx ivy ivz j = ivx `seq` ivy `seq` ivz `seq` if j > nbodies then do unsafeWrite b (i' vx) ivx unsafeWrite b (i' vy) ivy unsafeWrite b (i' vz) ivz else do let j' = (.|. shiftL j 3) jm <- unsafeRead b (j' m) dx <- liftM (ix-) (unsafeRead b (j' x)) dy <- liftM (iy-) (unsafeRead b (j' y)) dz <- liftM (iz-) (unsafeRead b (j' z)) let distance = sqrt (dx*dx+dy*dy+dz*dz) mag = 0.01 / (distance * distance * distance) addScaled3 (3 .|. (shiftL j 3)) ( im*mag) dx dy dz let a = -jm*mag ivx' = ivx+a*dx ivy' = ivy+a*dy ivz' = ivz+a*dz updateVel' ivx' ivy' ivz' $! (j+1) updateVel' ivx ivy ivz $! (i+1) addScaled (shiftL i 3) 0.01 (3 .|. (shiftL i 3)) updateVel (i+1) -- Helper functions addScaled i a j | i `seq` a `seq` j `seq` False = undefined -- stricitfy addScaled i a j = do set i1 =<< liftM2 scale (unsafeRead b i1) (unsafeRead b j1) set i2 =<< liftM2 scale (unsafeRead b i2) (unsafeRead b j2) set i3 =<< liftM2 scale (unsafeRead b i3) (unsafeRead b j3) where scale old new = old + a * new i1 = i; i2 = succ i1; i3 = succ i2; j1 = j; j2 = succ j1; j3 = succ j2; addScaled3 i a jx jy jz | i `seq` a `seq` jx `seq` jy `seq` jz `seq` False = undefined addScaled3 i a jx jy jz = do set i1 =<< liftM (scale jx) (unsafeRead b i1) set i2 =<< liftM (scale jy) (unsafeRead b i2) set i3 =<< liftM (scale jz) (unsafeRead b i3) where scale new old = a * new + old i1 = i; i2 = succ i1; i3 = succ i2;

joelkoerwer:

Thanks Chris. I was actually asking about analyzing Core output in general. I'm well aware of the problems we're having with the nbody entry. I'm convinced my list based version can go faster than it is now. That's why I was asking if Don could put together a few notes on how to optimize inner loops using -ddump-simpl and the resulting Core code.

Here's a brief introduction. I intend to write up (on the performance page on the wiki) a list of things we've done to improve the shootout entries. N.B we're now the 3rd *fastest* language, behind C and only a little behind D (a C varient) !! Consider the partial sums problem: wiki: http://www.haskell.org/hawiki/PartialSumsEntry site: http://shootout.alioth.debian.org/gp4/benchmark.php?test=partialsums&lang=ghc&id=2 What follows is a discussion of the steps I took to improve the performance of this code. Here's the naive translation of the Clean entry (which was fairly quick): Lots of math in a tight loop.

import System; import Numeric

main = do n <- getArgs >>= readIO . head let sums = loop 1 n 1 0 0 0 0 0 0 0 0 0 fn (s,t) = putStrLn $ (showFFloat (Just 9) s []) ++ "\t" ++ t mapM_ (fn :: (Double, String) -> IO ()) (zip sums names)

names = ["(2/3)^k", "k^-0.5", "1/k(k+1)", "Flint Hills", "Cookson Hills" , "Harmonic", "Riemann Zeta", "Alternating Harmonic", "Gregory"]

loop k n alt a1 a2 a3 a4 a5 a6 a7 a8 a9 | k > n = [ a1, a2, a3, a4, a5, a6, a7, a8, a9 ] | otherwise = loop (k+1) n (-alt) (a1 + (2/3) ** (k-1)) (a2 + k ** (-0.5)) (a3 + 1 / (k * (k + 1))) (a4 + 1 / (k*k*k * sin k * sin k)) (a5 + 1 / (k*k*k * cos k * cos k)) (a6 + 1 / k) (a7 + 1 / (k*k)) (a8 + alt / k) (a9 + alt / (2 * k - 1))

Compiled with "-O2". However, the performance is _really_ bad :/ Somewhere greater than 128M heap, in fact eventually running out of memory on my laptop. A classic space leak. (2) So look at the generated core. "ghc -o naive Naive.hs -O2 -ddump-simpl | less" And we find that our loop has the following type:

$sloop_r2U6 :: GHC.Prim.Double# -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Prim.Double# -> [GHC.Float.Double]

Hmm. Ok, I certainly don't want boxed doubles in such a tight loop. (3) My next step is to encourage GHC to unbox this loop, by providing some strictness annotations. Now the loop looks like this:

loop k n alt a1 a2 a3 a4 a5 a6 a7 a8 a9 | () !k !n !False = undefined | k > n = [ a1, a2, a3, a4, a5, a6, a7, a8, a9 ] | otherwise = loop (k+1) n (-alt) (a1 + (2/3) ** (k-1)) (a2 + k ** (-0.5)) (a3 + 1 / (k * (k + 1))) (a4 + 1 / (k*k*k * sin k * sin k)) (a5 + 1 / (k*k*k * cos k * cos k)) (a6 + 1 / k) (a7 + 1 / (k*k)) (a8 + alt / k) (a9 + alt / (2 * k - 1)) where x ! y = x `seq` y

I've played a little game here, using ! for `seq`, reminiscent of the new !-pattern proposal for strictness. Let's see how this compiles. Here's the Core:

$sloop_r2Vh :: GHC.Prim.Double# -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Float.Double -> GHC.Prim.Double# -> [GHC.Float.Double]

Ok, so it unboxed one extra argument. Let's see if we can get them all unboxed. Strictify all args, and GHC produces an inner loop of:

$sloop_r2WS :: GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> GHC.Prim.Double# -> [GHC.Float.Double]

Ah! perfect. Let's see how that runs: $ ghc Naive.hs -O2 -no-recomp $ time ./a.out 2500000 3.000000000 (2/3)^k 3160.817621887 k^-0.5 0.999999600 1/k(k+1) 30.314541510 Flint Hills 42.995233998 Cookson Hills 15.309017155 Harmonic 1.644933667 Riemann Zeta 0.693146981 Alternating Harmonic 0.785398063 Gregory ./a.out 2500000 4.45s user 0.02s system 99% cpu 4.482 total (4) Not too bad. No space leak and quite zippy. But let's see what more can be done. First, double arithmetic usually benefits from -fexcess-precision, and cranking up the flags to gcc: paprika$ ghc Naive.hs -O2 -fexcess-precision -optc-O3 -optc-ffast-math -no-recomp paprika$ time ./a.out 2500000 3.000000000 (2/3)^k 3160.817621887 k^-0.5 0.999999600 1/k(k+1) 30.314541510 Flint Hills 42.995233998 Cookson Hills 15.309017155 Harmonic 1.644933667 Riemann Zeta 0.693146981 Alternating Harmonic 0.785398063 Gregory ./a.out 2500000 3.71s user 0.01s system 99% cpu 3.726 total Even better. Now, let's dive into the Core to see if there are any optimisation opportunites that GHC missed. So add -ddump-simpl and peruse the output. (5) Looking at the Core, I see firstly that some of the common subexpressions haven't been factored out: case [GHC.Float.Double] GHC.Prim./## 1.0 (GHC.Prim.*## (GHC.Prim.*## (GHC.Prim.*## (GHC.Prim.*## sc10_s2VS sc10_s2VS) sc10_s2VS) (GHC.Prim.sinDouble# sc10_s2VS)) (GHC.Prim.sinDouble# sc10_s2VS)) Multiple calls to sin. Hmm :/ And similar for cos, as well as k*k. Not sure why GHC isn't removing these (SimonM?), so let's do that by hand, and the inner loop looks like:

loop k n alt a1 a2 a3 a4 a5 a6 a7 a8 a9 | () !k !n !alt !a1 !a2 !a3 !a4 !a5 !a6 !a7 !a8 !a9 !False = undefined | k > n = [ a1, a2, a3, a4, a5, a6, a7, a8, a9 ] | otherwise = loop (k+1) n (-alt) (a1 + (2/3) ** (k-1)) (a2 + k ** (-0.5)) (a3 + 1 / (k * (k + 1))) (a4 + 1 / (k3 * sk * sk)) (a5 + 1 / (k3 * ck * ck)) (a6 + 1 / k) (a7 + 1 / k2) (a8 + alt / k) (a9 + alt / (2 * k - 1)) where sk = sin k ; ck = cos k; k2 = k * k; k3 = k2 * k; x ! y = x `seq` y

looking at the Core shows the sins are now allocated and shared:

let a9_s2MI :: GHC.Prim.Double# a9_s2MI = GHC.Prim.sinDouble# sc10_s2Xa

So the common expressions are floated out, and it now runs: paprika$ time ./a.out 2500000 3.000000000 (2/3)^k 3160.817621887 k^-0.5 0.999999600 1/k(k+1) 30.314541510 Flint Hills 42.995233998 Cookson Hills 15.309017155 Harmonic 1.644933667 Riemann Zeta 0.693146981 Alternating Harmonic 0.785398063 Gregory ./a.out 2500000 3.29s user 0.00s system 99% cpu 3.290 total Faster. So we gained 12% by floating out those common expressions. (6) Finally, another trick. When I checked the C entry, it used an integer for the k parameter to the loop, and cast it to a double for the math each time around, so perhaps we can make it an Int parameter. And secondly, the alt parameter only has it's sign flipped each time, so perhaps we can factor out the alt / k arg (it's either 1 / k or -1 on k), saving a division. So the final loop looks like:

loop i n alt a1 a2 a3 a4 a5 a6 a7 a8 a9 | () !i !n !alt !a1 !a2 !a3 !a4 !a5 !a6 !a7 !a8 !a9 !False = undefined | k > n = [ a1, a2, a3, a4, a5, a6, a7, a8, a9 ] | otherwise = loop (i+1) n (-alt) (a1 + (2/3) ** (k-1)) (a2 + k ** (-0.5)) (a3 + 1 / (k * (k + 1))) (a4 + 1 / (k3 * sk * sk)) (a5 + 1 / (k3 * ck * ck)) (a6 + dk) (a7 + 1 / k2) (a8 + alt * dk) (a9 + alt / (2 * k - 1)) where k = fromIntegral i dk = 1/k sk = sin k ck = cos k k2 = k * k k3 = k2 * k x ! y = x `seq` y

Checking the generated C code shows that the same operations are generated as the C entry. And it runs: $ time ./a.out 2500000 3.000000000 (2/3)^k 3160.817621887 k^-0.5 0.999999600 1/k(k+1) 30.314541510 Flint Hills 42.995233998 Cookson Hills 15.309017155 Harmonic 1.644933667 Riemann Zeta 0.693146981 Alternating Harmonic 0.785398063 Gregory ./a.out 2500000 3.17s user 0.01s system 99% cpu 3.206 total This entry in fact runs faster than the original (though not the new vectorised entry) optimised C entry (and faster than all other languages): http://shootout.alioth.debian.org/gp4/benchmark.php?test=partialsums&lang=all So, by carefully tweaking things, we first squished a space leak, and then gained another 30%. In summary: * Check the Core that is generated * Watch out for optimisations that are missed * Read the generated C for the tight loops. * Make sure tight loops are unboxed * Use -fexcess-precision and -optc-ffast-math for doubles This is roughly the process I used for the other shootout entries. Cheers, Don

Hello Donald, Wednesday, February 01, 2006, 8:00:04 AM, you wrote: DBS> Here's a brief introduction. I intend to write up (on the performance page on DBS> the wiki) a list of things we've done to improve the shootout entries. N.B DBS> we're now the 3rd *fastest* language, behind C and only a little behind D (a C DBS> varient) !! 3rd fastest or 3rd overall, counting program lines and so on? -- Best regards, Bulat mailto:bulatz@HotPOP.com

--- Chris Kuklewicz wrote: -snip,snip-

It is 3rd fastest. Looking at Just Memory Use, Haskell is 8th Looking at Just Lines Of Code, Haskell is 1st Lookat at the 1:1:1 even balance Haskell is 1st

Programmer skill and effort really does matter ;-) Congratulations. __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com

--- Ketil Malde wrote:

Isaac Gouy writes:

...

Programmer skill and effort really does matter ;-)

Yes, more so, than any inherent language disadvantage, perhaps, which happens to be the general lesson from the ICFP contests as well. Any idea if other languages have seen similar efforts?

FreePascal and Smart Effiel, somewhat - and there have been excellent individual efforts with Lua and Tcl and ... imo the Haskell Cafe discussions and wiki have been a more open and shared learning experience than we usually see, and maybe some of the success stems from that collaboration and competition. __________________________________________________ Do You Yahoo!? Tired of spam? Yahoo! Mail has the best spam protection around http://mail.yahoo.com

Joel Koerwer wrote:

Don, that's a great little mini tutorial, exactly what I was hoping for. I'm looking forward to learning more tricks.

On an unrelated note, I have an STUArray nbody. I haven't really looked closely at the chris+dons version, but I suspect they amount to doing the same thing. I get commensurate runtimes at least. But I'll post it on the wiki in a while in case there is some optimization I missed.

On an even more unrelated note, I get slower runtimes with -optc-O3 and -optc-ffast-math than without. Individually or in tandem. Odd. This is ghc 6.4.1 and gcc 4.0.3 on a (Banias) Pentium M.

More architecture benchmarking: On a powerbook G4, Joel Koerwer's entry on http://haskell.org/hawiki/NbodyEntry runs faster than dons+chris. And I edited it to make it smaller, but by hoisting 'size' and 'dt' I also made it run quite a bit faster. It now takes 1.7x less time than dons+chris. Don, could you check the speed on your architecture? -- Chris

Hey this is great. Chris your improvements are awesome. I mean the speed is nice, but you really cleaned up the code. There's an extraneous call to energy in the second runST block, but it should be insignificant. Also, -fglasgow-exts is necessary for the left-hand-side type declarations of size and dt. One question, in: calcMomentum (i+1) $! (px+vx*m,py+vy*m,pz+vz*m) that $! doesn't actually do much for a tuple, does it? Of course, there's not much point in further optimizing the initialization routine, as we'd never be able to detect the difference in runtime. The shootout has been a great learning tool for me :-) Thanks to Chris, Don, and the rest of the Haskell community. Joel