c - 具有结构数组的快速选择具有非线性运行时间

Question

已更新，行下方的原始问题:

我需要计算中位数，并想使用 O(N) 快速选择算法。然而事实证明，当数组不再是 double 平面数组，而是结构数组(其中一个元素是用于中值计算的元素)时，运行时间不再与 O(N) 成比例。

以下平面阵列版本具有近似线性的运行时间:

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

#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;

double quickselect(unsigned long k, unsigned long n, double *arr)
{
  unsigned long i, ir, j, l, mid;
  double a, temp;

  l=1;
  ir=n-1;
  for (;;) {
    if (ir <= l+1) {
      if (ir == l+1 && arr[ir] < arr[l]) {
        SWAP(arr[l],arr[ir])
      }
      return arr[k];
    } else {
      mid=(l+ir) >> 1;
      SWAP(arr[mid],arr[l+1])
      if (arr[l] > arr[ir]) {
        SWAP(arr[l],arr[ir])
      }
      if (arr[l+1] > arr[ir]) {
        SWAP(arr[l+1],arr[ir])
      }
      if (arr[l] > arr[l+1]) {
        SWAP(arr[l],arr[l+1])
      }
      i=l+1;
      j=ir;
      a=arr[l+1];
      for (;;) {
        do i++; while (arr[i] < a);
        do j--; while (arr[j] > a);
        if (j < i) break;
        SWAP(arr[i],arr[j])
      }
      arr[l+1]=arr[j];
      arr[j]=a;
      if (j >= k) ir=j-1;
      if (j <= k) l=i;
    }
  }
}

int main()
{
  unsigned long i, j, k, l, m;
  unsigned long ntest = 1e2;
  unsigned long N[5] = {1e3, 1e4, 1e5, 1e6, 1e7};
  clock_t start, diff;

  int seed = 215342512; //time(NULL);
  srand(seed);

  double *arr = (double*) malloc(N[4] * sizeof(double));

  for (i=0; i<5; i++)
  {
    start = clock();

    for (j=0; j<ntest; j++)
    {
      for (k=0; k<N[i]; k++)
      {
        arr[k] = (double) rand() / (double) RAND_MAX;
      }

      quickselect(N[i] / 2, N[i], arr);
    }

    diff = clock() - start;

    printf("%lu %.5f\n", N[i], (double) diff / CLOCKS_PER_SEC);
  }
}

给予:

1000 0.00228
10000 0.02014
100000 0.19868
1000000 2.01272
10000000 20.41286

但是以下带有结构的版本具有非线性运行时间:

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

#define SWAP(a,b) temp=(a);(a)=(b);(b)=temp;

typedef struct {
  double x;
  double y;
  double z;
  int id;
} point_t;

point_t* quickselect(unsigned long k, unsigned long n, point_t **arr)
{
  unsigned long i, ir, j, l, mid;
  point_t *a, *temp;

  l=1;
  ir=n-1;
  for (;;) {
    if (ir <= l+1) {
      if (ir == l+1 && arr[ir]->x < arr[l]->x) {
        SWAP(arr[l],arr[ir])
      }
      return arr[k];
    } else {
      mid=(l+ir) >> 1;
      SWAP(arr[mid],arr[l+1])
      if (arr[l]->x > arr[ir]->x) {
        SWAP(arr[l],arr[ir])
      }
      if (arr[l+1]->x > arr[ir]->x) {
        SWAP(arr[l+1],arr[ir])
      }
      if (arr[l]->x > arr[l+1]->x) {
        SWAP(arr[l],arr[l+1])
      }
      i=l+1;
      j=ir;
      a=arr[l+1];
      for (;;) {
        do i++; while (arr[i]->x < a->x);
        do j--; while (arr[j]->x > a->x);
        if (j < i) break;
        SWAP(arr[i],arr[j])
      }
      arr[l+1]=arr[j];
      arr[j]=a;
      if (j >= k) ir=j-1;
      if (j <= k) l=i;
    }
  }
}

int main()
{
  unsigned long i, j, k, l, m;
  unsigned long ntest = 1e2;
  unsigned long N[5] = {1e3, 1e4, 1e5, 1e6, 1e7};
  clock_t start, diff;

  int seed = 215342512; //time(NULL);
  srand(seed);

  point_t **ap, *a;

  ap = (point_t**) malloc(N[4] * sizeof(point_t*));

  if (ap == NULL) printf("Error in ap\n");

  a = (point_t*) malloc(N[4] * sizeof(point_t));

  if (a == NULL) printf("Error in a\n");

  for (i=0; i<N[4]; i++)
  {
    ap[i] = a+i;
  }

  for (i=0; i<5; i++)
  {
    start = clock();

    for (j=0; j<ntest; j++)
    {
      for (k=0; k<N[i]; k++)
      {
        ap[k]->x = (double) rand() / (double) RAND_MAX;
      }

      quickselect(N[i] / 2, N[i], ap);
    }

    diff = clock() - start;

    printf("%lu %.5f\n", N[i], (double) diff / CLOCKS_PER_SEC);
  }
}

