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

[apximhd]

Рабочая программа для моделирования эволюции тройной звездной системы с учетом общей теории относительности непрерывно улучшалась и модифицировалась в течение последних десяти лет. В результате ее код превратился в такую кашу, что дальнейшая его поддержка стала невозможной.
Пришлось заняться рефакторингом, а заодно и оптимизацией.
В итоге код сократился в полтора раза, а быстродействие выросло в два раза. Интересно, что основной прирост быстродействия удалось получить, заменив везде где только можно деление умножением. Например, простая замена выражения a/(b/c+d/e) выражением a*c*e/(b*e+c*d) дает выигрыш в скорости выполнения в три раза.

Upd.: оптимизация компилятора отключена по причинам, которые обсуждаются далее в комментариях.

Upd2.: впрочем, для случая замены выражения a/(b/c+d/e) выражением a*c*e/(b*e+c*d) включение или выключение оптимизации не играет роли.

Comments

[al-pas.livejournal.com]

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

[v1adis1av.livejournal.com]

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

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

[al-pas.livejournal.com]

Попробую вкратце рассказать, как работает алгоритм. У меня используется модифицированный симплектический метод, но в качестве независимой переменной используется не время, а 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}})]. Целые члены отвечают за поправки к потенциалу, полуцелые — за гравитационное излучение.
Включение оптимизации приводит к тому, что итоговый результат отличается. Понятно, что из-за потерь точности, результат после 10 миллионов оборотов всё равно не соответствует точному, но при отключенной оптимизации я, по крайней мере, способен контролировать код, а копаться каждый раз в ассемблерном коде, порожденном оптимизатором, у меня просто нет времени и сил.

[v1adis1av.livejournal.com]

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

[al-pas.livejournal.com]

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