所以,在那种情况下,我开始玩测试用例.我正在研究的其中一个是对经典的4阶Runge-Kutta方法的天真,直接的实现,适用于相当微不足道的IVP dy / dt = -y; y(0)= 1,得到y = e ^ -t.我在Haskell和C中编写了一个完全直接的实现(稍后我将发布). Haskell版本非常简洁,当我看到它时,内部给了我温暖的模糊,但C版本(实际上解析起来并不可怕)的速度超过了两倍.
我意识到比较两种不同语言的表现并非百分之百;直到我们都死C的那一天,很可能总是拥有作为表演之王的王冠,特别是手动优化的C代码.我不是试图让我的Haskell实现运行与我的C实现一样快.但我非常肯定,如果我更清楚自己在做什么,那么我可以从这个特定的Haskell实现中获得更高的速度.
Haskell版本在OS X 10.8.4的GHC 7.6.3下用-02编译,C版本用Clang编译,我没有给它标记. Haskell版本在跟踪时间时平均约为0.016秒,而C版本约为0.006秒.
这些时间考虑了二进制文件的整个运行时间,包括输出到stdout,这显然占了一些开销,但我确实通过使用-prof -auto-all重新编译并运行了一些来对GHC二进制文件进行一些分析. RTS -p还用RTS -s查看GC统计数据.我并没有真正理解我所看到的所有内容,但似乎我的GC并没有失控,尽管可能会受到一点点控制(5%,生产率约为93%用户,~85已经过去的百分比总和)并且大部分生产时间都花在了函数iterateRK上,我知道在写这个函数时会很慢但是对于我来说如何清理它并不是很明显.我意识到我在使用List时可能会受到惩罚,无论是在持续存在还是将结果倾倒到stdout中的懒惰.
我究竟做错了什么?什么库函数或Monadic魔法我悲惨地不知道我可以用来清理iterateRK?什么是学习如何成为GHC剖析摇滚明星的好资源?
RK.hs
rk4 :: (Double -> Double -> Double) -> Double -> Double -> Double -> Double
rk4 y' h t y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
where k1 = y' t y
k2 = y' (t + h/2) (y + ((h/2) * k1))
k3 = y' (t + h/2) (y + ((h/2) * k2))
k4 = y' (t + h) (y + (h * k3))
iterateRK y' h t0 y0 = y0:(iterateRK y' h t1 y1)
where t1 = t0 + h
y1 = rk4 y' h t0 y0
main = do
let y' t y = -y
let h = 1e-3
let y0 = 1.0
let t0 = 0
let results = iterateRK y' h t0 y0
(putStrLn . show) (take 1000 results)
RK.c
#include<stdio.h>
#define ITERATIONS 1000
double rk4(double f(double t, double x), double h, double tn, double yn)
{
double k1, k2, k3, k4;
k1 = f(tn, yn);
k2 = f((tn + h/2), yn + (h/2 * k1));
k3 = f((tn + h/2), yn + (h/2 * k2));
k4 = f(tn + h, yn + h * k3);
return yn + (h/6) * (k1 + 2*k2 + 2*k3 + k4);
}
double expDot(double t, double x)
{
return -x;
}
int main()
{
double t0, y0, tn, yn, h, results[ITERATIONS];
int i;
h = 1e-3;
y0 = 1.0;
t0 = 0.0;
yn = y0;
for(i = 0; i < ITERATIONS; i++)
{
results[i] = yn;
yn = rk4(expDot, h, tn, yn);
tn += h;
}
for(i = 0; i < ITERATIONS; i++)
{
printf("%.10lf", results[i]);
if(i != ITERATIONS - 1)
printf(", ");
else
printf("\n");
}
return 0;
}
使用大小增加的程序会导致堆栈溢出:
Stack space overflow: current size 8388608 bytes. Use `+RTS -Ksize -RTS' to increase it.
这可能是由于太多的懒惰造成的.
查看按类型细分的堆配置文件,您将获得以下内容:
(注意:我修改了你的程序,左撇子指出)
这看起来不太好.您不应该为算法要求线性空间.您似乎持有的Double值超过了要求.严格解决问题:
{-# LANGUAGE BangPatterns #-}
iterateRK :: (Double -> Double -> Double) -> Double -> Double -> Double -> [Double]
iterateRK y' !h !t0 !y0 = y0:(iterateRK y' h t1 y1)
where t1 = t0 + h
y1 = rk4 y' h t0 y0
通过此修改,新堆配置文件如下所示:
这看起来好多了,内存使用率要低得多. -sstderr`还确认修改后我们只花费垃圾收集器总时间的2.5%:
%GC time 2.5% (2.9% elapsed)
现在,haskell版本仍然比C版本慢40%(使用用户时间):
$time ./tesths; time ./testc 2.47e-321 ./tesths 0,73s user 0,01s system 86% cpu 0,853 total 2.470328e-321 ./testc 0,51s user 0,01s system 95% cpu 0,549 total
增加迭代次数并使用堆分配的数组为C中的结果存储再次降低了差异:
time ./tesths; time ./testc 2.47e-321 ./tesths 18,25s user 0,04s system 96% cpu 19,025 total 2.470328e-321 ./testc 16,98s user 0,14s system 98% cpu 17,458 total
这仅相差约9%.
但我们仍然可以做得更好.使用stream-fusion软件包,我们可以完全消除列表,同时仍保持解耦.以下是包含该优化的完整代码:
{-# LANGUAGE BangPatterns #-}
import qualified Data.List.Stream as S
rk4 :: (Double -> Double -> Double) -> Double -> Double -> Double -> Double
rk4 y' !h !t !y = y + (h/6) * (k1 + 2*k2 + 2*k3 + k4)
where k1 = y' t y
k2 = y' (t + h/2) (y + ((h/2) * k1))
k3 = y' (t + h/2) (y + ((h/2) * k2))
k4 = y' (t + h) (y + (h * k3))
iterateRK :: (Double -> Double -> Double) -> Double -> Double -> Double -> [Double]
iterateRK y' h = curry $S.unfoldr $\(!t0, !y0) -> Just (y0, (t0 + h, rk4 y' h t0 y0))
main :: IO ()
main = do
let y' t y = -y
let h = 1e-3
let y0 = 1.0
let t0 = 0
let results = iterateRK y' h t0 y0
print $S.head $(S.drop (pred 10000000) results)
我跟着:
$ghc -O2 ./test.hs -o tesths -fllvm
以下是时间安排:
$time ./tesths; time ./testc 2.47e-321 ./tesths 15,85s user 0,02s system 97% cpu 16,200 total 2.470328e-321 ./testc 16,97s user 0,18s system 97% cpu 17,538 total
现在我们甚至比C快一点,因为我们没有分配.要对C程序进行类似的转换,我们必须将两个循环合并为一个并且松散好的抽象.即使这样,它也只有haskell一样快:
$time ./tesths; time ./testc 2.47e-321 ./tesths 15,86s user 0,01s system 98% cpu 16,141 total 2.470328e-321 ./testc 15,88s user 0,02s system 98% cpu 16,175 total
