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

[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

[vitus_wagner]

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

[al-pas.livejournal.com]

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

[v1adis1av.livejournal.com]

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

[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]

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

[dims12.livejournal.com]

Странно. Этим должен заниматься компилятор, а не программист.

[dims12.livejournal.com]

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

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

[al-pas.livejournal.com]

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

[al-pas.livejournal.com]

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

#include <iostream>
#include <time.h>
#include <stdio.h>
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;
}

[al-pas.livejournal.com]

См. http://al-pas.livejournal.com/164026.html?thread=1045434#t1045434

[back-in-usa.livejournal.com]

Так ведь ещё на первом курсе учат, что деление должно быть последним действием.

[caztd.livejournal.com]

Деление ВСЕГДА значительно медленнее (~ в 10..20 раз) чем умножение:
https://stackoverflow.com/questions/4125033/floating-point-division-vs-floating-point-multiplication

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

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