我的循环可以再次优化吗？

Question

下面是我最里面的循环,运行几千次,输入大小为20 - 1000或更多.这段代码占用了99-99.5%的执行时间.我能做些什么来帮助挤出更多的表现吗？

我不打算将此代码移动到类似于使用树代码(Barnes-Hut)的东西,而是为了优化内部发生的实际计算,因为在Barnes-Hut算法中进行相同的计算.

任何帮助表示赞赏!

编辑:我在Core 2 Duo T5850(2.16 GHz)上使用Visual Studio 2008版本运行Windows 7 64位

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));
}

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]混淆,并且这可能在每次循环迭代中不必要地强制寻址p [i].

4A.尝试用if替换if(...)tau_q =

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

许多编译器认为mininum函数是他们可以用谓词操作而不是真正的分支来完成的,从而避免了管道刷新.

4B.[编辑分裂4a和4b分开]您可以考虑存储的tau _ .. Q2代替作为tau_q,和对R2进行比较1/V 2,而不是R2*R2/V2*V2.然后你避免在内循环中为每次迭代做两次乘法,换一个平方运算来计算最后的tau..q2.要做到这一点,分别收集tau_q1和tau_q2(未平方)的最小值,并取最小值那些结果在一个单一的标量操作的环路结束]

5. [编辑:我建议如下,但实际上它对OP的代码无效,因为他在循环中更新的方式.]将两个循环折叠在一起.使用两个循环和足够大的粒子集,您可以颠倒缓存并强制从第二个循环中的那些初始值的非缓存中重新获取.折叠是微不足道的.

除此之外,你需要考虑a)循环展开,b)矢量化(使用SIMD指令;手动编码汇编程序或使用英特尔编译器,这应该是相当不错的[但我没有经验],和c )去多核(使用OpenMP).

Answer 2

这条线real r3i = Constants::G * pow(r2, -1.5);会受到伤害.使用平方根的任何类型的sqrt查找或平台特定帮助都会有所帮助.

如果你有simd能力,分解你的矢量减去并平方到它自己的循环并立即计算它们将有所帮助.您的质量/混蛋计算也是如此.

想到的是 - 你是否保持足够的精确度？把东西带到第4个电源和第4根电源会通过欠/溢出混合器真正地破坏你的可用位.我确信你的答案确实是你的答案.

除此之外,它是一个数学繁重的功能,需要一些CPU时间.汇编程序对此的优化不会比编译器已经为您做的多得多.

另一种想法.由于这似乎与引力有关,有没有办法根据距离检查来剔除繁重的数学？基本上,半径/半径平方检查以对抗循环的O(n ^ 2)行为.如果你消除了你的粒子的1/2,那么它将快速运行x4.

最后一件事.您可以将内部循环线程连接到多个处理器.您必须为每个线程创建一个单独的内部版本,以防止数据争用和锁定开销,但是一旦每个线程完成,您可以计算每个结构的质量/加加速度值.我没有看到任何会阻止这种情况的依赖,但到目前为止我不是这方面的专家:)