我的循环还能进行优化吗？

Question

以下是我的最内层循环，运行了数千次，输入大小为20-1000或更大。 这段代码占用了99-99.5％的执行时间。 有什么我可以做来帮助挤出更多性能吗？

我不想将此代码移动到使用树形代码（Barnes-Hut）之类的内容，而是朝着优化内部发生的实际计算的方向进行，因为Barnes-Hut算法中发生相同的计算。

感谢任何帮助！

编辑：我在Windows 7 64位上运行，使用Visual Studio 2008版在Core 2 Duo T5850（2.16 GHz）上运行。

typedef double real;

struct Particle
{
    Vector pos, vel, acc, jerk;
    Vector oldPos, oldVel, oldAcc, oldJerk;
    real mass;
};

class Vector
{
private:
    real vec[3];

public:
    // Operators defined here
};

real Gravity::interact(Particle *p, size_t numParticles)
{
    PROFILE_FUNC();
    real tau_q = 1e300;

    for (size_t i = 0; i < numParticles; i++)
    {
        p[i].jerk = 0;
        p[i].acc = 0;
    }

    for (size_t i = 0; i < numParticles; i++)
    {
        for (size_t j = i+1; j < numParticles; j++)
        {
            Vector r = p[j].pos - p[i].pos;
            Vector v = p[j].vel - p[i].vel;
            real r2 = lengthsq(r);
            real v2 = lengthsq(v);

            // Calculate inverse of |r|^3
            real r3i = Constants::G * pow(r2, -1.5);

            // da = r / |r|^3
            // dj = (v / |r|^3 - 3 * (r . v) * r / |r|^5
            Vector da = r * r3i;
            Vector dj = (v - r * (3 * dot(r, v) / r2)) * r3i;

            // Calculate new acceleration and jerk
            p[i].acc += da * p[j].mass;
            p[i].jerk += dj * p[j].mass;
            p[j].acc -= da * p[i].mass;
            p[j].jerk -= dj * p[i].mass;

            // Collision estimation
            // Metric 1) tau = |r|^2 / |a(j) - a(i)|
            // Metric 2) tau = |r|^4 / |v|^4
            real mij = p[i].mass + p[j].mass;
            real tau_est_q1 = r2 / (lengthsq(da) * mij * mij);
            real tau_est_q2 = (r2*r2) / (v2*v2);

            if (tau_est_q1 < tau_q)
                tau_q = tau_est_q1;
            if (tau_est_q2 < tau_q)
                tau_q = tau_est_q2;
        }
    }

    return sqrt(sqrt(tau_q));
}

- Mike Bailey

---

你知道什么是你最花费的东西吗？ - Keith Nicholas

@Keith：没有特别的，我的分析器只提供调用次数和时间。@Paul：我已经在帖子中添加了这些信息。 - Mike Bailey

有没有办法可以少调用它？如果你能少调用几次，它仍然会花费几乎所有的时间，但总时间会减少。 - Mike Dunlavey

我知道现在有点晚了，但是整数运算技巧怎么样？使用int64 vec [3]代替实数vec [3]。只有当您想要初始化模拟或显示当前迭代时，才需要将实数转换为整数向量数。决定您需要的小数位数。创建粒子时，只需将输入按1 ^ d缩放，其中d是小数位数。通过1 / 1 ^ d进行缩放以显示当前粒子。 - johnnycrash

@johnnycrash：我考虑过整数算术，但那会失去不可接受的大量精度。 - Mike Bailey

Answer 1

1. 将lengthsq()的调用内联。

2. 将pow(r2,-1.5)改为1/(r2*sqrt(r2))，以降低计算r^1.5的成本。

3. 在最内层循环中使用标量(p_i_acc等)，而不是p[i].acc来收集结果。编译器可能无法知道p[i]与p[j]没有别名(alias)，这可能会导致在每次循环迭代时需要寻址p[i]。

4a. 尝试用新的if(...) tau_q =替换原来的if(...) tau_q =

    tau_q=minimum(...,...)

许多编译器将最小函数识别为可以使用条件分支而不是真正的分支来执行的函数，从而避免了流水线刷新。

