c++ - 需要帮助优化代码(最小图像约定)

Question

我已经编写了一些模拟代码，并且正在使用“在 GDB 中随机中断”的调试方法。我发现我的程序 99.9% 的时间都花在了这个例程中(这是最小图像约定):

inline double distanceSqPeriodic(double const * const position1, double const * const position2, double boxWidth) {
        double xhw, yhw, zhw, x, y, z;                                                                                 

        xhw = boxWidth / 2.0;                                                                                      
        yhw = xhw;                                                                                      
        zhw = xhw;                                                                                      

        x = position2[0] - position1[0];                                                                               
        if (x > xhw) 
                x -= boxWidth;                                                                                     
        else if (x < -xhw)
                x += boxWidth;                                                                                     


        y = position2[1] - position1[1];                                                                               
        if (y > yhw) 
                y -= boxWidth;                                                                                     
        else if (y < -yhw)
                y += boxWidth;                                                                                     


        z = position2[2] - position1[2];                                                                               
        if (z > zhw) 
                z -= boxWidth;                                                                                     
        else if (z < -zhw)
                z += boxWidth;                                                                                     

        return x * x + y * y + z * z;                                                                                  
}

到目前为止我执行的优化(可能不是很重要):

• 返回距离的平方而不是平方根
• 内联
• 尽我所能
• 没有标准库膨胀
• 使用我能想到的每个 g++ 优化标志进行编译

我已经无能为力了。也许我可以使用花车而不是 double ，但我宁愿那是最后的选择。也许我可以以某种方式在这上面使用 SIMD，但我从来没有这样做过，所以我想这需要做很多工作。有什么想法吗？

谢谢

Answer 1

首先，您没有使用正确的算法。如果两个点之间的距离大于 boxWidth 怎么办？其次，如果你有多个粒子，调用一个函数来完成所有的距离计算并将结果放在输出缓冲区中会显着提高效率。内联有助于减少其中的一部分，但不是全部。任何预先计算(例如在您的算法中将框长度除以 2)都将在不需要时重复。

这里是一些用于计算的 SIMD 代码。您需要使用 -msse4 进行编译。在我的机器(macbook pro，llvm-gcc-4.2)上使用-O3，我的速度提高了大约 2 倍。这确实需要使用 32 位 float 而不是 double 算术。

SSE 实际上并没有那么复杂，只是看起来 很糟糕。例如不要写 a*b，你必须写笨重的 _mm_mul_ps(a,b)。

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <smmintrin.h>

// you can compile this code with -DDOUBLE to try using doubles vs. floats
// in the unoptimized code. The SSE code uses only floats.
#ifdef DOUBLE
typedef double real;
#else
typedef float real;
#endif

static inline __m128 loadFloat3(const float const* value) {
  // Load (x,y,z) into a SSE register, leaving the last entry
  // set to zero.
  __m128 x = _mm_load_ss(&value[0]);
  __m128 y = _mm_load_ss(&value[1]);
  __m128 z = _mm_load_ss(&value[2]);
  __m128 xy = _mm_movelh_ps(x, y);
  return _mm_shuffle_ps(xy, z, _MM_SHUFFLE(2, 0, 2, 0));
}

int fdistanceSqPeriodic(float* position1, float* position2, const float boxWidth,
                          float* out, const int n_points) {
  int i;
  __m128 r1, r2, r12, s12, r12_2, s, box, invBox;

  box = _mm_set1_ps(boxWidth);
  invBox = _mm_div_ps(_mm_set1_ps(1.0f), box);

  for (i = 0; i < n_points; i++) {
    r1 = loadFloat3(position1);
    r2 = loadFloat3(position1);
    r12 = _mm_sub_ps(r1, r2);

    s12 = _mm_mul_ps(r12, invBox);
    s12 = _mm_sub_ps(s12, _mm_round_ps(s12, _MM_FROUND_TO_NEAREST_INT));
    r12 = _mm_mul_ps(box, s12);

    r12_2 = _mm_mul_ps(r12, r12);

    // double horizontal add instruction accumulates the sum of
    // all four elements into each of the elements
    // (e.g. s.x = s.y = s.z = s.w =  r12_2.x + r12_2.y + r12_2.z + r12_2.w)
    s = _mm_hadd_ps(r12_2, r12_2);
    s = _mm_hadd_ps(s, s);

    _mm_store_ss(out++, s);
    position1 += 3;
    position2 += 3;
  }
  return 1;
}

inline real distanceSqPeriodic(real const * const position1, real const * const position2, real boxWidth) {
  real xhw, yhw, zhw, x, y, z;                                                                                 

  xhw = boxWidth / 2.0;                                                                                      
  yhw = xhw;                                                                                      
  zhw = xhw;                                                                                      

  x = position2[0] - position1[0];                                                                               
  if (x > xhw) 
    x -= boxWidth;                                                                                     
  else if (x < -xhw)
    x += boxWidth;                                                                                     


  y = position2[1] - position1[1];                                                                               
  if (y > yhw) 
    y -= boxWidth;                                                                                     
  else if (y < -yhw)
    y += boxWidth;                                                                                     


  z = position2[2] - position1[2];                                                                               
  if (z > zhw) 
    z -= boxWidth;                                                                                     
  else if (z < -zhw)
    z += boxWidth;                                                                                     

  return x * x + y * y + z * z;                                                                                  
}


int main(void) {
  real* position1;
  real* position2;
  real* output;
  int n_runs = 10000000;
  posix_memalign((void**) &position1, 16, n_runs*3*sizeof(real));
  posix_memalign((void**) &position2, 16, n_runs*3*sizeof(real));
  posix_memalign((void**) &output, 16, n_runs*sizeof(real));
  real boxWidth = 1.8;
  real result = 0;
  int i;
  clock_t t;

#ifdef OPT
  printf("Timing optimized SSE implementation\n");
#else
  printf("Timinig original implementation\n");
#endif

#ifdef DOUBLE
  printf("Using double precision\n");
#else
  printf("Using single precision\n");
#endif

  t = clock();
#ifdef OPT
  fdistanceSqPeriodic(position1, position2, boxWidth, output, n_runs);
#else
  for (i = 0; i < n_runs; i++) {
    *output = distanceSqPeriodic(position1, position2,  boxWidth);
    position1 += 3;
    position2 += 3;
    output++;
  }
#endif
  t = clock() - t;

  printf("It took me %d clicks (%f seconds).\n", (int) t, ((float)t)/CLOCKS_PER_SEC);
}