我需要一种快速获取64位整数中所有位为1的位置的方法。例如,给定x = 123703
,我想填充一个数组idx[] = {0, 1, 2, 4, 5, 8, 9, 13, 14, 15, 16}
。我们可以假设我们事先知道位数。这将被称为1012 - 1015倍,因此速度至关重要。到目前为止,我想到的最快的答案是以下的怪物代码,它使用64位整数的每个字节作为索引,并在其中存储每个字节中位为1的数量和位置的表:
int64_t x; // this is the input
unsigned char idx[K]; // this is the array of K bits that are set
unsigned char *dst=idx, *src;
unsigned char zero, one, two, three, four, five; // these hold the 0th-5th bytes
zero = x & 0x0000000000FFUL;
one = (x & 0x00000000FF00UL) >> 8;
two = (x & 0x000000FF0000UL) >> 16;
three = (x & 0x0000FF000000UL) >> 24;
four = (x & 0x00FF00000000UL) >> 32;
five = (x & 0xFF0000000000UL) >> 40;
src=tab0+tabofs[zero ]; COPY(dst, src, n[zero ]);
src=tab1+tabofs[one ]; COPY(dst, src, n[one ]);
src=tab2+tabofs[two ]; COPY(dst, src, n[two ]);
src=tab3+tabofs[three]; COPY(dst, src, n[three]);
src=tab4+tabofs[four ]; COPY(dst, src, n[four ]);
src=tab5+tabofs[five ]; COPY(dst, src, n[five ]);
COPY
是一个开关语句,用于复制最多 8 个字节,n
是一个字节数组,表示一个字节中设置的位数,tabofs
给出了偏移量,指向 tabX
中持有的 X-th 字节中设置的位的位置。这比在我的 Xeon E5-2609 上使用具有 __builtin_ctz()
的展开循环方法快约 3 倍。(见下文。)我目前正在按字典序遍历给定位数的 x
。
还有更好的方法吗?
编辑:添加了一个示例(我随后已经修复)。完整代码可在此处找到:http://pastebin.com/79X8XL2P。注意:GCC 使用 -O2 似乎会优化掉它,但英特尔编译器(我用来编写它)不会……
此外,让我提供一些额外的背景信息,以回答以下评论。目标是对可能的解释变量 N 中的每个可能子集执行统计测试;目前的特定目标是 N=41,但我可以看到某些项目需要 N 达到 45-50。测试基本上涉及将相应的数据子矩阵分解。伪代码如下:
double doTest(double *data, int64_t model) {
int nidx, idx[];
double submatrix[][];
nidx = getIndices(model, idx); // get the locations of ones in model
// copy data into submatrix
for(int i=0; i<nidx; i++) {
for(int j=0; j<nidx; j++) {
submatrix[i][j] = data[idx[i]][idx[j]];
}
}
factorize(submatrix, nidx);
return the_answer;
}
我已经编写了一个版本,在Intel Phi板上运行,应该能在大约15天内完成N=41的情况,其中约5-10%的时间花费在朴素的getIndices()
函数中,所以优化后的版本可以节省一天或更多时间。我正在为NVidia Kepler实现一个版本,但不幸的是,我的问题(大量小矩阵操作)并不完全适合硬件(大量大矩阵操作)。尽管如此,这篇论文提出了一种解决方案,似乎可以在与我相同大小的矩阵上实现数百GFLOPS / s的速度,通过积极展开循环并在寄存器中执行整个分解,但要注意, 矩阵的维度必须在编译时定义。(Phi版本中的循环展开也应有助于减少开销并提高矢量化性能,因此getIndices()
函数将变得更加重要!)所以现在我认为我的核心代码应该看起来更像:
double *data; // move data to GPU/Phi once into shared memory
template<unsigned int K> double doTestUnrolled(int *idx) {
double submatrix[K][K];
// copy data into submatrix
#pragma unroll
for(int i=0; i<K; i++) {
#pragma unroll
for(int j=0; j<K; j++) {
submatrix[i][j] = data[idx[i]][idx[j]];
}
}
factorizeUnrolled<K>(submatrix);
return the_answer;
}
Phi版本通过`cilk_for`循环解决模型,循环从model=0到2的N次方(或者为了测试而选定的一个子集)。但是,现在需要对每个K=1到41个比特设置的模型号进行按字典顺序迭代,以便为GPU批处理工作并分摊内核启动开销,正如doynax所指出的那样。
编辑2: 现在假期结束了,我使用icc 15版本在我的Xeon E5-2602上进行了一些基准测试。我用来进行基准测试的代码在这里:http://pastebin.com/XvrGQUat。我对恰好设置了K个比特位的整数执行位提取,因此在表格中"Base"列中测量了按字典顺序迭代的一些开销。这些操作将在N=48时重复执行2的30次方次(如有必要)。
"CTZ"是一个循环,使用gcc内置函数`__builtin_ctzll`来获取最低位设置。
for(int i=0; i<K; i++) {
idx[i] = __builtin_ctzll(tmp);
lb = tmp & -tmp; // get lowest bit
tmp ^= lb; // remove lowest bit from tmp
}
Mark是Mark的无分支循环:
for(int i=0; i<K; i++) {
*dst = i;
dst += x & 1;
x >>= 1;
}
- Tab1是我的基于表格的原始代码,其中包含以下拷贝宏:
#define COPY(d, s, n) \
switch(n) { \
case 8: *(d++) = *(s++); \
case 7: *(d++) = *(s++); \
case 6: *(d++) = *(s++); \
case 5: *(d++) = *(s++); \
case 4: *(d++) = *(s++); \
case 3: *(d++) = *(s++); \
case 2: *(d++) = *(s++); \
case 1: *(d++) = *(s++); \
case 0: break; \
}
Tab2是与Tab1相同的代码,但复制宏只是将8个字节作为单个复制移动(借鉴了doynax和Lưu Vĩnh Phúc的思路...但请注意,这并不能确保对齐):
#define COPY2(d, s, n) { *((uint64_t *)d) = *((uint64_t *)s); d+=n; }
以下是结果。我认为我的初步观点,即Tab1比CTZ快3倍,只适用于较大的K(我进行测试时)。马克的循环比我的原始代码更快,但是对于K > 8,在COPY2
宏中消除分支是最好的选择。
K Base CTZ Mark Tab1 Tab2
001 4.97s 6.42s 6.66s 18.23s 12.77s
002 4.95s 8.49s 7.28s 19.50s 12.33s
004 4.95s 9.83s 8.68s 19.74s 11.92s
006 4.95s 16.86s 9.53s 20.48s 11.66s
008 4.95s 19.21s 13.87s 20.77s 11.92s
010 4.95s 21.53s 13.09s 21.02s 11.28s
015 4.95s 32.64s 17.75s 23.30s 10.98s
020 4.99s 42.00s 21.75s 27.15s 10.96s
030 5.00s 100.64s 35.48s 35.84s 11.07s
040 5.01s 131.96s 44.55s 44.51s 11.58s