Оптимизация числодробилки

[al_pas]

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

Comments

[vitus_wagner]

Вообще такие вещи бы оптимизирующий компилятор должен делать.

[al_pas]

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

[v1adis1av]

А это проверялось? В смысле, если сравнить два расчёта с одинаковыми начальными условиями, один на программе, компилированной с оптимизацией, второй без, -- действительно появляется разница?

[al_pas]

Да, и очень существенная. У меня в программе используется коррекция, обеспечивающая сохранение полного момента и энергии с учетом потерь на излучение гравитационных волн. Шаг интегрирования тоже постоянно меняется, а про корректность алгоритма оптимизации gcc уже писал Линус Торвальдс.

[v1adis1av]

Линус вроде про последнюю версию gcc писал, а к более ранним претензий не возникало.

Я несколько раз свои числодробилки сравнивал с включенной и выключенной оптимизацией. Разница только в скорости счёта. Но у меня Монте-Карло, никакого долгого численного интегрирования.

[al_pas]

Попробую вкратце рассказать, как работает алгоритм. У меня используется модифицированный симплектический метод, но в качестве независимой переменной используется не время, а s=\int(\frac{A_1}{r_1}+\frac{A_2}{r_2})dt, где r_1 - расстояние между компонентами внутренней двойной (я считаю иерархическую систему), r_2 - расстояние третьего компонента от центра масс двойной, A_1 и A_2 - вычисляемые коэффициенты. Периодически вычисляется полная энергия и момент (которые являются суммой момента и энергии системы и момента и энергии, унесенных гравитационным излучением) и корректируются A_1 и A_2. Релятивистские термы считаются в приближении PN7/2, то есть: \frac{{{d^2}{\bf{x}}}}{{d{t^2}}} = \frac{m}{{{r^2}}}[{\bf{n}}( - 1 + {A_1} + {A_2} + {A_{5/2}} + {A_3} + {A_{7/2}}) + {\bf{v}}({B_1} + {B_2} + {B_{5/2}} + {B_3} + {B_{7/2}})]. Целые члены отвечают за поправки к потенциалу, полуцелые - за гравитационное излучение

[v1adis1av]

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

[al_pas]

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

[dims12]

А разве нельзя оставить только ту оптимизацию, которая ускоряет вычисление, но не влияет на точность?

В данном случае вообще непонятно, из-за чего происходит ускорение. Кажется, что деление на 2. и умножение на 0.5 должно занимать одинаковое время.

[al_pas]

Деление на 2 неудачный пример - тут действительно достаточно изменить на единицу экспоненту. Я убрал упоминание об этом из основного поста. Но для выражений типа a/(b/c+d/e) замена на a*c*e/(b*e+c*d) (при отключенной оптимизации) дает ускорение в три раза.

[al_pas]

Вот пример, когда оптимизация не спасает:

#include
#include
#include
using namespace std;

double secnds()
{
return clock()*1./CLOCKS_PER_SEC;
}

void div()
{
double a=1,b=2,c=3,d=4,e=5,f=0;
for (int i=0;i<1000;i++)
for (int j=0;j<1000;j++)
for (int k=0;k<1000;k++) {
a+=0.00000001;
b+=0.00000002;
c+=0.00000003;
d+=0.00000004;
e+=0.00000005;
f+=a/(b/c+d/e)+i+j+k;
}
printf("f = %21.13e\n",f);
}

void mul()
{
double a=1,b=2,c=3,d=4,e=5,f=0;
for (int i=0;i<1000;i++)
for (int j=0;j<1000;j++)
for (int k=0;k<1000;k++) {
a+=0.00000001;
b+=0.00000002;
c+=0.00000003;
d+=0.00000004;
e+=0.00000005;
f+=a*c*e/(b*e+c*d)+i+j+k;
}
printf("f = %21.13e\n",f);
}

int main(int argc, char *argv[])
{
double t = secnds();
mul();
cout << "mul, t = " << secnds() - t << "\n";
t = secnds();
div();
cout << "div, t = " << secnds() - t << "\n";
return 0;
}