给予:

1000 0.00224
10000 0.02587
100000 0.37574
1000000 7.18962
10000000 96.34863

两个版本都是用 gcc -O2 编译的(但 -O0 给出了相同的缩放比例)。

这种缩放变化从何而来，如何解决？

请注意，虽然我可以更改结构，但我不能只更改中位数 y因为我还需要知道与中点对应的其他参数。此外，我需要结果数组的快速选择行为(例如 a.y <= m.y 用于 a 左侧的所有 m 和 b.y > m.y 用于 b 右侧的所有 m)。

---

我需要计算中位数，并想使用 O(N) 快速选择算法。然而事实证明，当数组不再是 double 平面数组，而是结构数组(其中一个元素是用于中值计算的元素)时，运行时间不再与 O(N) 成比例。

我使用以下实现:

#define SWAP(a,b) temp=(a); (a)=(b); (b)=temp;

typedef struct point_t point_t;

struct point_t {
  double y;

//  unsigned long something;
//
//  double *something_else;
//
//  double yet_another thing;
//
//  point_t* again_something;
};

void median(point_t *arr, unsigned long n)
{
  unsigned long k = n / 2;
  unsigned long i, ir, j, l, mid;
  point_t a, temp;

  l=0;
  ir=n-1;
  for (;;)
  {
    if (ir <= l+1)
    {
      if (ir == l+1 && arr[ir].y < arr[l].y)
      {
        SWAP(arr[l], arr[ir])
      }
      return arr + k;
    }
    else
    {
      mid = (l + ir) >> 1;
      SWAP(arr[mid], arr[l+1])
      if (arr[l].y > arr[ir].y)
      {
        SWAP(arr[l], arr[ir])
      }
      if (arr[l+1].y > arr[ir].y)
      {
        SWAP(arr[l+1], arr[ir])
      }
      if (arr[l].y > arr[l+1].y)
      {
        SWAP(arr[l], arr[l+1])
      }
      i = l+1;
      j = ir;
      a = arr[l+1];
      for (;;)
      {
        do i++; while (arr[i].y < a.y);
        do j--; while (arr[j].y > a.y);
        if (j < i) break;
        SWAP(arr[i],  arr[j])
      }
      arr[l+1] = arr[j];
      arr[j] = a;
      if (j >= k) ir = j-1;
      if (j <= k) l = i;
    }
  }
}

与 -O2结构被优化掉了(我认为，至少缩放看起来与普通数组相同)并且缩放是线性的。但是，当取消注释该结构的其他组件时，缩放比例不再是线性的。怎么会这样？如何解决这个问题？

请注意，虽然我可以更改结构，但我不能只更改中位数 y因为我还需要知道与中点对应的其他参数。此外，我需要结果数组的快速选择行为(例如 a.y <= m.y 用于 a 左侧的所有 m 和 b.y > m.y 用于 b 右侧的所有 m)。

Answer 1

我认为内存缓存错误可以解释执行时间的非线性增长。在我的 x86_64 架构 PC (Linux + gcc) 中，sizeof(double) 为 8 而 sizeof(point_t) 为 32，因此适合缓存内存的元素更少。但非线性增长的更大原因是通过代码中的指针数组对 point_t 结构的内存访问将很快高度随机化，因此会发生越来越多的缓存未命中...

我修改了如下代码:

--- test2.c
+++ test3.c
         -80,14 +80,12
     if (a == NULL)
         printf("Error in an");

