Можно ли оптимизировать мой цикл?

Question

Ниже приведен мой внутренний цикл, который выполняется несколько тысяч раз, с размерами входных данных от 20 до 1000 и более. Этот фрагмент кода занимает 99 - 99,5% времени выполнения. Есть ли что-нибудь, что я могу сделать, чтобы выжать из этого больше производительности?

Я не собираюсь перемещать этот код к чему-то вроде использования древовидных кодов (Barnes-Hut), но к оптимизации реальных вычислений, происходящих внутри, так как те же вычисления происходят в алгоритме Barnes-Hut.

Любая помощь приветствуется!

Редактировать: я использую 64-разрядную версию Windows 7 с выпуском Visual Studio 2008 на Core 2 Duo T5850 (2,16 ГГц)

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

Вы можете заменить любое вхождение:

a = b/c
d = e/f

с

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

, если вы знаете, что icf не приведет к выходу за пределы диапазона, и если ваше оборудование может выполнять 3 умножения быстрее, чем деление. Вероятно, не стоит объединять больше разделений, если у вас нет действительно старого оборудования с очень медленным разделением.

Вы избежите меньшего количества временных шагов, если будете использовать другие интеграционные схемы (например, Рунге-Кутта), но я подозреваю, что вы уже знаете это.

Answer 2

Вы должны повторно использовать реалы и векторы, которые вы всегда используете. Стоимость создания Vector или Real может быть тривиальной ... но не в том случае, если numParticles очень велика, особенно с вашим, казалось бы, циклом O ((n ^ 2) / 2).

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

Answer 3

То, что я ищу, это ветвление, они, как правило, убивают производительность.

Вы можете использовать разматывание петли.

также, запомните несколько с меньшими частями проблемы: -

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

- это почти то же самое, что выполнение всего цикла одним циклом, и вы можете получить ускорение за счет развертывания цикла и лучшего попадания в кэш

Вы можете сделать это, чтобы лучше использовать несколько ядер

у вас есть некоторые дорогостоящие вычисления, которые вы можете сократить, особенно если кальки в итоге вычисляют одно и то же, могут использовать кэширование и т. Д.

но на самом деле нужно знать, где это вам дороже всего

Answer 4

Да. Попробуйте посмотреть на вывод сборки. Это может дать подсказки относительно того, где компилятор делает это неправильно.

Теперь всегда всегда сначала применяйте оптимизацию алгоритма, и только если нет более быстрого алгоритма, вы должны выполнять частичную оптимизацию по сборке. А затем сначала выполните внутренние циклы.

Возможно, вы захотите получить профиль, чтобы увидеть, действительно ли это первое узкое место.