我应该如何使用递归扫描算法实现 Cilk 并行性？答案

【问题标题】：How Should I Implement Cilk Parallelism with a Recursive Scan Algorithm?我应该如何使用递归扫描算法实现 Cilk 并行性？
【发布时间】：2014-07-16 19:55:00
【问题描述】：

我实现了一个递归扫描（前缀总和）算法，我将其包含在下面。在这里，我简单地生成大小为 2 到 27 次方的随机列表，检查简单的顺序扫描的准确性。它有效。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <mkl.h>

int *pscan(int *x, int n, int z, int chunk_size);
int reduce(int *x, int n);

int main(int argc, char **argv)
{
  int n;
  int i, j, k;
  int *x, *seq, *r;
  double begin, end;

  srand48(time(0));

  /* Randomly generate array of size n. */

  for (k = 2; k < 28; k++) {
    n = (int) pow(2, k);

    seq = (int *) malloc(sizeof(int) * n);
    x = (int *) malloc(sizeof(int) * n);

    for (i = 0; i < n; i++) {
      x[i] = lrand48() % 100 - 50;
      seq[i] = x[i];
    }

    /* Parallel scan. */

    begin = dsecnd();

    r = pscan(x, n, 0, 2);

    end = dsecnd();

    printf("%d %lf\n", n, end - begin);

    /* Sequential check. */

    for (i = 1; i < n; i++) {
      seq[i] = seq[i - 1] + seq[i];
    }

    for (i = 0; i < n; i++) {
      if (r[i] != seq[i]) {
        fprintf(stderr, "AGGGHHH!!! ERROR.  Found with vector: \n");
        for (j = 0; j < n; j++) {
          printf("%d ", x[i]);
        }
        printf("\n");
        exit(1);
      }
    }
    free(r);
    free(x);
    free(seq);
  }

  return 0;
}

/* Perform parallel scan. */
int *pscan(int *x, int n, int z, int chunk_size)
{
  int i, j;
  int *sums, *sumscan, *scan, **fsum, *rv;

  /* Base case, serially scan a chunk. */
  if (n <= chunk_size) {
    scan = (int *) malloc(sizeof(int) * n);

    scan[0] = x[0] + z;
    for (i = 1; i < n; i++) {
      scan[i] = x[i] + scan[i - 1];
    }

    return scan;
  }

  sums = (int *) malloc(sizeof(int) * (n / chunk_size));

  /* Reduce each chunk of the array. */
  for (i = 0; i < n / chunk_size; i++) {
    sums[i] = reduce(&x[i * chunk_size], chunk_size);
  }

  /* Perform a scan on the sums. */
  sumscan = pscan(sums, n / chunk_size, 0, chunk_size);

  free(sums);

  fsum = (int **) malloc(sizeof(int *) * (n / chunk_size));

  /* Perform a recursive scan on each chunk, using
     the appropriate offset from the sums scan. */
  for (i = 0; i < n / chunk_size; i++) {
    if (i > 0) {
      fsum[i] = pscan(&x[i * chunk_size], chunk_size, sumscan[i - 1], chunk_size);
    } else {
      fsum[i] = pscan(&x[i * chunk_size], chunk_size, 0, chunk_size);
    }
  }

  free(sumscan);

  rv = (int *) malloc(sizeof(int) * n);

  /* Join the arrays. */
  for (i = 0; i < n / chunk_size; i++) {
    for (j = 0; j < chunk_size; j++) {
      rv[i * chunk_size + j] = fsum[i][j];
    }
  }

  for (i = 0; i < n / chunk_size; i++) {
    free(fsum[i]);
  }

  free(fsum);

  return rv;
}

/* Serial reduction. */
int reduce(int *x, int n)
{
  int i;
  int sum;

  sum = 0;

  for (i = 0; i < n; i++) {
    sum += x[i];
  }

  return sum;
}

现在，我想并行化它。因为我觉得有点时髦，所以我破解了 Cilk 的实现。我只是替换了两个主要的 for 循环来并行化 1）减少和 2）每个块的递归扫描，使用块减少的适当扫描作为偏移量。看起来是这样的。

#include <stdio.h>
#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <math.h>
#include <cilk/cilk.h>
#include <mkl.h>

int *pscan(int *x, int n, int z, int chunk_size);
int reduce(int *x, int n);

int main(int argc, char **argv)
{
  int n;
  int i, j, k;
  int *x, *seq, *r;
  double begin, end;

  srand48(time(0));

  /* Randomly generate array of size n. */

  for (k = 2; k < 28; k++) {
    n = (int) pow(2, k);

    seq = (int *) malloc(sizeof(int) * n);
    x = (int *) malloc(sizeof(int) * n);

    for (i = 0; i < n; i++) {
      x[i] = lrand48() % 100 - 50;
      seq[i] = x[i];
    }

    /* Parallel scan. */

    begin = dsecnd();

    r = pscan(x, n, 0, 2);

    end = dsecnd();

    printf("%d %lf\n", n, end - begin);

    /* Sequential check. */

    for (i = 1; i < n; i++) {
      seq[i] = seq[i - 1] + seq[i];
    }

    for (i = 0; i < n; i++) {
      if (r[i] != seq[i]) {
        fprintf(stderr, "AGGGHHH!!! ERROR.  Found with vector: \n");
        for (j = 0; j < n; j++) {
          printf("%d ", x[i]);
        }
        printf("\n");
        exit(1);
      }
    }
    free(r);
    free(x);
    free(seq);
  }

  return 0;
}

