Haskell - 优化微分方程求解器

Question

我正在学习 Haskell，并尝试尽可能快地编写 C 代码。对于这个练习，我正在编写 Euler integrator对于一个简单的一维物理系统。

C 代码使用 GCC 4.5.4 和 -O3 编译。 .它在 1.166 秒内运行。

Haskell 代码使用 GHC 7.4.1 和 -O3 编译。 .它在 21.3 秒内运行。

如果我用 -O3 -fllvm 编译 Haskell ，它在 4.022 秒内运行。

那么，我是否缺少优化 Haskell 代码的内容？

PS.:我使用了以下参数:1e-8 5 .

C代码:

#include <stdio.h>
double p, v, a, t;

double func(double t) {
  return t * t;
}

void euler(double dt) {
  double nt = t + dt;
  double na = func(nt);
  double nv = v + na * dt;
  double np = p + nv * dt;

  p = np;
  v = nv;
  a = na;
  t = nt;
}

int main(int argc, char ** argv) {
  double dt, limit;
  sscanf(argv[1], "%lf", &dt);
  sscanf(argv[2], "%lf", &limit);
  p = 0.0;
  v = 0.0;
  a = 0.0;
  t = 0.0;

  while(t < limit) euler(dt);
  printf("%f %f %f %f\n", p, v, a, t);
  return 0;
}

haskell 代码:

import System.Environment (getArgs)

data EulerState = EulerState !Double !Double !Double !Double deriving(Show)
type EulerFunction = Double -> Double

main = do
  [dt, l] <- fmap (map read) getArgs
  print $ runEuler (EulerState 0 0 0 0) (**2) dt l

runEuler :: EulerState -> EulerFunction -> Double -> Double -> EulerState
runEuler s@(EulerState _ _ _ t) f dt limit = let s' = euler s f dt
                                             in case t `compare` limit of
                                                  LT -> s' `seq` runEuler s' f dt limit
                                                  _ -> s'

euler :: EulerState -> EulerFunction -> Double -> EulerState
euler (EulerState p v a t) f dt = (EulerState p' v' a' t')
    where t' = t + dt
          a' = f t'
          v' = v + a'*dt
          p' = p + v'*dt

Answer 1

关键点已经提过了by hammar和 by Philip JF .但是，让我收集它们并添加一些解释。

我会从上到下进行。

data EulerState = EulerState !Double !Double !Double !Double

您的类型具有严格的字段，因此每当该类型的对象被评估为 WHNF 时，其字段也被评估为 WHNF。在这种情况下，这意味着对象已被完全评估。这很好，但不幸的是，在我们的情况下还不够好。这种类型的对象仍然可以使用指向原始数据的指针来构造，而不是将原始数据解包到构造函数中，这就是加速度字段发生的情况(模循环不直接使用该类型，而是通过组件作为参数)。由于 euler 中未使用该字段， 你得到

Rec {
Main.$wrunEuler [Occ=LoopBreaker]
  :: GHC.Prim.Double#
     -> GHC.Prim.Double#
     -> GHC.Types.Double
     -> GHC.Prim.Double#
     -> Main.EulerFunction
     -> GHC.Prim.Double#
     -> GHC.Prim.Double#
     -> (# GHC.Types.Double,
           GHC.Types.Double,
           GHC.Types.Double,
           GHC.Types.Double #)

一个带有盒装参数的循环。这意味着在每次迭代中，一些 Double# s需要装箱，有的Double未装箱。装箱和拆箱不是非常昂贵的操作，但是在一个本来很紧的循环中，它们可能会消耗大量性能。相同装箱/拆箱问题的另一个实例与 EulerFunction 类型的参数有关。 ，稍后会详细介绍。 -funbox-strict-fields如suggested by Philp JF ，或 {-# UNPACK #-}至少加速场上的编译指示在这里有所帮助，但是只有在消除了功能评估的装箱/拆箱时，差异才变得相关。

print $ runEuler (EulerState 0 0 0 0) (**2) dt l

你正在通过 (** 2)这里作为一个论点。这与您在 C 中使用的函数不同，相应的 C 函数将是 return pow(t,2); ，并且使用我的 gcc，使用它几乎可以使 C 程序的运行时间增加一倍(不过，clang 没有区别)。最大的性能问题是 (**)是一个慢函数。由于(** 2)与 \x -> x*x 有不同的结果对于许多论点，没有重写规则，因此您确实可以使用 GHC 的 native 代码生成器获得该慢速功能(LLVM 似乎将其替换为 \x -> x*x 然而，由于两个 GHC 后端的巨大性能差异和 clang 结果)。如果您通过 (\x -> x*x)或 (^ 2)那里而不是 (** 2) ，你得到乘法(从 7.4 开始有 (^ 2) 的重写规则)。此时，在我的系统上，NCG 生成的代码和 LLVM 生成的代码的性能并没有太大的差异，但是 NCG 快了 10% 左右。

现在的大问题

runEuler :: EulerState -> EulerFunction -> Double -> Double -> EulerState
runEuler s@(EulerState _ _ _ t) f dt limit = let s' = euler s f dt
                                             in case t `compare` limit of
                                                  LT -> s' `seq` runEuler s' f dt limit
                                                  _ -> s'

runEuler是递归的，因此不能内联。这意味着传递的函数也不能在那里内联，dt和 limit每次迭代也会传递参数。函数不能被内联意味着在循环中，它的参数在传递给函数之前必须被装箱，然后它的结果必须被拆箱。那是昂贵的。这意味着在内联函数参数后无法进行任何优化。

如果您进行 worker /包装器转换和静态参数转换 suggested by hammar , runEuler可以内联，因此可以内联传递的函数，并且 - 在这种情况下 - 可以消除参数的装箱及其结果的拆箱。此外，甚至更大的影响，在这种情况下，可以消除函数调用并用一台机器操作代替。这导致了一个很好的紧密循环，如图所示

       174,208 bytes allocated in the heap
         3,800 bytes copied during GC

相比

16,000,174,912 bytes allocated in the heap
     1,475,432 bytes copied during GC

的原始。

总之，使用 native 代码生成器的 C 程序速度大约是 C 程序的一半，与我的机器上的 LLVM 后端的速度相同( native 代码生成器不是特别擅长优化循环，而 LLVM 是，因为循环非常在通过 LLVM 编译的许多语言中很常见)。