алгоритм — парное вычисление расстояния в переполнении стека

Преобразование из строки в удвоение затратно, переместите их за пределы цикла и кэшируйте их в отдельный кэш данных [],

typedef struct { double x, y, z; } position;
position datacache[data.size()];
for(int i = 0; i < data.size()-1; i++)
{
datacache[i].x = (double)atof(data[i][2].c_str());
datacache[i].y = (double)atof(data[i][3].c_str());
datacache[i].z = (double)atof(data[i][4].c_str());
}

Затем используйте datacache []. X, .y, .z в вашем цикле.

Вы можете использовать float вместо double, так как вы готовы приблизиться к 1% (таким образом, потеря дополнительных цифр точности, которые вы получаете от double, все еще в пределах 8-9 цифр точности).

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

Алгоритмическая оптимизация — разделите ваше трехмерное пространство на регионы, вычислите агрегирование по регионам, а затем агрегируйте эффекты от регионов.

1. Пересчитать все возможное вне цикла
2. Переместите преобразование из элементов массива из циклов в avoud, чтобы преобразовать любое значение более одного раза.

Код:

int colstart = 2;
int colend = 4;
double PE = 0;
double p_mass = 8.721e9 * 1.989e30; // mass of sim particle in kg
double mpc_to_m = 3.08567758e22; // meters per mpc
double G = 6.67384e-11; // grav. constant in mks units
double constVal= (-1 * G * pow(p_mass,2)) / mpc_to_m

double values[][3]= ... ;// allocate an array of size  data.size()*3
for(int i = 0; i < data.size()-1; i++){
value[i][0]=(double)atof(data[i][2].c_str());
value[i][1]=(double)atof(data[i][3].c_str());
value[i][2]=(double)atof(data[i][4].c_str());
}

// Calculating PE
for(int i = 0; i < data.size()-1; i++){
//cout << i << " " << data[i].size() << endl;
double xi=value[i][0]
double yi=value[i][1];
double yi=value[i][2];
for(int j = 0; j < i; j++){
double xDiff = xi - value[j][0] ;
double yDiff = yi - value[j][1] ;
double zDiff = zi - value[j][2];
PE += constVal / (xDiff*xDiff + yDiff*yDiff + zDiff*zDiff) ;
}
}

Но только профилирование может показать, если и насколько это увеличивает скорость.

Быстрый мультипольный метод была использована для ускорения, по-видимому, аналогичной задачи моделирования n-тел всех пар O (n ^ 2) с O (n). Я не знаком с деталями, за исключением того, что однажды увидел его в списке «Десять лучших алгоритмов, когда-либо изобретенных», но я думаю, что основная идея заключается в том, что большинство пар частиц представляют собой взаимодействия на больших расстояниях, которые являются слабыми и не слишком чувствительны к малым изменения в положении, что означает, что они могут быть хорошо аппроксимированы кластеризацией (я не уверен, как именно это делается), не теряя слишком много точности. Возможно, вы сможете применить аналогичные методы для вашей проблемы.

Вы можете использовать подход «разделяй и властвуй», как в случае с ближайшей парой: http://en.m.wikipedia.org/wiki/Closest_pair_of_points_problem. Вы также можете вычислить диаграмму Делоне или Вороного в O (n log n).

Одна вещь, которую вы должны сделать наверняка, по крайней мере, это переместить линии

   double x_i = (double)atof(data[i][2].c_str());
double y_i = (double)atof(data[i][3].c_str());
double z_i = (double)atof(data[i][4].c_str());

вне внутренней петли. Эти значения зависят только от i и не на jтак что вы определенно не хотите пересматривать их каждый раз. Затем есть несколько микрооптимизаций, которые могут заставить его работать немного быстрее. Наконец, если вы работаете на многопроцессорной машине, вы можете легко распараллелить это с openMP. Полуоптимизированная и распараллеленная версия кода может выглядеть так:

inline double squared(double x){
return x * x;
}

double compute_pe(vector<string *> data){
double PE = 0;
double p_mass = 8.721e9 * 1.989e30; // mass of sim particle in kg
double mpc_to_m = 3.08567758e22; // meters per mpc
double G = 6.67384e-11; // grav. constant in mks units

double PEarray[data.size()];

double numerator = (-1 * G * pow(p_mass,2))/ mpc_to_m;size_t i,j;

// Calculating PE
#pragma omp parallel for private(i, j)
for(i = 0; i < data.size()-1; i++) // -1 is for empty line at the end of every file
{
PEarray[i]=0;
double x_i = (double)atof(data[i][2].c_str());
double y_i = (double)atof(data[i][3].c_str());
double z_i = (double)atof(data[i][4].c_str());

//cout << i << " " << data[i].size() << endl;
for(j = 0; j < i; j++)
{

double x_j = (double)atof(data[j][2].c_str());
double y_j = (double)atof(data[j][3].c_str());
double z_j = (double)atof(data[j][4].c_str());

double dist_betw = sqrt(squared(x_i-x_j) + squared(y_i-y_j) + squared(z_i-z_j));

PEarray[i] += numerator / (dist_betw);
}
}

for(i = 0; i < data.size()-1; i++){
PE += PEarray[i];
}

return PE;
}

Вдобавок ко всем предыдущим советам, обычным делом является предварительное вычисление квадратов норм для всех векторов.

Учитывая, что || xy || ^ 2 = || x || ^ 2 + || y || ^ 2 — 2 * xy, можно избежать большого количества бесполезных умножений, если || x || ^ 2 вычисляется для всех х раз и навсегда в начале.

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