/* Perform parallel scan. */
int *pscan(int *x, int n, int z, int chunk_size)
{
  int i, j;
  int *sums, *sumscan, *scan, **fsum, *rv;

  /* Base case, serially scan a chunk. */
  if (n <= chunk_size) {
    scan = (int *) malloc(sizeof(int) * n);

    scan[0] = x[0] + z;
    for (i = 1; i < n; i++) {
      scan[i] = x[i] + scan[i - 1];
    }

    return scan;
  }

  sums = (int *) malloc(sizeof(int) * (n / chunk_size));

  /* Reduce each chunk of the array. */
  cilk_for (i = 0; i < n / chunk_size; i++) {
    sums[i] = reduce(&x[i * chunk_size], chunk_size);
  }

  /* Perform a scan on the sums. */
  sumscan = pscan(sums, n / chunk_size, 0, chunk_size);

  free(sums);

  fsum = (int **) malloc(sizeof(int *) * (n / chunk_size));

  /* Perform a recursive scan on each chunk, using
     the appropriate offset from the sums scan. */
  cilk_for (i = 0; i < n / chunk_size; i++) {
    if (i > 0) {
      fsum[i] = pscan(&x[i * chunk_size], chunk_size, sumscan[i - 1], chunk_size);
    } else {
      fsum[i] = pscan(&x[i * chunk_size], chunk_size, 0, chunk_size);
    }
  }

  free(sumscan);

  rv = (int *) malloc(sizeof(int) * n);

  /* Join the arrays. */
  for (i = 0; i < n / chunk_size; i++) {
    for (j = 0; j < chunk_size; j++) {
      rv[i * chunk_size + j] = fsum[i][j];
    }
  }

  for (i = 0; i < n / chunk_size; i++) {
    free(fsum[i]);
  }

  free(fsum);

  return rv;
}

/* Serial reduction. */
int reduce(int *x, int n)
{
  int i;
  int sum;

  sum = 0;

  for (i = 0; i < n; i++) {
    sum += x[i];
  }

  return sum;
}

而且它有效！好吧，它返回正确的结果。它没有达到我希望的性能。原来的表现是

4 0.000004
8 0.000001
16 0.000002
32 0.000003
64 0.000005
128 0.000011
256 0.000019
512 0.000035
1024 0.000068
2048 0.000130
4096 0.000257
8192 0.000512
16384 0.001129
32768 0.002262
65536 0.004519
131072 0.009065
262144 0.018297
524288 0.037416
1048576 0.078307
2097152 0.157448
4194304 0.313855
8388608 0.625689
16777216 1.251949
33554432 2.589439
67108864 5.084731
134217728 10.402186

对于单线程应用程序，但 Cilk 版本的性能更差，具有以下运行时

4 0.005383
8 0.000011
16 0.000009
32 0.000111
64 0.000055
128 0.000579
256 0.000339
512 0.000544
1024 0.000701
2048 0.001086
4096 0.001265
8192 0.001742
16384 0.002283
32768 0.003891
65536 0.005398
131072 0.009255
262144 0.020736
524288 0.058156
1048576 0.103893
2097152 0.215460
4194304 0.419988
8388608 0.749368
16777216 1.650938
33554432 2.960451
67108864 5.799836
134217728 11.294398

我有一台 24 核的机器，所以我们显然没有看到我们希望在这里得到的加速。我的第一个想法是 Cilk 对递归处理不当，导致超额订阅，但 Cilk 应该专门处理递归。有关如何正确实施此操作的任何提示？我尝试将 cilk_for 添加到 for 循环的底部（释放所有内容）和倒数第二组循环的内部 for 循环（加入数组），但这会进一步降低性能。

任何建议都值得赞赏。

但是，请不要告诉我切换到 here 讨论的 Belloch 的并行扫描算法。我已经在 Cilk 中实现了它，而且效果很好。我想看看我是否可以将它的性能与这个递归解决方案相匹配。

【问题讨论】：

只看代码，你使用的是chunk_size of 2吗？这不会导致缓存中的错误共享吗？
我认真考虑过 -1 试图将编程与时髦主义联系起来......
@twalberg 我遵循一个既定的传统，即把强调功能的新语言称为“时髦”。我认为 Cilk 的基于任务的模型给了它一种功能性的感觉，类似于 Ruby 中的闭包。 Cilk 是 90 年代 MIT 的一个旧项目这一事实进一步增加了时髦主义。
@Leeor 感谢您的建议。我觉得答案可能在于块大小，特别是因为我正在测试多个问题大小，这些问题大小可能都有不同的最佳块大小。如果我学到任何东西，我会做更多的测试并回来。

标签： performance algorithm recursion cilk

【解决方案1】：

我通过为每个问题找到最佳块大小来解决我的性能问题。在该块大小下，（相同）并行版本的性能优于顺序版本。

总之，我的一般方法，尤其是两个块的大小都有一些问题：

我的基准测试方法。在带有调优参数的代码中，使用相同的调优参数值来绘制运行时间与问题大小的关系没有多大意义，因为最佳值取决于问题大小。
块大小为 2 可能存在问题，因为它在最大化并行性的同时，也最大化了递归级别的数量，同样，也最大化了随之而来的开销。
块大小为 2 会阻止矢量化。
正如 Leeor 所建议的，大小为 2 的块也可能导致缓存中的错误共享。

感谢 Leeor 引导我朝着正确的方向前进。

【讨论】：