代码之家 › 专栏 › 技术社区 › Mike Bailey

我的循环还能再优化吗?

physics loops optimization math c++

Mike Bailey · 技术社区 · 15 年前

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

我不打算把这段代码移到像使用树代码(Barnes-Hut)那样的地方,而是要优化内部发生的实际计算,因为在Barnes-Hut算法中也会发生同样的计算。

编辑:我运行在Windows7 64位的VisualStudio2008版本上的核心2双核T5850(2.16GHz)

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

14 回复 | 直到 14 年前

mana 15 年前

内联对lengthsq()的调用。
将pow(r2,-1.5)改为1/(r2*sqrt(r2)),以降低计算r^1.5的成本
在inner-most循环中使用标量(p\i\u acc等),而不是p[i].acc来收集结果。编译器可能不知道p[i]不是p[j]的别名,这可能会在每个循环迭代中不必要地强制p[i]的寻址。

4a级。尝试将if(…)tau_q=替换为

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

第4b条[编辑分开4a和4b分开]你可能会考虑存储Tauiq.Q2而不是Tauqq,并与R2/V2进行比较,而不是R2*R2/V2*V2。然后避免在内环中对每个迭代执行两次乘法,以换取在最后计算tau..q2的一次平方运算。为此,分别收集tau_q1和tau_q2的最小值(不是平方),并在循环完成时在单个标量运算中取这些结果的最小值]

[编辑:我建议如下,但事实上这对OP的代码无效,因为他在循环中更新的方式。]将两个循环折叠在一起。有了这两个循环和足够大的粒子集,就可以重击缓存,并强制从第二个循环中的那些初始值的非缓存中进行重取。折叠是琐碎的事。

除此之外,还需要考虑一个循环展开,b)矢量化(使用SIMD指令;要么手工编写汇编程序,要么使用英特尔编译器,这应该是相当不错的(但我没有经验),和c)多核(使用OpenMP)。

Michael Dorgan 15 年前

real r3i = Constants::G * pow(r2, -1.5); 会痛的。任何类型的sqrt查找或具有平方根的特定于平台的帮助都会有所帮助。

如果你有单指令多数据处理能力,把向量减法和平方分解成自己的循环,一次计算出来会有点帮助。质量/挺举计算也一样。

除此之外,它是一个数学繁重的函数,需要一些CPU时间。汇编程序对此的优化不会产生比编译器已经能为您做的更多的结果。

另一个想法。由于这似乎与重力有关,有没有办法根据距离检查剔除你繁重的数学题?基本上,一个radius/radius平方检查来对抗循环的O(n^2)行为。如果你去掉1/2的粒子,它的速度会快4倍。

最后一件事。您可以将内部循环线程连接到多个处理器。为了防止数据争用和锁定开销,您必须为每个线程创建一个单独的内部版本,但是一旦每个线程都完成了,您就可以统计每个结构的mass/jerk值。我没有看到任何依赖关系可以阻止这种情况,但我在这方面还不是专家:)

Paul R 15 年前

你可能会考虑是否可以使用浮动而不是双打。
如果您使用的是gcc,那么请确保您使用的是 -O2 -O3 .
再次假设这是(英特尔)x86,如果你有一个64位的CPU,然后建立一个64位的可执行文件,如果你还没有-额外的寄存器可以造成明显的差异(约30%)。

Emile Cormier 15 年前

如果这是为了视觉效果,而粒子位置/速度只需要近似值,则可以尝试替换 sqrt Taylor 系列。下一个未使用项的大小表示近似值的误差幅度。

celion 15 年前

首先,简单的事情是:将所有“旧”变量移到不同的数组中。您从未在主循环中访问它们,因此您所接触的内存是实际需要的两倍(从而导致两倍的缓存未命中)。以下是最近一篇关于这一主题的博文: http://msinilo.pl/blog/?p=614 prefetch 前面的一些粒子,例如p[j+k],其中k是需要一些实验的常数。

如果你把质量也移出去,你可以存储这样的东西:

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

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

然后,从新数据更新旧粒子数据就变成了从当前粒子到旧粒子的单个大内存。

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字节对齐的。不过,你肯定要花几分钟的时间才能清醒过来:)

希望有帮助,让我知道如果你需要一个关于使用转置表示的参考(虽然我也不知道它在这里到底有多少帮助)。

Adrien 15 年前

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

对于许多粒子,杀死你的当然是双倍粒子 for 循环。我不知道你需要多精确地计算粒子系统的时间演化,但是如果你不需要非常精确的计算,你可以忽略相距太远的粒子之间的相互作用(你必须设置一个截止距离)。一个非常有效的方法是使用 neighbour lists 只有在需要时才使用缓冲区来更新这些列表。

Ade Miller 15 年前

上面都是好东西。我一直在做类似的事情,以二阶(蛙跳)积分。在考虑了上面建议的许多改进之后,接下来我做的两件事是开始使用SSE intrinsics来利用矢量化,并使用一种新的算法来并行化代码,该算法避免了竞争条件,并利用了缓存局部性。

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,我在这里得到了一些合理的性能改进。

AshleysBrain 15 年前

pow() 是我在循环体中看到的唯一重量级函数。可能很慢。你能预先计算它或者把它去掉,或者用更简单的东西代替它吗?

什么 real

除此之外,您还必须转向MMX/SSE/装配优化。

Escualo 15 年前

你会从著名的 fast inverse square root" algorithm ?

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的一个相当精确的表示(牛顿方法的第一次迭代,带有一个巧妙的初始猜测)。它被广泛应用于计算机图形学和游戏开发。

StarLight 15 年前

任何修剪粒子结构大小的操作都将有助于改进缓存局部性。你好像没有用这里的老会员。如果它们能被移除,这将有可能产生重大的影响。

Joshua 15 年前

您可能想先分析一下这是否真的是瓶颈。

Keith Nicholas 15 年前

我寻找的是分支,它们往往是性能杀手。

可以使用循环展开。

另外,请记住问题的多个小部分:-

  for (size_t i = 0; i < numParticles; i++)
    {
        for (size_t j = i+1; j < numParticles; j++)
        {

这和让一个循环做任何事情差不多,你可以通过循环展开和更好地命中缓存来提高速度

您可以通过线程来更好地利用多核

你有一些昂贵的计算,你可能可以减少,特别是如果计算结束计算相同的东西,可以使用缓存等。。。。

但我真的需要知道你在哪里花费最多

Puppy 15 年前

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

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

            // da = r / |r|^3
            // dj = (v / |r|^3 - 3 * (r . v) * r / |r|^5
            da = r * r3i;
            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
            mij = p[i].mass + p[j].mass;
            tau_est_q1 = r2 / (lengthsq(da) * mij * mij);
            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;
        }
    }

sigfpe 15 年前

a = b/c
d = e/f

icf = 1/(c*f)
a = bf*icf
d = ec*icf

如果您使用其他集成方案(例如Runge-Kutta),您将获得更少的时间步,但我怀疑您已经知道这一点。