4b. [编辑以分离4a和4b] 您可能希望将 tau_..q2 存储为 tau_q，并与 r2/v2 进行比较，而不是 r2*r2/v2*v2。这样，您就可以避免在内部循环的每次迭代中进行两次乘法运算，而是换取单个平方运算，在结束时计算 tau..q2。为此，请分别收集 tau_q1 和 tau_q2（未平方）的最小值，并在完成循环时对这些结果进行单个标量操作。

5. [编辑：我建议以下内容，但实际上它对于 OP 的代码并不有效，因为他在循环中更新方式不同。] 将两个循环合并。如果使用两个循环和足够大的粒子集，则会使缓存失效，并从第二个循环中强制重新获取那些初始值。折叠是很容易做到的。

除此之外，您还需要考虑 a) 循环展开，b) 向量化（使用 SIMD 指令；手工编写汇编语言或使用英特尔编译器，据说非常擅长此项工作[但我没有经验]），以及 c) 多核（使用 OpenMP）。

- Ira Baxter

---

更好的选择是使用 r2*r2*inverse_sqrt(r2)（根据您所用的平台选择适当的函数），完全避免了除法。或者写成 r2*r2/sqrt(r2)，使用编译器的宽松数学设置，检查汇编输出或进行基准测试。 - Potatoswatter

@Potatoswatter: 不错的想法，但这是 pow(r2, +1.5) - 应该是 pow(r2, -1.5)。 - Paul R

@Paul，啊哦！嗯，那就用inverse_sqrt(r2)/r2吧。 - Potatoswatter

@Ira：什么是循环折叠？我做了一些谷歌搜索，但我找不到你在这里谈论的内容。 - Mike Bailey

@Sagekilla：我想我错过了你在循环中更新p[j].jerk的事实，因为我误读为p[i].jerk。如果是这样，你是正确的，你不能那样做。你仍然应该在循环中使用标量技巧。 - Ira Baxter

Answer 2

这行代码 real r3i = Constants::G * pow(r2, -1.5); 可能会导致性能问题。如果有sqrt查找或平台特定的平方根帮助，将有所帮助。

如果您具备simd能力，将向量减法和平方分解成自己的循环，并一次计算所有值，将有所帮助。质量/加速度计算同理。

另一个需要考虑的是-您的计算是否保持足够的精度？进行4次方和4次根运算会消耗大量位数。确保在完成时您的答案确实是您想要的答案。

除此之外，这是一个需要大量计算的函数，需要一些CPU时间。使用汇编优化不会比编译器为您做的更多。

另一个想法。由于这似乎与重力有关，是否有任何方法可以基于距离检查来消除繁重的数学计算？基本上，半径/半径平方检查以对抗您的循环的O(n^2)行为。如果您消除了一半的粒子，它将运行约4倍快。

最后一件事。您可以将内部循环线程化到多个处理器。您必须为每个线程创建一个单独的版本以防止数据争用和锁定开销，但是一旦每个线程完成，您就可以从每个结构中计算质量/加速度值。我没有看到任何会阻止此操作的依赖项，但在这个领域我绝不是专家 :)

- Michael Dorgan

---

是的，我的计算非常精确。tau_q的计算方式实际上是必要的，因为需要进行尺寸分析。tau_est_q1和tau_est_q2的计算方式是这样设计的，以便我得到正确的秒数单位的四次方，最终会转换成秒。 - Mike Bailey

粒子根据距离进行剔除，是否可行？ - Michael Dorgan

很遗憾，剔除是不可能的，因为这是一个纯重力模拟。所有成对交互都必须被考虑在内。我将尝试使用Barnes-Hut算法来实现O(n logn)，但我的主要问题仍然存在：循环内的计算是否仍然可以优化？ - Mike Bailey

Answer 3

首先，将所有“旧”的变量移动到不同的数组中。您在主循环中从未访问它们，因此实际上使用了两倍的内存（从而获得了两倍的缓存未命中）。以下是有关此主题的最近博客文章：http://msinilo.pl/blog/?p=614。当然，您可以提前预取一些粒子，例如 p[j+k]，其中 k 是需要进行一些实验的常数。

