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

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 を使用します。

double の代わりに float を使用できます。これは、1% 以内で概算することをいとわないためです (したがって、double から得られる余分な桁数の精度を失うことは、依然として 8 ～ 9 桁の精度の範囲内です)。

もう 1 つの効率改善 - 値の範囲を決定する固定小数点整数演算 (距離用) を使用することを検討することもできます。値 (距離の計算から除外されます。

アルゴリズムの最適化 - 3D 空間を領域に分割し、領域全体の集計を計算してから、領域からの効果を集計します。

于 2013-10-08T17:00:11.890 に答える
1

少なくとも確実にすべきことの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 に答える