C中代码的优化

Question

我正在尝试优化 C 中的代码，特别是占用总执行时间近 99.99% 的关键循环。这是那个循环：

#pragma omp parallel shared(NTOT,i) num_threads(4)
{
  # pragma omp for private(dx,dy,d,j,V,E,F,G) reduction(+:dU) nowait
  for(j = 1; j <= NTOT; j++){
    if(j == i) continue;
    dx = (X[j][0]-X[i][0])*a;
    dy = (X[j][1]-X[i][1])*a;
    d = sqrt(dx*dx+dy*dy);
    V = (D/(d*d*d))*(dS[0]*spin[2*j-2]+dS[1]*spin[2*j-1]);
    E = dS[0]*dx+dS[1]*dy;
    F = spin[2*j-2]*dx+spin[2*j-1]*dy;
    G = -3*(D/(d*d*d*d*d))*E*F;
    dU += (V+G);
  }
}

所有变量都是本地的。 NTOT=3600 的循环需要 0.7 秒，这是一个很长的时间，尤其是当我必须在整个程序中执行 500,000 次时，导致在此循环中花费了 97 小时。我的问题是这个循环中是否还有其他需要优化的地方？

我的电脑的处理器是 Intel core i5，带有 4 个 CPU(4X1600Mhz) 和 3072K L3 缓存。

Answer 1

针对硬件还是软件进行优化？

柔和：

摆脱除零等耗时的异常：

d = sqrt(dx*dx+dy*dy   + 0.001f   );
V = (D/(d*d*d))*(dS[0]*spin[2*j-2]+dS[1]*spin[2*j-1]);

您也可以试试 John Carmack、Terje Mathisen 和 Gary Tarolli 的"Fast inverse square root"

   D/(d*d*d)

部分。你也摆脱了分裂。

  float qrsqrt=q_rsqrt(dx*dx+dy*dy + easing);
  qrsqrt=qrsqrt*qrsqrt*qrsqrt * D;

牺牲一些精度。

还有一个部门也需要去掉：

(D/(d*d*d*d*d))

比如

 qrsqrt_to_the_power2 * qrsqrt_to_the_power3 * D

这里是快速逆sqrt：

float Q_rsqrt( float number )
{
    long i;
    float x2, y;
    const float threehalfs = 1.5F;

    x2 = number * 0.5F;
    y  = number;
    i  = * ( long * ) &y;                       // evil floating point bit level hacking
    i  = 0x5f3759df - ( i >> 1 );               // what ?
    y  = * ( float * ) &i;
    y  = y * ( threehalfs - ( x2 * y * y ) );   // 1st iteration
//  y  = y * ( threehalfs - ( x2 * y * y ) );   // 2nd iteration, this can be removed

    return y;
}

为了克服大数组的非缓存行为，您可以在较小的补丁/组中进行计算，尤其是在多对多 O(N*N) 算法的情况下。如：

get 256 particles.
compute 256 x 256 relations.
save 256 results on variables.

select another 256 particles as target(saving the first 256 group in place)
do same calculations but this time 1st group vs 2nd group.

save first 256 results again.

move to 3rd group

repeat.

do same until all particles are versused against first 256 particles.

Now get second group of 256. 

iterate until all 256's are complete.

您的 CPU 具有大缓存，因此您可以直接尝试 32k 粒子与 32k 粒子。但是 L1 可能不大，所以如果我是你，我会坚持使用 512 与 512（或 500 与 500 以避免缓存线 ---> 这将取决于架构）。

困难：

SSE、AVX、GPGPU、FPGA .....

正如@harold 评论的那样，SSE 应该是比较的起点，您应该向量化或至少通过 4-packed 向量指令并行化，这些指令具有最佳内存获取能力和流水线的优势。当您需要 3 到 10 倍的性能（在使用所有内核的 SSE 版本之上）时，您将需要一个符合 opencl/cuda 的 gpu（与 i5 同等价格）和 opencl（或 cuda）api，或者您也可以学习 opengl，但看起来更难（也许 directx 更容易）。

