通过已知索引、收集、分散重新调整的数组缓存友好复制
Cache-friendly copying of an array with readjustment by known index, gather, scatter
假设我们有一个数据数组和另一个带索引的数组。
data = [1, 2, 3, 4, 5, 7]
index = [5, 1, 4, 0, 2, 3]
我们想在 index 的位置从 data 的元素创建一个新数组。结果应该是
[4, 2, 5, 7, 3, 1]
朴素算法适用于 O(N)，但它执行随机内存访问。
你能推荐 CPU 具有相同复杂度的缓存友好算法吗？
PS
在我的特定情况下，数据数组中的所有元素都是整数。
PPS
数组可能包含数百万个元素。
PPPS 我同意 SSE/AVX 或任何其他特定于 x64 的优化
将索引和数据合并到一个数组中。然后使用一些 cache-friendly 排序算法对这些对进行排序（按索引）。然后摆脱索引。 （您可以将 merging/removing 索引与排序算法的 first/last 遍结合起来以稍微优化一下）。
对于 cache-friendly O(N) 排序，使用足够小的基数排序 radix（最多 CPU 缓存中缓存行数的一半）。
这里是 radix-sort-like 算法的 C 实现：
void reorder2(const unsigned size)
{
const unsigned min_bucket = size / kRadix;
const unsigned large_buckets = size % kRadix;
g_counters[0] = 0;
for (unsigned i = 1; i <= large_buckets; ++i)
g_counters[i] = g_counters[i - 1] + min_bucket + 1;
for (unsigned i = large_buckets + 1; i < kRadix; ++i)
g_counters[i] = g_counters[i - 1] + min_bucket;
for (unsigned i = 0; i < size; ++i)
{
const unsigned dst = g_counters[g_index[i] % kRadix]++;
g_sort[dst].index = g_index[i] / kRadix;
g_sort[dst].value = g_input[i];
__builtin_prefetch(&g_sort[dst + 1].value, 1);
}
g_counters[0] = 0;
for (unsigned i = 1; i < (size + kRadix - 1) / kRadix; ++i)
g_counters[i] = g_counters[i - 1] + kRadix;
for (unsigned i = 0; i < size; ++i)
{
const unsigned dst = g_counters[g_sort[i].index]++;
g_output[dst] = g_sort[i].value;
__builtin_prefetch(&g_output[dst + 1], 1);
}
}
它在两个方面不同于基数排序：（1）它不进行计数遍历，因为所有计数器都是预先知道的； (2) 它避免使用 power-of-2 值作为基数。
This C++ code was used for benchmarking（如果你想在 32 位系统上 运行 它，稍微减少 kMaxSize 常量）。
以下是基准测试结果（在具有 6Mb 缓存的 Haswell CPU 上）：
很容易看出即使对于朴素的算法，小数组（低于 ~2 000 000 个元素）也是 cache-friendly。此外，您可能会注意到排序方法在图表的最后一点开始为 cache-unfriendly（L3 缓存中的 size/radix 接近 0.75 缓存行）。在这些限制之间，排序方法比朴素算法更有效。
理论上（如果我们仅将这些算法所需的内存带宽与 64 字节缓存行和 4 字节值进行比较）排序算法应该快 3 倍。实际上，我们的差异要小得多，大约 20%。如果我们为 data 数组使用更小的 16 位值（在这种情况下，排序算法大约快 1.5 倍），这可能会得到改善。
排序方法的另一个问题是当 size/radix 接近某些 power-of-2 时它的 worst-case 行为。这可能会被忽略（因为没有那么多 "bad" 大小）或通过使该算法稍微复杂一些来解决。
如果我们将通道数增加到 3，则所有 3 个通道主要使用 L1 缓存，但内存带宽增加了 60%。我用这段代码得到了实验结果：TL; DR。在（通过实验）确定最佳基数值后，对于大于 4 000 000 的大小（其中 2-pass 算法使用 L3 缓存一次）我得到了更好的结果，但对于较小的数组（其中 2-pass 算法使用 L2缓存两次）。正如预期的那样，16 位数据的性能更好。
结论：性能差异远小于算法复杂度差异，因此天真的方法几乎总是更好；如果性能非常重要并且只使用 2 或 4 字节值，则排序方法更可取。
data = [1, 2, 3, 4, 5, 7]
index = [5, 1, 4, 0, 2, 3]
We want to create a new array from elements of data at position from
index. Result should be
result -> [4, 2, 5, 7, 3, 1]
单线程，一次通过
我认为，对于几百万个元素和 单线程 ，朴素 方法可能是最好的方法。
data和index都是顺序访问（读取）的，这对于CPU缓存来说已经是最优的了。这留下了随机写入，但是写入内存并不像从内存中读取那样对缓存友好。
这只需要一次顺序 传递数据和索引。并且有可能一些（有时很多）写入也已经是 cache-friendly。
为 result 使用多个块 - 多个线程
我们可以为结果分配或使用 cache-friendly 大小的块（块是 result array 中的区域），并多次循环 index 和 data（而他们留在缓存中）。
在每个循环中，我们只将 result 中 适合当前 result-block 的元素写入。这对于写入也是 'cache friendly'，但需要多个循环（循环的数量甚至可能变得相当高 - 即 size of data / size of result-block）。
以上可能是使用 多线程 时的一个选项：data 和 index，即 read-only，将由所有人共享缓存中某个级别的核心（取决于缓存架构）。每个线程中的 result 块将是完全独立的（一个核心永远不必等待另一个核心的结果，或同一区域中的写入）。例如：1000 万个元素 - 每个线程都可以处理一个独立的结果块，比如 500.000 个元素（数字应该是 2 的幂）。
将值组合成对并首先对它们进行排序：这已经比简单的选项花费更多的时间（而且缓存也不友好）。
此外，如果 只有 几百万个元素（整数），也不会产生太大影响。如果我们要谈论数十亿或不适合内存的数据，其他策略可能更可取（例如内存映射结果集，如果它不适合内存）。
如果您的问题处理的数据比此处显示的要多得多，那么最快的方法（可能也是对缓存最友好的方法）就是进行大范围的合并排序操作。
所以你会将输入数据分成合理的块，并让一个单独的线程对每个块进行操作。此操作的结果将是两个与输入非常相似的数组（一个数据和一个目标索引），但是索引将被排序。然后你会有一个最终线程对数据进行合并操作到最终输出数组中。
只要分段选择得当，这应该是一个非常适合缓存的算法。明智的意思是让不同线程使用的数据映射到（您选择的处理器的）不同缓存行，以避免缓存抖动。
如果您有大量数据，这确实是瓶颈，您将需要使用基于块的算法，尽可能从相同的块读取和写入。最多需要遍历 2 次数据以确保新数组完全填充，并且需要适当地设置块大小。伪代码如下。
def populate(index,data,newArray,cache)
blockSize = 1000
for i = 0; i < size(index); i++
//We cached this value earlier
if i in cache
newArray[i] = cache[i]
remove(cache,i)
else
newIndex = index[i]
newValue = data[i]
//Check if this index is in our block
if i%blockSize != newIndex%blockSize
//This index is not in our current block, cache it
cache[newIndex] = newValue
else
//This value is in our current block
newArray[newIndex] = newValue
cache = {}
newArray = []
populate(index,data,newArray,cache)
populate(index,data,newArray,cache)
分析
天真的解决方案按顺序访问索引和数据数组，但按随机顺序访问新数组。由于新数组是随机访问的，因此您最终得到 O(N^2)，其中 N 是数组中的块数。
基于块的解决方案不会从一个块跳到另一个块。它依次读取索引、数据和新数组，以读取和写入相同的块。如果一个索引将在另一个块中，它会被缓存并在它所属的块出现时被检索，或者如果该块已经被传递，它将在第二次传递时被检索。第二遍不会有任何伤害。这是 O(N)。
唯一需要注意的是处理缓存。这里有很多发挥创意的机会，但一般来说，如果大量读取和写入最终位于不同的块上，缓存将会增长，这不是最佳选择。这取决于您的数据构成、这种情况发生的频率以及您的缓存实现。
让我们想象一下，缓存中的所有信息都存在于一个块中，并且适合内存。假设缓存有 y 个元素。天真的方法会随机访问至少 y 次。基于块的方法将在第二遍中获得这些。
我注意到您的索引完全覆盖了该域，但顺序是随机的。
如果您要对索引进行排序，同时对索引数组和数据数组应用相同的操作，则数据数组将成为您想要的结果。
select 有很多排序算法可以满足您的缓存友好标准。但它们的复杂性各不相同。我会考虑快速排序或合并排序。
如果你对这个答案感兴趣，我可以用伪代码详细说明。
我担心这可能不是一个成功的模式。
我们有一段代码表现良好，我们通过删除副本对其进行了优化。
结果是性能很差（由于缓存问题）。我看不出你如何产生一个解决问题的单通道算法。使用 OpenMP，可能允许在多个线程之间共享这将导致的停顿。
我假设重新排序只会以同样的方式发生一次。如果它发生多次，那么事先创建一些更好的策略（通过适当的排序算法）将提高性能
我编写了以下程序来实际测试将目标简单拆分为 N 个块是否有帮助，我的发现是：
a) 即使在最坏的情况下，单线程性能（使用分段写入）也不可能不超过原始策略，并且通常至少差 2 倍
b) 但是，某些细分（可能取决于处理器）和数组大小的性能接近统一，因此表明它实际上会提高 multi-core 性能
这样做的结果是：是的，它比不细分更 "cache-friendly"，但是对于单个线程（并且只有一个重新排序），这对您没有任何帮助。
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
void main(char **ARGS,int ARGC) {
int N=1<<26;
double* source = malloc(N*sizeof(double));
double* target = malloc(N*sizeof(double));
int* idx = malloc(N*sizeof(double));
int i;
for(i=0;i<N;i++) {
source[i]=i;
target[i]=0;
idx[i] = rand() % N ;
};
struct timeval now,then;
gettimeofday(&now,NULL);
for(i=0;i<N;i++) {
target[idx[i]]=source[i];
};
gettimeofday(&then,NULL);
printf("%f\n",(0.0+then.tv_sec*1e6+then.tv_usec-now.tv_sec*1e6-now.tv_usec)/N);
gettimeofday(&now,NULL);
int j;
int targetblocks;
int M = 24;
int targetblocksize = 1<<M;
targetblocks = (N/targetblocksize);
for(i=0;i<N;i++) {
for(j=0;j<targetblocks;j++) {
int k = idx[i];
if ((k>>M) == j) {
target[k]=source[i];
};
};
};
gettimeofday(&then,NULL);
printf("%d,%f\n",targetblocks,(0.0+then.tv_sec*1e6+then.tv_usec-now.tv_sec*1e6-now.tv_usec)/N);
};
假设我们有一个数据数组和另一个带索引的数组。
data = [1, 2, 3, 4, 5, 7]
index = [5, 1, 4, 0, 2, 3]
我们想在 index 的位置从 data 的元素创建一个新数组。结果应该是
[4, 2, 5, 7, 3, 1]
朴素算法适用于 O(N)，但它执行随机内存访问。
你能推荐 CPU 具有相同复杂度的缓存友好算法吗？
PS 在我的特定情况下，数据数组中的所有元素都是整数。
PPS 数组可能包含数百万个元素。
PPPS 我同意 SSE/AVX 或任何其他特定于 x64 的优化
将索引和数据合并到一个数组中。然后使用一些 cache-friendly 排序算法对这些对进行排序（按索引）。然后摆脱索引。 （您可以将 merging/removing 索引与排序算法的 first/last 遍结合起来以稍微优化一下）。
对于 cache-friendly O(N) 排序，使用足够小的基数排序 radix（最多 CPU 缓存中缓存行数的一半）。
这里是 radix-sort-like 算法的 C 实现：
void reorder2(const unsigned size)
{
const unsigned min_bucket = size / kRadix;
const unsigned large_buckets = size % kRadix;
g_counters[0] = 0;
for (unsigned i = 1; i <= large_buckets; ++i)
g_counters[i] = g_counters[i - 1] + min_bucket + 1;
for (unsigned i = large_buckets + 1; i < kRadix; ++i)
g_counters[i] = g_counters[i - 1] + min_bucket;
for (unsigned i = 0; i < size; ++i)
{
const unsigned dst = g_counters[g_index[i] % kRadix]++;
g_sort[dst].index = g_index[i] / kRadix;
g_sort[dst].value = g_input[i];
__builtin_prefetch(&g_sort[dst + 1].value, 1);
}
g_counters[0] = 0;
for (unsigned i = 1; i < (size + kRadix - 1) / kRadix; ++i)
g_counters[i] = g_counters[i - 1] + kRadix;
for (unsigned i = 0; i < size; ++i)
{
const unsigned dst = g_counters[g_sort[i].index]++;
g_output[dst] = g_sort[i].value;
__builtin_prefetch(&g_output[dst + 1], 1);
}
}
它在两个方面不同于基数排序：（1）它不进行计数遍历，因为所有计数器都是预先知道的； (2) 它避免使用 power-of-2 值作为基数。
This C++ code was used for benchmarking（如果你想在 32 位系统上 运行 它，稍微减少 kMaxSize 常量）。
以下是基准测试结果（在具有 6Mb 缓存的 Haswell CPU 上）：
很容易看出即使对于朴素的算法，小数组（低于 ~2 000 000 个元素）也是 cache-friendly。此外，您可能会注意到排序方法在图表的最后一点开始为 cache-unfriendly（L3 缓存中的 size/radix 接近 0.75 缓存行）。在这些限制之间，排序方法比朴素算法更有效。
理论上（如果我们仅将这些算法所需的内存带宽与 64 字节缓存行和 4 字节值进行比较）排序算法应该快 3 倍。实际上，我们的差异要小得多，大约 20%。如果我们为 data 数组使用更小的 16 位值（在这种情况下，排序算法大约快 1.5 倍），这可能会得到改善。
排序方法的另一个问题是当 size/radix 接近某些 power-of-2 时它的 worst-case 行为。这可能会被忽略（因为没有那么多 "bad" 大小）或通过使该算法稍微复杂一些来解决。
如果我们将通道数增加到 3，则所有 3 个通道主要使用 L1 缓存，但内存带宽增加了 60%。我用这段代码得到了实验结果：TL; DR。在（通过实验）确定最佳基数值后，对于大于 4 000 000 的大小（其中 2-pass 算法使用 L3 缓存一次）我得到了更好的结果，但对于较小的数组（其中 2-pass 算法使用 L2缓存两次）。正如预期的那样，16 位数据的性能更好。
结论：性能差异远小于算法复杂度差异，因此天真的方法几乎总是更好；如果性能非常重要并且只使用 2 或 4 字节值，则排序方法更可取。
data = [1, 2, 3, 4, 5, 7]
index = [5, 1, 4, 0, 2, 3]We want to create a new array from elements of data at position from index. Result should be
result -> [4, 2, 5, 7, 3, 1]
单线程，一次通过
我认为，对于几百万个元素和 单线程 ，朴素 方法可能是最好的方法。
data和index都是顺序访问（读取）的，这对于CPU缓存来说已经是最优的了。这留下了随机写入，但是写入内存并不像从内存中读取那样对缓存友好。
这只需要一次顺序 传递数据和索引。并且有可能一些（有时很多）写入也已经是 cache-friendly。
为 result 使用多个块 - 多个线程
我们可以为结果分配或使用 cache-friendly 大小的块（块是 result array 中的区域），并多次循环 index 和 data（而他们留在缓存中）。
在每个循环中，我们只将 result 中 适合当前 result-block 的元素写入。这对于写入也是 'cache friendly'，但需要多个循环（循环的数量甚至可能变得相当高 - 即 size of data / size of result-block）。
以上可能是使用 多线程 时的一个选项：data 和 index，即 read-only，将由所有人共享缓存中某个级别的核心（取决于缓存架构）。每个线程中的 result 块将是完全独立的（一个核心永远不必等待另一个核心的结果，或同一区域中的写入）。例如：1000 万个元素 - 每个线程都可以处理一个独立的结果块，比如 500.000 个元素（数字应该是 2 的幂）。
将值组合成对并首先对它们进行排序：这已经比简单的选项花费更多的时间（而且缓存也不友好）。
此外，如果 只有 几百万个元素（整数），也不会产生太大影响。如果我们要谈论数十亿或不适合内存的数据，其他策略可能更可取（例如内存映射结果集，如果它不适合内存）。
如果您的问题处理的数据比此处显示的要多得多，那么最快的方法（可能也是对缓存最友好的方法）就是进行大范围的合并排序操作。
所以你会将输入数据分成合理的块，并让一个单独的线程对每个块进行操作。此操作的结果将是两个与输入非常相似的数组（一个数据和一个目标索引），但是索引将被排序。然后你会有一个最终线程对数据进行合并操作到最终输出数组中。
只要分段选择得当，这应该是一个非常适合缓存的算法。明智的意思是让不同线程使用的数据映射到（您选择的处理器的）不同缓存行，以避免缓存抖动。
如果您有大量数据，这确实是瓶颈，您将需要使用基于块的算法，尽可能从相同的块读取和写入。最多需要遍历 2 次数据以确保新数组完全填充，并且需要适当地设置块大小。伪代码如下。
def populate(index,data,newArray,cache)
blockSize = 1000
for i = 0; i < size(index); i++
//We cached this value earlier
if i in cache
newArray[i] = cache[i]
remove(cache,i)
else
newIndex = index[i]
newValue = data[i]
//Check if this index is in our block
if i%blockSize != newIndex%blockSize
//This index is not in our current block, cache it
cache[newIndex] = newValue
else
//This value is in our current block
newArray[newIndex] = newValue
cache = {}
newArray = []
populate(index,data,newArray,cache)
populate(index,data,newArray,cache)
分析
天真的解决方案按顺序访问索引和数据数组，但按随机顺序访问新数组。由于新数组是随机访问的，因此您最终得到 O(N^2)，其中 N 是数组中的块数。
基于块的解决方案不会从一个块跳到另一个块。它依次读取索引、数据和新数组，以读取和写入相同的块。如果一个索引将在另一个块中，它会被缓存并在它所属的块出现时被检索，或者如果该块已经被传递，它将在第二次传递时被检索。第二遍不会有任何伤害。这是 O(N)。
唯一需要注意的是处理缓存。这里有很多发挥创意的机会，但一般来说，如果大量读取和写入最终位于不同的块上，缓存将会增长，这不是最佳选择。这取决于您的数据构成、这种情况发生的频率以及您的缓存实现。
让我们想象一下，缓存中的所有信息都存在于一个块中，并且适合内存。假设缓存有 y 个元素。天真的方法会随机访问至少 y 次。基于块的方法将在第二遍中获得这些。
我注意到您的索引完全覆盖了该域，但顺序是随机的。
如果您要对索引进行排序，同时对索引数组和数据数组应用相同的操作，则数据数组将成为您想要的结果。
select 有很多排序算法可以满足您的缓存友好标准。但它们的复杂性各不相同。我会考虑快速排序或合并排序。
如果你对这个答案感兴趣，我可以用伪代码详细说明。
我担心这可能不是一个成功的模式。
我们有一段代码表现良好，我们通过删除副本对其进行了优化。
结果是性能很差（由于缓存问题）。我看不出你如何产生一个解决问题的单通道算法。使用 OpenMP，可能允许在多个线程之间共享这将导致的停顿。
我假设重新排序只会以同样的方式发生一次。如果它发生多次，那么事先创建一些更好的策略（通过适当的排序算法）将提高性能
我编写了以下程序来实际测试将目标简单拆分为 N 个块是否有帮助，我的发现是：
a) 即使在最坏的情况下，单线程性能（使用分段写入）也不可能不超过原始策略，并且通常至少差 2 倍
b) 但是，某些细分（可能取决于处理器）和数组大小的性能接近统一，因此表明它实际上会提高 multi-core 性能
这样做的结果是：是的，它比不细分更 "cache-friendly"，但是对于单个线程（并且只有一个重新排序），这对您没有任何帮助。
#include <stdlib.h>
#include <stdio.h>
#include <sys/time.h>
void main(char **ARGS,int ARGC) {
int N=1<<26;
double* source = malloc(N*sizeof(double));
double* target = malloc(N*sizeof(double));
int* idx = malloc(N*sizeof(double));
int i;
for(i=0;i<N;i++) {
source[i]=i;
target[i]=0;
idx[i] = rand() % N ;
};
struct timeval now,then;
gettimeofday(&now,NULL);
for(i=0;i<N;i++) {
target[idx[i]]=source[i];
};
gettimeofday(&then,NULL);
printf("%f\n",(0.0+then.tv_sec*1e6+then.tv_usec-now.tv_sec*1e6-now.tv_usec)/N);
gettimeofday(&now,NULL);
int j;
int targetblocks;
int M = 24;
int targetblocksize = 1<<M;
targetblocks = (N/targetblocksize);
for(i=0;i<N;i++) {
for(j=0;j<targetblocks;j++) {
int k = idx[i];
if ((k>>M) == j) {
target[k]=source[i];
};
};
};
gettimeofday(&then,NULL);
printf("%d,%f\n",targetblocks,(0.0+then.tv_sec*1e6+then.tv_usec-now.tv_sec*1e6-now.tv_usec)/N);
};