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

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

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] на каждой итерации цикла.

4а. Попробуйте заменить if (...) tau_q = на

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

Многие компиляторы распознают функцию mininum как функцию, которую они могут выполнять с предикатными операциями, а не с реальными ветвями, избегая сбрасывания конвейера.

4b. [РЕДАКТИРОВАТЬ, чтобы разделить 4a и 4b на части]. Вы можете вместо этого хранить tau _ .. q2 как tau_q и сравнивать с r2 / v2, а не с r2 * r2 / v2 * v2. Тогда вы избегаете делать два умножения для каждой итерации во внутреннем цикле, в обмен на одну операцию возведения в квадрат для вычисления tau..q2 в конце. Чтобы сделать это, соберите минимумы tau_q1 и tau_q2 (не в квадрате) отдельно и возьмите минимум этих результатов в одной скалярной операции по завершении цикла]

5. [РЕДАКТИРОВАТЬ: я предложил следующее, но на самом деле это не подходит для кода OP, из-за способа, которым он обновляет в цикле.] Сложите два цикла вместе. Имея два цикла и достаточно большой набор частиц, вы очищаете кэш и производите принудительное повторное извлечение из не кеша этих начальных значений во втором цикле. Сгиб тривиально сделать.

Помимо этого, вам необходимо рассмотреть а) развертывание цикла, б) векторизацию (с использованием инструкций SIMD; либо на ассемблере ручного кодирования, либо с помощью компилятора Intel, что должно быть довольно неплохо в этом [но у меня нет опыта работы с ним] и c) многоядерность (с использованием OpenMP).

Answer 2

Эта строка real r3i = Constants::G * pow(r2, -1.5); будет больно. Любой вид справки sqrt или конкретной платформы с квадратным корнем поможет.

Если у вас есть способности simd, то вам поможет разбивка векторных вычитаний и квадратов в отдельный цикл и их одновременное вычисление. То же самое для вашей массы / рывков.

Что-то, что приходит на ум, - вы держите достаточно точности с вашим кальцием? Подводя вещи к 4-й степени и 4-му корню, вы получаете доступ к битам через блендер Under / Overflow. Я был бы уверен, что ваш ответ действительно будет вашим ответом, когда он будет завершен.

Кроме того, это сложная математическая функция, которая потребует некоторого времени процессора. Оптимизация ассемблером этого не даст слишком много, чем компилятор уже может сделать для вас.

Еще одна мысль. Поскольку это связано с гравитацией, есть ли способ отбросить вашу тяжелую математику на основе проверки расстояния? По сути, проверка радиуса / радиуса в квадрате для борьбы с O (n ^ 2) поведением вашей петли. Если вы удалите половину своих частиц, они будут бегать примерно в 4 раза быстрее.

И последнее. Вы можете подключить свой внутренний цикл к нескольким процессорам. Вы должны были бы сделать отдельную версию своих внутренних компонентов для каждого потока, чтобы предотвратить конфликт данных и накладные расходы на блокировку, но как только каждый поток был завершен, вы могли бы подсчитать значения массы / рывка из каждой структуры. Я не видел никаких зависимостей, которые могли бы помешать этому, но я далеко не эксперт в этой области:)

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 байт. Хотя определенно требуется несколько минут, чтобы обернуть голову:)

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

Answer 4

Если это для визуальных эффектов, и ваше положение / скорость частиц должны быть только приблизительными, тогда вы можете попробовать заменить sqrt на первые несколько членов соответствующей серии Taylor . Величина следующего неиспользованного члена представляет предел погрешности вашего приближения.

Answer 5

• Сначала вам нужно профилировать код. Метод для этого будет зависеть от того, какой процессор и ОС вы используете.

• Вы можете подумать, можете ли вы использовать числа с плавающей запятой вместо двойных.

• Если вы используете gcc, убедитесь, что вы используете -O2 или, возможно, -O3.

• Возможно, вы захотите попробовать хороший компилятор, такой как Intel ICC (при условии, что он работает на x86?).

• Опять же, если предположить, что это (Intel) x86, если у вас 64-битный процессор, то создайте 64-битный исполняемый файл, если вы этого еще не сделали - дополнительные регистры могут внести заметные изменения (около 30%) .

Answer 6

Все хорошие вещи выше. Я делал подобные вещи с интегратором 2-го порядка (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, у меня здесь есть некоторые разумные улучшения производительности.

Answer 7

Мой первый совет - взглянуть на литературу по молекулярной динамике. Люди в этой области рассмотрели множество оптимизаций в области систем частиц. Взгляните, например, на GROMACS .

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

Answer 8

Рассмотрим также вытащить умножение констант :: G из цикла.Если вы можете изменить семантическое значение хранимых векторов, чтобы они эффективно сохраняли фактическое значение / G, вы можете сделать мультипликацию гравитационной постоянной по мере необходимости.

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

Рассмотрите возможность разбиения нашей структуры частиц на пару структур.Ваш первый цикл по данным для сброса всех значений acc и jerk может быть эффективным memset, если вы сделаете это.Тогда у вас будет по существу два массива (или вектора), в которых частица n находится в индексе n каждого из массивов.

Answer 9

Не могли бы вы воспользоваться известным алгоритмом " быстрый обратный квадратный корень" ?

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 (первая итерация метода Ньютона с умным начальным предположением). Он широко используется для компьютерной графики и разработки игр.

Answer 10

Помимо простого сложения / вычитания / деления / умножения, pow() - единственная тяжелая функция, которую я вижу в теле цикла Это, вероятно, довольно медленно. Можете ли вы предварительно вычислить его, или избавиться от него, или заменить его чем-то более простым?

Что такое real? Это может быть поплавок?

Кроме того, вам придется обратиться к оптимизации MMX / SSE / сборки.