---

如果您还将质量移出，则可以这样存储：

struct ParticleData
{
    Vector pos, vel, acc, jerk;
};

ParticleData* currentParticles = ...
ParticleData* oldParticles = ...
real* masses = ...

然后，将旧粒子数据从新数据中更新变为将当前粒子数据一次性复制到旧粒子数据中。

---

如果您愿意让代码变得有点丑陋，通过以“转置”格式存储内容，可能可以获得更好的SIMD优化效果，例如

struct ParticleData
{
    // data_x[0] == pos.x, data_x[1] = vel.x, data_x[2] = acc.x, data_x[3] = jerk.x
    Vector4 data_x;

    // data_y[0] == pos.y, data_y[1] = vel.y, etc.
    Vector4 data_y;

    // data_z[0] == pos.z, data_y[1] = vel.z, etc.
    Vector4 data_z;
};

Vector4是一个单精度或两个双精度SIMD向量。在射线追踪中，这种格式很常见，因为它可以同时测试多条射线；它让你更有效地执行操作，如点积（无需洗牌），并且还意味着你的内存加载可以16字节对齐。不过，理解起来需要花费几分钟 :)

希望这能帮到你，如果你需要使用转置表示的参考资料，请告诉我（虽然我不确定它在这里能提供多少帮助）。

- celion

---

谢谢！这是一个有趣的想法，我得试一试。虽然我的代码没有明确说明，但是oldData在我的集成函数中被使用，我需要跟踪旧数据。我将尝试拆分我的数据并进行memcpy操作。 - Mike Bailey

memcpy是一个小优化，只有O（N）的复杂度。在内部循环中减少缓存未命中才是你将看到的主要收益（我希望）。 - celion

Answer 4

• 首先，您需要对代码进行性能剖析。这将取决于您正在运行的 CPU 和操作系统。

• 您可能要考虑是否可以使用浮点数而不是双精度浮点数。

• 如果您正在使用 gcc，则确保您正在使用 -O2 或可能的话使用 -O3。

• 您可能还想尝试一个好编译器，如英特尔的 ICC（假设正在运行在 x86 上？）。

• 再次假设这是（英特尔）x86，如果您拥有64位CPU，则请构建64位可执行文件（如果尚未）。额外的寄存器可以产生明显的差异（大约30%）。

- Paul R

---

我已经研究了使用浮点数的可能性，并决定通过提供真实的typedef来考虑它。在这种特殊情况下，我需要高速和精度，利用完整的机器精度（到1E-15）。 - Mike Bailey

作为一个快速测试 - 尝试使用 long double，以便与 double 的性能进行比较。你永远不知道，它可能会更快。 - Steve Jessop

Answer 5

如果这是用于视觉效果，并且您的粒子位置/速度只需要是近似的，那么您可以尝试用其相应泰勒级数的前几项替换sqrt。下一个未使用的项的大小表示您近似的误差范围。

- Emile Cormier

Answer 6

我的第一个建议是查看分子动力学文献，这个领域的人已经考虑了许多粒子系统的优化。例如，可以看一下 GROMACS。

对于许多粒子，最耗费时间的当然是双重循环。我不知道您需要多准确地计算粒子系统的时间演化，但如果您不需要非常精确的计算，可以简单地忽略太远的粒子之间的相互作用（必须设置截止距离）。一种非常有效的方法是使用带有缓冲区域的邻居列表，仅在需要更新这些列表时才进行更新。

- Adrien

Answer 7

以上都是好东西。我一直在做类似于二阶（Leapfrog）积分器的事情。在考虑了上述许多改进措施之后，我接下来做的两件事是开始使用SSE指令集以利用向量化，并使用一种新颖的算法并行化代码，该算法避免竞争条件并利用缓存局部性。

SSE示例：

http://bitbucket.org/ademiller/nbody/src/tip/NBody.DomainModel.Native/LeapfrogNativeIntegratorImpl.cpp

小说缓存算法，解释和示例代码：

