4

我有一段代码可以从“double **times”计算一个值。假设“时间”的维度为 [nsims][N]（使用 malloc.. 创建），其中 int N=40 和 int nsims=50000。

结果存储在“double **moments”中。所以我们有 5 个嵌套的 for 循环。

然而问题在于速度，因为这段代码需要运行大约 100 万次。

我已经在使用线程（此处未显示）将最里面的 for 循环拆分为 10 个并行线程，这已经节省了大量时间。

有没有人看到其他优化的可能性，特别是关于不同的数据结构或类似的东西？

即使我没有“interm = ...”公式，它仍然需要太多时间。

for(j=2;j<=N;j++) {     
    for(k=j;k<=N;k++) {
        moment=0;
        for(i=2;i<=N;i++) {
            for(l=i;l<=N;l++) {
                if(strcmp(mmethod, "emp")==0) {
                    for(a=0;a<nsims;a++) {
                        interm=interm + (double) times[a][k] *
                                        times[a][j]*times[a][i] *
                                        times[a][l];    
                    }
                    interm = (double) interm/nsims;
                    moment = moment + (interm*i*l);
                    interm=0;
                }
            }
        }
        if(!(changed_times[k]==0
             && changed_times[j]==0
             && changed_times[l]==0
             && changed_times[i]==0))
        {
            moments[0][pcount]=(double) moment;
        } else {
            moments[0][pcount]=moments[0][pcount];
        }
        pcount++;
    }
}
4

3 回答 3

3

请注意，在您的内部循环中，您每次都在查找和相乘times[a][k]*times[a][j]*times[a][i]，即使该表达式对于 . 的每个值都是相同的a。对于乘法和内存查找来说，它可能很昂贵。（也许编译器足够聪明，可以将其优化掉，我不知道。）您可以尝试在内部循环中缓存这些值，如下所示：

  ...
  double akji[nsims];
  for (a = 0; a < nsims; ++a) { akji[a] = times[a][k]*times[a][j]*times[a][i]; }
  for(l=i;l<=N;l++) {
interm=0;
for(a=0;a<nsims;a++) {
  interm += akji[a]*times[a][l]; 
}
moment += (interm*l);
  }
  moment = moment * i / nsims;
  ...
于 2013-03-15T14:14:45.990 回答
1

我想应该从对问题的更高层次的描述开始。

但作为次要选项，我建议交换数组索引，以便更轻松地编写一个结合了四个（可能是不同的）向量的超快 SSE 内循环：

 double times[N+1][nsims], *tend = times[N+1];
 double *j,*k,*i,*l;
 for (j=times[2];j<tend;j+=nsims)
  for (k=j;k<tend;k+=nsims)
   if (strcmp( )) ... /* One _can_ move this elsewhere, but why bother? */
   for (i=times[2];i<tend;i+=nsims)
    for (l=i;l<tend;l+=nsims) {
      interm = efficient_sse_implementation(j,k,i,l, nsims);
      ...
    }

对于少于 4 个不同数组的情况，也可以通过编写不同的内核来实现微小的优化。（在这种情况下，可以跳过每步一个内存操作。）

编辑

在这种情况下，模式的结构for(j=2;j<=N;j++) for (k=j;k<=N;k++)重复了两次，仅此一项就暗示了更高级别优化的可能性——执行的操作是什么？虽然为此苦苦挣扎，但这种模式仍然提出了另一种方法：缓存 780 个（左右）子产品，但同时执行循环阻塞。这种方法不应该与我对先生的评论有同样的问题。格贝尼森

 for (A=0;A<50000;A+=100) {
    int k=0;
    for (i=2;i<=N;i++)
     for (j=i;j<=N;j++,k++)
       for (a=0;a<100;a++) precalc[k][a]=times[i][A+a]*times[j][A+a];

    for (i=0;i<k;i++)  // Now i loops from 0..779 or so
     for (j=0;j<k;j++) {
       for (a=0;a<100;a++) partial_product+=precalc[i][a]*precalc[j][a];
       // accumulate also partial_product to moment
     }
 }

免责声明：这是未经尝试的，但存在一些最佳的块大小（不一定是 100）（它可能比以前的情况更糟）。另请注意，这种方法会为预先计算的表使用大量内存。（选择 100 的块大小需要 624000 字节的内存，这听起来不错。要低于 256k，块长度只能是 42）。