尝试 SSE 是最简单的，应该比我上面提到的快速反向快 3 倍。一个同等价格的 gpu 应该至少为数千个粒子提供另外 3 倍的 SSE。对于超过 100k 的粒子，当你对它进行足够的优化（使其对主存的依赖程度降低）时，对于这种算法，整个 gpu 可以实现单核 cpu 的 80 倍性能。 Opencl 提供了寻址缓存以保存数组的能力。因此，您可以在其中使用 TB/秒的带宽。

Answer 2

我会一直这样做random pausing 确定哪些线路成本最高。 然后，在修复某些东西后，我会再做一次，以找到另一个修复，依此类推。

也就是说，有些事情看起来很可疑。 人们会说编译器的优化器应该修复这些，但如果我能提供帮助，我从不依赖它。

• X[i]、X[j]、spin[2*j-1(and 2)] 看起来像是指针的候选对象。不需要做这个索引计算然后希望优化器可以去掉它。

• 你可以定义一个变量d2 = dx*dx+dy*dy，然后说d = sqrt(d2)。然后无论你有d*d，你都可以写d2。

• 我怀疑很多样本会落在sqrt 函数中，所以我会尝试想办法解决这个问题。

• 我确实想知道是否可以在此循环之外的单独展开循环中计算其中的一些量，例如 (dS[0]*spin[2*j-2]+dS[1]*spin[2*j-1])。在某些情况下，如果编译器可以保存一些寄存器，则两个循环可能比一个循环快。

Answer 3

我不敢相信 O(1) 循环的 3600 次迭代可能需要 0.7 秒。也许您的意思是具有 3600 * 3600 次迭代的双循环？否则我可以建议检查是否启用了优化，以及线程产生需要多长时间。

一般

您的内部循环非常简单，它只包含几个操作。请注意，除法和平方根大约比加法、减法和乘法慢 15-30 倍。你在做三个，所以大部分时间都被他们吃掉了。

首先，您可以在一次运算中计算平方根的倒数，而不是计算平方根，然后得到它的倒数。其次，您应该保存结果并在必要时重新使用它（现在您除以d 两次）。这将导致每次迭代产生一个有问题的操作，而不是三个。

invD = rsqrt(dx*dx+dy*dy);
V = (D * (invD*invD*invD))*(...);
...
G = -3*(D * (invD*invD*invD*invD*invD))*E*F;
dU += (V+G);

为了进一步减少rsqrt 花费的时间，我建议将其矢量化。我的意思是：使用 SSE 一次为两个或四个输入值计算 rsqrt。根据您的参数大小和所需的结果精度，您可以从this question 获取例程之一。请注意，它包含一个包含所有实现的小型 GitHub 项目的链接。

确实，您可以更进一步，使用 SSE（甚至 AVX）对整个循环进行矢量化，这并不难。

OpenCL

如果你准备使用一些大框架，那么我建议使用 OpenCL。您的循环非常简单，因此将它移植到 OpenCL 不会有任何问题（除了对 OpenCL 的一些初始适应）。

然后您可以使用 OpenCL 的 CPU 实现，例如来自英特尔或 AMD。他们俩都会自动使用多线程。此外，它们可能会自动矢量化您的循环（例如，参见 this article）。最后，如果您使用native_rsqrt 函数或类似的函数，他们有可能会自动找到rsqrt 的良好实现。

此外，您还可以在 GPU 上运行您的代码。如果你使用单精度，它可能会导致显着的加速。如果您使用双精度，那么就不太清楚了：现代消费级 GPU 在双精度下通常很慢，因为它们缺乏必要的硬件。

Answer 4

小优化：

• (d * d * d) 计算两次。存储 d*d 并将其用于 d^3 和 d^5
• 将 2 * x 修改为 x