http://software.intel.com/en-us/articles/a-cute-technique-for-avoiding-certain-race-conditions/

http://bitbucket.org/ademiller/nbody/src/tip/NBody.DomainModel.Native.Ppl/LeapfrogNativeParallelRecursiveIntegratorImpl.cpp

您可能会发现我在西雅图代码营上做的以下演示文稿也很有趣：

http://www.ademiller.com/blogs/tech/2010/04/seattle-code-camp/

你的第四阶积分器比较复杂，在双核系统上并行化的收益有限，但我肯定建议尝试使用SSE，我在这里得到了一些合理的性能提升。

- Ade Miller

---

非常感谢！那些链接非常有用。我已经有了另一个版本，它是并行化的，我循环了n^2次（而不是(n^2-n)/2），这样我就可以使用第j个粒子的力来累积第i个粒子。由于我的编写方式，该版本没有竞争条件，也不需要锁定。 - Mike Bailey

我刚刚看了你的NBody代码（你也在读《计算科学艺术》吗？），我有一个小建议：尝试将你的积分器编写为预测和校正步骤，并提取出力计算。然后，你可以将其编写为PEC方案（预测、评估、校正）。 - Mike Bailey

这就是我的编写方式。新算法使你能够回到（n ^ 2 - n）/2的并行代码，并且还能让你将问题块适合L1缓存，减少缓存未命中。这对我来说产生了显著的差异-接近2倍。在某个点上，缓存未命中成为远远超过您的最大开销的部分。 - Ade Miller

是的，我已经读了《计算机科学艺术》的部分内容。我应该重新开始阅读。目前所有的东西都使用全局时间步长。我应该再多尝试一些并转移到个体时间步长。我有点被CUDA这个快速但难以编码的东西所分心。我在C#中有一个Hermite积分器，也许我会尝试一下。谢谢！ - Ade Miller

Answer 8

你会从著名的 "快速反平方根算法" 受益吗？

float InvSqrt(float x)
{
    union {
        float f;
        int i;
    } tmp;
    tmp.f = x;
    tmp.i = 0x5f3759df - (tmp.i >> 1);
    float y = tmp.f;
    return y * (1.5f - 0.5f * x * y * y);
}

它返回1/r ** 2的相对准确表示（牛顿法的第一次迭代，具有聪明的初始猜测）。 它在计算机图形和游戏开发中被广泛使用。

- Escualo

---

这在几乎任何现有架构上都不是优化。如果您需要快速的近似倒数平方根，请使用专用硬件指令来获取它（例如，在具有SSE的x86上使用“rsqrtss”，在具有NEON的ARM上使用“vrsqrte.f32”，在具有AltiVec的PowerPC上使用“vrsqrtefp”）。它们既更快，又更准确。 - Stephen Canon

我可以验证Stephen所说的。我对invsqrt和sqrt然后除法进行了一些计时，结果invsqrt不仅速度略快，而且更不准确。为了达到类似的精度（通过添加更多的牛顿迭代步骤），invsqrt会更慢。 - Mike Bailey

Answer 9

除了简单的加减乘除，我在循环体内只看到pow()这个大型函数。它可能非常慢。你可以预先计算它、摆脱它或用更简单的东西替换它吗？

“real”是什么？它可以是浮点数吗？

除此之外，你将不得不转向MMX / SSE /汇编优化。

- AshleysBrain

---

在这段代码中，real 是 double 的同义词。请参见顶部的 typedef。 - kurige

Answer 10

同时考虑将常量::G的乘法操作移出循环。如果您可以更改存储向量的语义含义，使其有效地存储实际值/G，那么您可以根据需要进行万有引力常数的乘法运算。

任何可以减小Particle结构大小的操作都有助于提高缓存局部性。您似乎没有在这里使用old*成员。如果它们可以被删除，那可能会产生显着的差异。

考虑将我们的particle结构拆分为一对结构体。如果这样做，您第一次循环遍历数据以重置所有acc和jerk值的操作可以是一个高效的memset操作。然后，您将基本上拥有两个数组（或向量），其中每个数组的第n个粒子存储在各自数组的索引'n'处。

- StarLight