I am new to functional programming, and now learn Haskell. As an exercise I decided to implement the explicit Euler method for 1D linear diffusion equation. While the code below works correctly, I am not happy about its performance. In fact, I am concerned with memory consumption. I believe that it is related to lazy evaluation, but cannot figure out how I can reduce its memory usage.
The idea of the algorithm is really simple, to make it clear in imperative terms: it takes an `array', and to every inner point it adds a value, which is calculated as a combination of the values in the point itself and in its neighbors. Boundary points are special cases.
So, this is my Euler1D.hs module:
module Euler1D
( stepEuler
, makeu0
) where
-- impose zero flux condition
zeroflux :: (Floating a) => a -> [a] -> [a]
zeroflux mu (boundary:inner:xs) = [boundary+mu*2*(inner-boundary)]
-- one step of integration
stepEuler :: (Floating a) => a -> [a] -> [a]
stepEuler mu u@(x:xs) = (applyBC . (diffused mu)) u
where
diffused mu (left:x:[]) = [] -- ignore outer points
diffused mu (left:x:right:xs) = -- integrate inner points
(x+mu*(left+right-2*x)) : diffused mu (x:right:xs)
applyBC inner = (lbc u') ++ inner ++ (rbc u') -- boundary conditions
where u' = [head u] ++ inner ++ [last u]
lbc = zeroflux mu -- left boundary
rbc = (zeroflux mu) . reverse -- right boundary
-- initial condition
makeu0 :: Int -> [Double]
makeu0 n = [ ((^2) . sin . (pi*) . xi) x | x <- [0..n]]
where xi x = fromIntegral x / fromIntegral n
And a simple Main.hs:
module Main where
import System ( getArgs )
import Euler1D
main = do
args <- getArgs
let n = read $ head args :: Int
let u0 = makeu0 n
let un = stepEuler 0.5 u0
putStrLn $ show $ sum un
For comparison, I also wrote a pure C implementation.
Now, if I try to run Haskell implementation for a sufficiently large size of the array n
, I have:
$ time ./eulerhs 200000
100000.00000000112
real 0m3.552s
user 0m3.304s
sys 0m0.128s
For comparison, C version is faster by almost two orders of magnitude:
$ time ./eulerc 200000
100000
real 0m0.088s
user 0m0.048s
sys 0m0.008s
EDIT: This comparison is not really fair, because Haskell version is compiled with profiling flags, and C is not. If I compile both programs with
-O2
and both without profiling flags, I can increasen
. In this casetime ./eulerhs 1000000
takes 0m2.236s, whiletime ./eulerc 1000000
takes only 0m0.293s. So the problem still remains with all optimizations and without profiling, it is only offset.I would like also to note, that memory allocation of the Haskell program seems to grow lineary with
n
. This is probably OK.
But the worst are memory requirements. My Haskell version requires over 100MB (my estimation of the bare minimum in C is 4MB). I think this may be the source of the problem. According to profiling report the program spends 85% of time in GC, and
total time = 0.36 secs (18 ticks @ 20 ms)
total alloc = 116,835,180 bytes (excludes profiling overheads)
COST CENTRE MODULE %time %alloc
makeu0 Euler1D 61.1 34.9
stepEuler Euler1D 33.3 59.6
CAF:sum Main 5.6 5.5
I was surprized to see that makeu0
is so expensive. I decided that this is due to its lazy evaluation (if its thunks remain in the memory until the end of stepEuler
).
I tried this change in Main.hs
:
let un = u0 `seq` stepEuler 0.5 u0
but did not notice any difference. I have no idea how to reduce memory usage in stepEuler
. So, my questions are:
- Is there a way in Haskell to build lists / do list comprehensions strictly? In this case there is no benefit to keep it lazy.
- How can I reduce overall memory usage in this case? I suppose, I have to make something strict, but fail to see what. In other words, if I have to put some
seq
s and bangs, where and why? - And finally, the most important, what is the best strategy to identify such costly constructs?
I did read a chapter on profiling and optimization in Real World Haskell, but it remains unclear how exactly I can decide what should be strict and what not.
Please forgive me such a long post.
EDIT2: As suggested by A. Rex in comments, I tried running both programs in valgrind. And this is what I observed. For Haskell program (
n
=200000) it found:malloc/free: 33 allocs, 30 frees, 84,109 bytes allocated. ... checked 55,712,980 bytes.
And for C program (after a small fix):
malloc/free: 2 allocs, 2 frees, 3,200,000 bytes allocated.
So, it appears that while Haskell allocates much smaller memory blocks, it does it often, and due to delay in garbage collection, they accumulate and remain in memory. So, I have another question:
- Is it possible to avoid a lot of small allocations in Haskell? Basically, to declare, that I need to process the whole data structure rather than only its fragments on demand.