编辑 2：

// 请注意，EDIT_1 中的循环同时计算P[2][a]*P[3][a]和P[3][a]*P[2][a]。

    for (i=0;i<k;i++)  // Now i loops from 0..779 or so, but... we can limit the
     for (j=i;j<k;j++) { // calculation to the upper triangle of the 780^2 matrix
       for (a=0;a<100;a++) partial_product+=precalc[i][a]*precalc[j][a];
       moment[i]+=partial_product;
       moment[lower_triangle(i)]+=partial_product;  // <-- 50% speed increase
     }

编辑3：这里有一些尝试：

gcc -O4 -DCACHELEVEL=2 -DPOPULATE=1 -DARRAY_OPT=1 && time ./a.out
  • POPULATE初始化数组（假设非零内容很重要）
  • ARRAY_OPT=1将数组索引切换到（也许）更好的顺序
  • CACHELEVEL=2或 3 次切换缓存中间结果。
  • STRCMP可以在源代码中找到测试 memcmp vs. strcmp vs '1'

NOT TODO 1 : LOOP_BLOCKING with cached values - 降低性能
TODO 2 : 只计算上三角形
TODO 3 : 找出changed_times[n]and的含义moments[0][p]
- 正如它现在所表明的那样，没有任何计算被保存！

#include <stddef.h>
#define N 40
#define nsims 8000

#if ARRAY_OPT
#define TIMES(n,a) times[n][a]
double times[N+1][nsims]; // [nsims];
#else
#define TIMES(n,a) times[a][n]
double times[nsims][N+1];
#endif

#define STRCMP 1
// vs.
// #define STRCMP1 strcmp(mmethod, "emp")==0

void init()
{
#ifdef POPULATE
    int i,a;
    for (i=0;i<=N;i++)
      for (a=0;a<nsims;a++)
         TIMES(i,a) = (double)((i^a)&7) - 3.5;
#endif
}

  double moments[4000] = { 0 };
  double cache1[nsims];
  double cache2[nsims];

int main()
{
  int j,k,i,l,a, pcount=0;
  init();
  int changed_times[N+1]={0};
  char *mmethod="emp";

  double moment,interm;
  for(j=2;j<=N;j++) {     
    for(k=j;k<=N;k++) {
#if CACHELEVEL == 2
        for (a=0;a<nsims;a++) cache1[a]=TIMES(j,a)*TIMES(k,a);
#endif
        moment=0;
        for(i=2;i<=N;i++) {
#if CACHELEVEL == 3
            for (a=0;a<nsims;a++)         cache2[a]=TIMES(j,a)*TIMES(k,a)*TIMES(i,a);
#else
            for (a=0;a<nsims;a++) cache2[a]=cache1[a]*TIMES(i,a);
#endif
            for(l=i;l<=N;l++) {
                if(STRCMP) { 
                    for(a=0;a<nsims;a++) {
#if CACHELEVEL >= 2
                        interm += (double) cache2[a]*TIMES(l,a);
#else
                        interm=interm + (double) TIMES(k,a) * TIMES(j,a) * TIMES(i,a) * TIMES(l,a);
#endif
                    }
                    interm = (double) interm/(double)nsims;
                    moment = moment + (interm*i*l);
                    interm=0;
                }
            }
        }
        //if(!(changed_times[k]==0
        //     && changed_times[j]==0
        //     && changed_times[l]==0
        //     && changed_times[i]==0))
        //{
        //    moments[0][pcount]=(double) moment;
        //      changed_times[k]++;changed_times[j]++; /* or what? */
        //      changed_times[l]++;changed_times[i]++;
        //} else {
        //    moments[0][pcount]=moments[0][pcount];
        //}
        pcount++;
    }
  }
  printf("%d %f\n",pcount, moment);
}
于 2013-03-15T14:02:40.763 回答
0

第一个明显的优化是移出strcmp()循环。

字符串比较可能需要相当长的时间（实际上并不多，但是如此多次重复此调用会产生很大的不同）。此外，编译器可能永远不会优化此调用，而其结果在整个处理过程中是恒定的。因此，在进入嵌套循环之前将结果存储在一个临时布尔变量中，并且只测试循环内的布尔值。

此外，在尝试优化一段代码时，请确保使用发布目标（没有调试信息）进行编译并打开所有可能的编译器优化。

于 2013-03-15T14:10:50.120 回答