2

我正在计算 C++ 中大量（~1e5）粒子的势能。为了做到这一点，我正在运行一个双循环，在其中我计算成对距离，并从这些距离计算系统的总势能。下面是相关的代码片段（它没有准备好复制/粘贴，因为需要定义数据，并且有些东西脱离了上下文；该方法仍然有效，这就是我要在这里展示的内容） ：

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

// Calculating PE
for(int i = 0; i < data.size()-1; i++) // -1 is for empty line at the end of every file
  {
    //cout << i << " " << data[i].size() << endl;
    for(int j = 0; j < i; j++)
     {

       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());

       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(pow((x_i-x_j),2) + pow(y_i-y_j,2) + pow(z_i-z_j,2)) * mpc_to_m;

       PE += (-1 * G * pow(p_mass,2)) / (dist_betw);  

     }
  }

有没有更快的方法来进行这种类型的计算？我对也涉及近似值的建议持开放态度，也就是说，如果它将总势能返回到大约 1% 左右。

谢谢！

4

7 回答 7

5

一些潜在的微观优化：

  • 乘法可能比使用pow()平方距离更快；
  • 公因数-1 * G * pow(p_mass,2) / mpc_to_m可以设为常数；
  • 求和可能会稍微快一点1.0 / dist_betw，然后乘以最后的公因数。
  • 您可能（或可能不）能够足够准确地近似平方根，比sqrt(). 有很多近似值可以尝试。
  • 将每个坐标转换为double一次，将值存储在另一个数组中，而不是在内循环的每次迭代中。

从算法上讲，您可以丢弃离当前点太远而无法对能量做出重大贡献的粒子。一个简单的修改可以在昂贵的平方根之前检查平方距离，如果它太大则继续。通过将粒子存储在空间感知数据结构（如Octree或只是一个简单的网格）中，您可能会获得进一步的改进，这样您甚至不需要查看远处的点；但如果粒子频繁移动，这可能不切实际。

于 2013-10-08T16:54:09.117 回答
4

从 string 到 double 的转换很昂贵，将它们移出循环并将它们缓存到单独的 datacache[] 中，

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。

您可以使用浮点数而不是双精度数，因为您愿意在 1% 以内进行近似（因此失去从双精度数获得的额外精度仍然在 8-9 位精度范围内）。

另一个效率改进 - 您也可以考虑使用定点整数算术（对于距离），您可以决定值的范围，而不是像浮点/双精度那样使用显式小数点存储，而是通过隐式缩放定点数值（这会影响距离计算。

算法优化 - 将您的 3D 空间分成多个区域，计算区域上的聚合，然后聚合来自区域的效果。

于 2013-10-08T17:00:11.890 回答
4
  1. 预先计算循环外的一切可能
  2. 将数组元素的转换移出循环，避免多次转换任何值。

编码：

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) ;  
     }
}

但只有分析才能显示这是否以及在多大程度上提高了速度。

于 2013-10-08T17:08:08.833 回答
2

您可以使用分而治之的方法，例如最近对问题：http ://en.m.wikipedia.org/wiki/Closest_pair_of_points_problem 。您还可以在 O(n log n) 中计算 delaunay 或 voronoi 图。

于 2013-10-08T17:03:16.887 回答
1

快速多极方法已用于加速看似与 O(n) 相似的 O(n^2) 全对 n 体模拟问题。除了在“曾经发明的十大算法”列表中看到它之外，我不熟悉细节，但我认为基本思想是大多数粒子对都是远程相互作用，它们很弱，对小粒子不太敏感位置的变化，这意味着它们可以通过聚类很好地近似（我不确定这是如何完成的）而不会失去太多的准确性。您也许可以将类似的技术应用于您的问题。

于 2013-10-08T17:07:57.157 回答
1

至少你应该做的一件事是移动线条

   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;
}
于 2013-10-08T17:23:52.030 回答
1

在所有先前的建议之上，一个常见的技巧是预先计算所有向量的平方范数。

鉴于||xy||^2 = ||x||^2 + ||y||^2 - 2*xy，如果||x||^2 计算所有x 一劳永逸地开始。

因此，这个成对距离问题变成了一个点积问题，这是一个基本的线性代数运算，可以由各种优化的库根据您的硬件进行计算。

于 2018-07-31T08:18:57.313 回答