Haskell - 优化微分方程求解器

Question

我正在学习Haskell,并且正在尝试尽可能快地编写代码.在本练习中,我正在为一个简单的一维物理系统编写一个Euler积分器.

• C代码是用GCC 4.5.4和-O3.编译的.它运行在1.166秒.
• Haskell代码使用GHC 7.4.1和-O3.编译.它运行时间为21.3秒.
• 如果我编译Haskell -O3 -fllvm,它运行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

Hammar和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正如Philp JF所建议的那样,或{-# UNPACK #-}至少在加速度字段上的一个pragma在这里有所帮助,但只有当功能评估的装箱/拆箱也被消除时,差异才变得相关.

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的本机代码生成器获得了那么慢的功能(似乎LLVM取而代之\x -> x*x的是两个GHC后端的巨大性能差异和铿锵的结果).如果你通过(\x -> x*x)或(^ 2)那里而不是(** 2),你得到乘法((^ 2)从7.4开始有一个重写规则).此时,在我的系统上,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参数.函数不能内联意味着在循环中,它的参数必须在传递给函数之前被装箱,然后必须取消装箱其结果.那很贵.这意味着在内联函数参数之后不能进行任何优化.

如果你使用hammar建议的worker/wrapper转换和静态参数转换,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

原来的.

总之,使用本机代码生成器实现了C程序速度的一半,与我的盒子上的LLVM后端速度相同(本机代码生成器在优化循环方面不是特别好,而LLVM是,因为循环非常在通过LLVM编译的许多语言中很常见).

Answer 2

通过应用worker-wrapper转换,我获得了很好的提升runEuler.

runEuler :: EulerState -> EulerFunction -> Double -> Double -> EulerState
runEuler s f dt limit = go s
  where go s@(EulerState _ _ _ t) = if t < limit then go (euler s f dt) else s

这有助于f内联循环(这可能也发生在C版本中),摆脱了大量的开销.

Answer 3

我目前没有一个正常运行的LLVM,但是我在两倍之内

1. 使用-O2而不是-O3GHC(虽然我怀疑它很重要,但-O3没有维护)
2. 运用 -funbox-strict-fields
3. 使用x*x而不是x ** 2(与C代码相同)
4. 将您的"euler函数"移动为与C中相同的独立函数.

即

  func :: EulerFunction
  func x = x * x

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

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

你可以推进它(或者像Dons这样的Haskell性能专家会出现一个解决方案),我甚至没有看过这个生成的核心,但一般来说,使Haskell代码"像C一样快"的方法是"用C写它并使用FFI."