-    for (i = 0; i < N[4]; i++) {
-        ap[i] = a + i;
-    }
- 
     for (i = 0; i < 5; i++) {
         start = clock();

         for (j = 0; j < ntest; j++) {
            for (k = 0; k < N[i]; k++) {
+                ap[k] = a + k;
                 ap[k]->x = (double) rand() / (double) RAND_MAX;
             }

并且执行时间增长更加线性。

带有 double 数组的原始 quicselect():

    1000 0.00000
   10000 0.04000
  100000 0.22000
 1000000 1.98000
10000000 20.73000

带有 point_t * 数组的原始 quickselect():

    1000 0.01000
   10000 0.02000
  100000 0.71000
 1000000 8.64000
10000000 157.77000

与上面的 point_t * 数组完全相同的 quickselect()，但在调用 quickselect() 之前确保数组中的指针按顺序排列 通过应用上面的补丁:

    1000 0.00000
   10000 0.02000
  100000 0.40000
 1000000 4.71000
10000000 49.80000

请注意，即使修改后的版本在计时循环中进行了额外的排序，它仍然更快。

我正在运行 3.2GHz Pentium(R) 双核 CPU E6700、64 位 Linux、gcc 4.6、优化 -O2。 (我的机器没有闲置，所以我的基准数据有一些波动 - 如果我要进行更严格的基准计算，我也会考虑使用 clock_gettime(CLOCK_PROCESS_CPUTIME_ID, ...) 来提高 Linux 系统的准确性排除内核未调度基准进程运行的时间。)

更新:例如，valgrind(如果您的平台支持)可用于分析缓存命中的影响。我修改了程序以接受两个参数，数组大小(对应于数组 N[] 的元素)和测试计数(对应于 ntest)。没有 valgrind 的执行时间，其中 test2 本质上是问题中列出的未修改程序，test4 是重新排列 ap 的修改版本[]调用quickselect()函数前的数组:

bash$ ./test2 10000000 100
Array of 10000000 elements, 100 times: 154.40000 seconds
bash$ ./test4 10000000 100
Array of 10000000 elements, 100 times: 48.45000 seconds

这是使用 cachegrind 工具运行 valgrind 的结果:

bash$ valgrind --tool=cachegrind ./test2 10000000 100
==23563== Cachegrind, a cache and branch-prediction profiler
==23563== Copyright (C) 2002-2010, and GNU GPL'd, by Nicholas Nethercote et al.
==23563== Using Valgrind-3.6.1-Debian and LibVEX; rerun with -h for copyright info
==23563== Command: ./test2 10000000 100
==23563== 
Array of 10000000 elements, 100 times: 1190.24000 seconds
==23563== 
==23563== I   refs:      80,075,384,594
==23563== I1  misses:             1,091
==23563== LLi misses:             1,087
==23563== I1  miss rate:           0.00%
==23563== LLi miss rate:           0.00%
==23563== 
==23563== D   refs:      36,670,292,139  (25,550,518,762 rd   + 11,119,773,377 wr)
==23563== D1  misses:     4,218,722,190  ( 3,223,975,942 rd   +    994,746,248 wr)
==23563== LLd misses:     4,190,889,241  ( 3,198,934,125 rd   +    991,955,116 wr)
==23563== D1  miss rate:           11.5% (          12.6%     +            8.9%  )
==23563== LLd miss rate:           11.4% (          12.5%     +            8.9%  )
==23563== 
==23563== LL refs:        4,218,723,281  ( 3,223,977,033 rd   +    994,746,248 wr)
==23563== LL misses:      4,190,890,328  ( 3,198,935,212 rd   +    991,955,116 wr)
==23563== LL miss rate:             3.5% (           3.0%     +            8.9%  )

和

bash$ valgrind --tool=cachegrind ./test4 10000000 100
==24436== Cachegrind, a cache and branch-prediction profiler
==24436== Copyright (C) 2002-2010, and GNU GPL'd, by Nicholas Nethercote et al.
==24436== Using Valgrind-3.6.1-Debian and LibVEX; rerun with -h for copyright info
==24436== Command: ./test4 10000000 100
==24436== 
Array of 10000000 elements, 100 times: 862.89000 seconds
==24436== 
==24436== I   refs:      82,985,384,485
==24436== I1  misses:             1,093
==24436== LLi misses:             1,089
==24436== I1  miss rate:           0.00%
==24436== LLi miss rate:           0.00%
==24436== 
==24436== D   refs:      36,640,292,192  (24,530,518,829 rd   + 12,109,773,363 wr)
==24436== D1  misses:     2,814,232,350  ( 2,189,229,679 rd   +    625,002,671 wr)
==24436== LLd misses:     2,796,287,872  ( 2,171,294,250 rd   +    624,993,622 wr)
==24436== D1  miss rate:            7.6% (           8.9%     +            5.1%  )
==24436== LLd miss rate:            7.6% (           8.8%     +            5.1%  )
==24436== 
==24436== LL refs:        2,814,233,443  ( 2,189,230,772 rd   +    625,002,671 wr)
==24436== LL misses:      2,796,288,961  ( 2,171,295,339 rd   +    624,993,622 wr)
==24436== LL miss rate:             2.3% (           2.0%     +            5.1%  )

参见 Valgrind manual如何阅读这些统计数据。一个要点是:“在现代机器上，L1 未命中通常会花费大约 10 个周期，而 LL 未命中可能花费多达 200 个周期”。计算两种情况之间的 LLd(最后一级数据 缓存)错误成本差异(每个错误差异乘以 3.2GHz CPU 每 3.2e9 周期/秒假设的 200 个周期)给出

bash$ echo $(( (4190889241 - 2796287872) * 200 / 3200000000 )) seconds
87 seconds

D1 失误在这里贡献很小(如果 D1 失误成本与 LLd 失误成本无关，则总计 91 秒)；由于我们所有的错误(最值得注意的是这台计算机中 LLd 错误的实际成本)，D1 错误可以忽略不计。

test2 和 test4 的运行时间差异约为 106 秒，与上述 86 秒相当接近。这一切都可以变得更准确，但这似乎已经证明了缓存未命中在测试安排中的影响。

附言valgrind 写了一个日志文件，我可以在其中验证它似乎正确检测了 L2 和 L1 缓存大小和类型。