Мой цикл может больше оптимизироваться?

Сначала простая вещь: переместите все «старые» переменные в другой массив. Вы никогда не обращаетесь к ним в своем основном цикле, поэтому вы касаетесь вдвое большего объема памяти, чем вам действительно нужно (и, таким образом, получаете вдвое больше промахов кеша). Вот недавнее сообщение в блоге на эту тему: http://msinilo.pl/blog/?p=614 . И, конечно же, вы можете выполнить предварительную выборку на несколько частиц впереди, например p [j + k], где k - некоторая константа, требующая некоторого экспериментирования.

Если вы также переместите массу, вы можете сохранить что-то вроде этого:

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

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

тогда обновление старых данных о частицах из новых данных становится одной большой memcpy от текущих частиц до старых частиц.

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

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

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

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

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

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

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

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

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

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

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

Вы можете использовать потоки это для более эффективного использования нескольких ядер

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

но действительно необходимо знаю, где это вам дороже всего

Вам следует повторно использовать действительные числа и векторы, которые вы всегда используете. Стоимость построения вектора или вещественного числа может быть тривиальной ... но не в том случае, если число 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;
        }
    }

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

a = b/c
d = e/f

на

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

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

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

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

Что настоящее ? Может это быть с плавающей точкой?

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

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

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

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

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

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

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

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

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

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

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

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

Все хорошее выше . Я делал аналогичные вещи с интегратором 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- sure-race-conditions/

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

Вам также может показаться интересной следующая колода, которую я дал в Seattle Code Camp:

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

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

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

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

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

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

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

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

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