位向量和浮点向量的快速点积
位向量和浮点向量的快速点积
提问于 2013-04-17 12:15:44
我试图在i7上以最有效的方式计算浮点数和位向量之间的点积。实际上，我正在对128或256维向量执行此操作，但为了说明起见，让我编写64维的代码来说明问题：
// a has 64 elements. b is a bitvector of 64 dimensions.
float dot(float *restrict a, uint64_t b) {
float sum = 0;
for(int i=0; b && i<64; i++, b>>=1) {
if (b & 1) sum += a[i];
}
return sum;
}这当然是有效的，但问题是，这是整个程序的时间关键点(消耗了50分钟运行的95%的CPU时间)，所以我迫切需要让它更快。
我的猜测是上面的分支是游戏杀手(防止无序执行，导致错误的分支预测)。我不确定向量指令是否可以在这里使用和帮助。使用带有-std=c99 -march=native -mtune=native -Ofast -funroll-loops的gcc 4.8，我目前得到以下输出
movl $4660, %edx
movl $5, %ecx
xorps %xmm0, %xmm0
.p2align 4,,10
.p2align 3
.L4:
testb $1, %cl
je .L2
addss (%rdx), %xmm0
.L2:
leaq 4(%rdx), %rax
shrq %rcx
testb $1, %cl
je .L8
addss 4(%rdx), %xmm0
.L8:
shrq %rcx
testb $1, %cl
je .L9
addss 4(%rax), %xmm0
.L9:
shrq %rcx
testb $1, %cl
je .L10
addss 8(%rax), %xmm0
.L10:
shrq %rcx
testb $1, %cl
je .L11
addss 12(%rax), %xmm0
.L11:
shrq %rcx
testb $1, %cl
je .L12
addss 16(%rax), %xmm0
.L12:
shrq %rcx
testb $1, %cl
je .L13
addss 20(%rax), %xmm0
.L13:
shrq %rcx
testb $1, %cl
je .L14
addss 24(%rax), %xmm0
.L14:
leaq 28(%rax), %rdx
shrq %rcx
cmpq $4916, %rdx
jne .L4
ret编辑可以对数据进行置换(只要所有参数的置换都是相同的)，顺序并不重要。
我想知道有没有什么东西可以比Chris Dodd的SSE2代码快3倍以上。
新注释: AVX/AVX2代码也是受欢迎的！
编辑2给定一个位向量，我必须将它与128 (或256位，如果是256位)不同的浮点向量相乘(所以一次涉及多个浮点向量也是可以的)。这就是整个过程。任何可以加速整个过程的东西都是受欢迎的！
回答 8
Stack Overflow用户
发布于 2013-04-17 14:50:17
最好的选择是使用SSE ps指令，一次操作4个浮点数。您可以利用float 0.0全为0位的事实，使用andps指令屏蔽掉不需要的元素：
#include <stdint.h>
#include <xmmintrin.h>
union {
uint32_t i[4];
__m128 xmm;
} mask[16] = {
{ 0, 0, 0, 0 },
{ ~0, 0, 0, 0 },
{ 0, ~0, 0, 0 },
{ ~0, ~0, 0, 0 },
{ 0, 0, ~0, 0 },
{ ~0, 0, ~0, 0 },
{ 0, ~0, ~0, 0 },
{ ~0, ~0, ~0, 0 },
{ 0, 0, 0, ~0 },
{ ~0, 0, 0, ~0 },
{ 0, ~0, 0, ~0 },
{ ~0, ~0, 0, ~0 },
{ 0, 0, ~0, ~0 },
{ ~0, 0, ~0, ~0 },
{ 0, ~0, ~0, ~0 },
{ ~0, ~0, ~0, ~0 },
};
float dot(__m128 *a, uint64_t b) {
__m128 sum = { 0.0 };
for (int i = 0; i < 16; i++, b>>=4)
sum += _mm_and_ps(a[i], mask[b&0xf].xmm);
return sum[0] + sum[1] + sum[2] + sum[3];
}如果你希望在掩码中有很多0，那么缩短0可能会更快：
for (int i = 0; b; i++, b >>= 4)
if (b & 0xf)
sum += _mm_and_ps(a[i], mask[b&0xf].xmm);但如果b是随机的，这将是较慢的。
Stack Overflow用户
发布于 2013-10-12 14:43:34
为了扩展Aki Suihkonen的答案，重塑位串对有条件地移动浮点数很有帮助。在下面的解决方案中，使用SSE指令PMOVMASKB和PSHUFB以及指令BLENDVPS的两级比特置换在Core 2 Duo 2.26 the上实现了1.25个元素的处理/周期，这是我的参考C代码的20倍。
编辑:已添加AVX2实现。性能是未知的，因为我不能自己测试它，但预计会有两倍的速度。
这是我的实现和测试平台，解释如下。
SSE4.1 (旧版，适用于AVX2见下文)
代码
/* Includes */
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <smmintrin.h> /* SSE 4.1 */
#include <time.h>
/* Defines */
#define ALIGNTO(n) __attribute__((aligned(n)))
#define USE_PINSRW 1
#define NUM_ITERS 2260000
/**
* Bit mask shuffle.
*
* This version uses a loop to store eight u16 and reloads them as one __m128i.
*/
__m128 bitMaskShuffleStoreAndReload(__m128i mask){
const __m128i perm ALIGNTO(16) = _mm_set_epi8(15, 7, 14, 6, 13, 5, 12, 4,
11, 3, 10, 2, 9, 1, 8, 0);
int i;
uint16_t interMask[8] ALIGNTO(16);
/* Shuffle bitmask */
/* Stage 1 */
for(i=7;i>=0;i--){
interMask[i] = _mm_movemask_epi8(mask);
mask = _mm_slli_epi32(mask, 1);
}
/* Stage 2 */
return _mm_castsi128_ps(
_mm_shuffle_epi8(
_mm_load_si128((const __m128i*)interMask),
perm)
);
}
/**
* Bit mask shuffle.
*
* This version uses the PINSTRW instruction.
*/
__m128 bitMaskShufflePINSRW(__m128i mask){
const __m128i perm ALIGNTO(16) = _mm_set_epi8(15, 7, 14, 6, 13, 5, 12, 4,
11, 3, 10, 2, 9, 1, 8, 0);
__m128i imask ALIGNTO(16);
/* Shuffle bitmask */
/* Stage 1 */
imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 7);
mask = _mm_slli_epi16(mask, 1);
imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 6);
mask = _mm_slli_epi16(mask, 1);
imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 5);
mask = _mm_slli_epi16(mask, 1);
imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 4);
mask = _mm_slli_epi16(mask, 1);
imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 3);
mask = _mm_slli_epi16(mask, 1);
imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 2);
mask = _mm_slli_epi16(mask, 1);
imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 1);
mask = _mm_slli_epi16(mask, 1);
imask = _mm_insert_epi16(imask, _mm_movemask_epi8(mask), 0);
/* Stage 2 */
return _mm_castsi128_ps(
_mm_shuffle_epi8(
imask,
perm)
);
}
/**
* SSE 4.1 implementation.
*/
float dotSSE41(__m128 f[32], unsigned char maskArg[16]){
int i, j, k;
__m128i mask ALIGNTO(16) = _mm_load_si128((const __m128i*)maskArg);
__m128 shufdMask ALIGNTO(16);
__m128 zblended ALIGNTO(16);
__m128 sums ALIGNTO(16) = _mm_setzero_ps();
float sumsf[4] ALIGNTO(16);
/* Shuffle bitmask */
#if USE_PINSRW
shufdMask = bitMaskShufflePINSRW(mask);
#else
shufdMask = bitMaskShuffleStoreAndReload(mask);
#endif
/* Dot product */
for(i=1;i>=0;i--){
for(j=1;j>=0;j--){
for(k=7;k>=0;k--){
zblended = _mm_setzero_ps();
zblended = _mm_blendv_ps(zblended, f[i*16+j+k*2], shufdMask);
sums = _mm_add_ps(sums, zblended);
shufdMask = _mm_castsi128_ps(_mm_slli_epi32(_mm_castps_si128(shufdMask), 1));
}
}
}
/* Final Summation */
_mm_store_ps(sumsf, sums);
return sumsf[0] + sumsf[1] + sumsf[2] + sumsf[3];
}
/**
* Reference C implementation
*/
float dotRefC(float f[128], unsigned char mask[16]){
float sum = 0.0;
int i;
for(i=0;i<128;i++){
sum += ((mask[i>>3]>>(i&7))&1) ? f[i] : 0.0;
}
return sum;
}
/**
* Main
*/
int main(void){
/* Variables */
/* Loop Counter */
int i;
/* Data to process */
float data[128] ALIGNTO(16);
unsigned char mask[16] ALIGNTO(16);
float refCSum, sseSum;
/* Time tracking */
clock_t t1, t2, t3;
double refCTime, sseTime;
/* Initialize mask and float arrays with some random data. */
for(i=0;i<128;i++){
if(i<16)
mask[i]=rand();
data[i] = random();
}
/* RUN TESTS */
t1 = clock();
for(i=0;i<NUM_ITERS;i++){
refCSum = dotRefC(data, mask);
}
t2 = clock();
for(i=0;i<NUM_ITERS;i++){
sseSum = dotSSE41((__m128*)data, mask);
}
t3 = clock();
/* Compute time taken */
refCTime = (double)(t2-t1)/CLOCKS_PER_SEC;
sseTime = (double)(t3-t2)/CLOCKS_PER_SEC;
/* Print out results */
printf("Results:\n"
"RefC: Time: %f Value: %f\n"
"SSE: Time: %f Value: %f\n",
refCTime, refCSum,
sseTime, sseSum);
return 0;
}解释
BLENDVPS使用128位寄存器XMM0的所有四个32位通道中的顶部位来确定是否将其源操作数的相应通道中的值移动到其目标操作数中。当使用MOVAPS加载数据时，会得到4个连续的浮点数:例如，第8、9、10和11个浮点数。当然，它们的选择或取消选择必须由相应的位集控制:例如，位串中的第8、9、10和11位。
问题是，当第一次加载掩码时，这些组的位彼此相邻(在第8、9、10和11个位置)，而实际上它们应该相隔32位；请记住，在某些情况下，它们必须占据每个通道的第31位位置( XMM0寄存器中的第31、63、95和127个位置)。
理想情况下发生的是比特转置，它带来比特0，4，8，12，...在通道0，比特1，5，9，13，...在第一通道，比特2，6，10，14，...在第二道和第3，7，11，15位...在三号车道。因此，先前连续的所有4比特的集合现在间隔32比特，在四个32比特通道中的每一个中一个。然后，它所需要的是一个循环，该循环迭代32次，每次将新的4比特的集合移位到每个通道的顶部比特位置。
不幸的是，x86没有很好的位操作指令，因此，由于缺乏一种干净的方法来进行完美的转置，这里是一个合理的折衷方案。
在掩码中，128位
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127被八条PMOVMASKB和八条PSLLW指令置换，首先到
0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 120
1 9 17 25 33 41 49 57 65 73 81 89 97 105 113 121
2 10 18 26 34 42 50 58 66 74 82 90 98 106 114 122
3 11 19 27 35 43 51 59 67 75 83 91 99 107 115 123
4 12 20 28 36 44 52 60 68 76 84 92 100 108 116 124
5 13 21 29 37 45 53 61 69 77 85 93 101 109 117 125
6 14 22 30 38 46 54 62 70 78 86 94 102 110 118 126
7 15 23 31 39 47 55 63 71 79 87 95 103 111 119 127并通过单个PSHUFB指令
0 8 16 24 32 40 48 56 4 12 20 28 36 44 52 60
64 72 80 88 96 104 112 120 68 76 84 92 100 108 116 124
1 9 17 25 33 41 49 57 5 13 21 29 37 45 53 61
65 73 81 89 97 105 113 121 69 77 85 93 101 109 117 125
2 10 18 26 34 42 50 58 6 14 22 30 38 46 54 62
66 74 82 90 98 106 114 122 70 78 86 94 102 110 118 126
3 11 19 27 35 43 51 59 7 15 23 31 39 47 55 63
67 75 83 91 99 107 115 123 71 79 87 95 103 111 119 127。我们现在迭代四个“运行”，每个“运行”包含八个四位的集合，以32位的间隔分布(正如我们所希望的)，使用这些集合作为BLENDVPS的掩码控制。位混洗固有的笨拙导致了dotSSE41()中看起来很笨拙的三重嵌套循环，但在
clang -Ofast -ftree-vectorize -finline-functions -funroll-loops -msse4.1 -mssse3 dot.c -o dottest无论如何，循环都会展开。内部循环迭代由16次重复组成
blendvps 0x90(%rsi),%xmm1
addps %xmm4,%xmm1
pslld $0x1,%xmm2
movdqa %xmm2,%xmm0
xorps %xmm4,%xmm4。
顺便说一句，我无法确定我的两个bit shuffle版本中哪一个是最快的，所以我在答案中给出了两个实现。
AVX2 (新的，但未经测试)
代码
/* Includes */
#include <stdio.h>
#include <stdint.h>
#include <stdlib.h>
#include <immintrin.h> /* AVX2 */
#include <time.h>
/* Defines */
#define ALIGNTO(n) __attribute__((aligned(n)))
#define NUM_ITERS 2260000
/**
* Bit mask shuffle.
*
* This version uses the PINSTRW instruction.
*/
__m256 bitMaskShufflePINSRW(__m256i mask){
__m256i imask ALIGNTO(32);
/* Shuffle bitmask */
imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 7);
mask = _mm256_slli_epi32(mask, 1);
imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 6);
mask = _mm256_slli_epi32(mask, 1);
imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 5);
mask = _mm256_slli_epi32(mask, 1);
imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 4);
mask = _mm256_slli_epi32(mask, 1);
imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 3);
mask = _mm256_slli_epi32(mask, 1);
imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 2);
mask = _mm256_slli_epi32(mask, 1);
imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 1);
mask = _mm256_slli_epi32(mask, 1);
imask = _mm256_insert_epi32(imask, _mm256_movemask_epi8(mask), 0);
/* Return bitmask */
return _mm256_castsi256_ps(imask);
}
/**
* AVX2 implementation.
*/
float dotAVX2(__m256 f[16], unsigned char maskArg[16]){
int i, j, k;
/* Use _mm_loadu_si128 */
__m256i mask ALIGNTO(32) = _mm256_castsi128_si256(_mm_load_si128((const __m128i*)maskArg));
__m256 shufdMask ALIGNTO(32);
__m256 zblended ALIGNTO(32);
__m256 sums ALIGNTO(32) = _mm256_setzero_ps();
float sumsf[8] ALIGNTO(32);
/* Shuffle bitmask */
shufdMask = bitMaskShufflePINSRW(mask);
shufdMask = _mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(shufdMask), 16));
/* Dot product */
for(i=15;i>=0;i--){
zblended = _mm256_setzero_ps();
/* Replace f[i] with _mm256_loadu_ps((float*)&f[i]) */
zblended = _mm256_blendv_ps(zblended, f[i], shufdMask);
sums = _mm256_add_ps(sums, zblended);
shufdMask = _mm256_castsi256_ps(_mm256_slli_epi32(_mm256_castps_si256(shufdMask), 1));
}
/* Final Summation */
_mm256_store_ps(sumsf, sums);
return sumsf[0] + sumsf[1] + sumsf[2] + sumsf[3] + sumsf[4] + sumsf[5] + sumsf[6] + sumsf[7];
}
/**
* Reference C implementation
*/
float dotRefC(float f[128], unsigned char mask[16]){
float sum = 0.0;
int i;
for(i=0;i<128;i++){
sum += ((mask[i>>3]>>(i&7))&1) ? f[i] : 0.0;
}
return sum;
}
/**
* Main
*/
int main(void){
/* Variables */
/* Loop Counter */
int i;
/* Data to process */
float data[128] ALIGNTO(32);
unsigned char mask[16] ALIGNTO(32);
float refCSum, sseSum;
/* Time tracking */
clock_t t1, t2, t3;
double refCTime, sseTime;
/* Initialize mask and float arrays with some random data. */
for(i=0;i<128;i++){
if(i<16)
mask[i]=rand();
data[i] = random();
}
/* RUN TESTS */
t1 = clock();
for(i=0;i<NUM_ITERS;i++){
refCSum = dotRefC(data, mask);
}
t2 = clock();
for(i=0;i<NUM_ITERS;i++){
sseSum = dotAVX2((__m256*)data, mask);
}
t3 = clock();
/* Compute time taken */
refCTime = (double)(t2-t1)/CLOCKS_PER_SEC;
sseTime = (double)(t3-t2)/CLOCKS_PER_SEC;
/* Print out results */
printf("Results:\n"
"RefC: Time: %f Value: %f\n"
"SSE: Time: %f Value: %f\n",
refCTime, refCSum,
sseTime, sseSum);
return 0;
}解释
使用与SSE4.1相同的概念。不同之处在于，我们现在一次处理8个浮点数，并利用AVX2的256位寄存器和来自ymm寄存器的PMOVMASKB (收集256/8 = 32位)。正因为如此，我们现在有了一个更简单的位掩码混洗和一个更简单的循环。
在掩码中，256位
0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31
32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47
48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63
64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79
80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95
96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111
112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127
128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143
144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159
160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175
176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191
192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207
208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223
224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239
240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255使用8条PMOVMASKB和8条PSLLW指令进行置换，以
0 8 16 24 32 40 48 56 64 72 80 88 96 104 112 120
128 136 144 152 160 168 176 184 192 200 208 216 224 232 240 248
1 9 17 25 33 41 49 57 65 73 81 89 97 105 113 121
129 137 145 153 161 169 177 185 193 201 209 217 225 233 241 249
2 10 18 26 34 42 50 58 66 74 82 90 98 106 114 122
130 138 146 154 162 170 178 186 194 202 210 218 226 234 242 250
3 11 19 27 35 43 51 59 67 75 83 91 99 107 115 123
131 139 147 155 163 171 179 187 195 203 211 219 227 235 243 251
4 12 20 28 36 44 52 60 68 76 84 92 100 108 116 124
132 140 148 156 164 172 180 188 196 204 212 220 228 236 244 252
5 13 21 29 37 45 53 61 69 77 85 93 101 109 117 125
133 141 149 157 165 173 181 189 197 205 213 221 229 237 245 253
6 14 22 30 38 46 54 62 70 78 86 94 102 110 118 126
134 142 150 158 166 174 182 190 198 206 214 222 230 238 246 254
7 15 23 31 39 47 55 63 71 79 87 95 103 111 119 127
135 143 151 159 167 175 183 191 199 207 215 223 231 239 247 255。对于128个元素的带浮点数的位乘积，我们然后在8组16个元素上并行迭代。通过在32个元素的集合上迭代，可以很容易地将此实现扩展到256个元素的DPs。现在只需要一个循环。
具体地说，要将其更改为适用于256个元素的点积，您需要
- 将函数参数的大小加倍。使用
= _mm256_load_si256((const __m256i*)maskArg);. - Delete将掩码加载(
= _mm256_castsi128_si256(_mm_load_si128((const __m128i*)maskArg));)的补偿移位左移16，使其恰好位于bitMaskShufflePINSRW调用的下方。从set 31开始向下运行循环，而不是15：for(i=31;i>=0;i--)
：
因为我的CPU是SSE4.1，所以我既不能测试代码，也不能运行代码，但可以用
clang -Ofast -ftree-vectorize -finline-functions -funroll-loops -mavx2 -msse4.1 -mssse3 dotavx2.c -o dottest干净地编译(无需展开就可以获得更快的代码)，生成以下代码：
(gdb) disas dotAVX2
Dump of assembler code for function dotAVX2:
0x0000000100001730 <+0>: push %rbp
0x0000000100001731 <+1>: mov %rsp,%rbp
0x0000000100001734 <+4>: vmovdqa (%rsi),%xmm0
0x0000000100001738 <+8>: vpslld $0x1,%ymm0,%ymm1
0x000000010000173d <+13>: vpslld $0x1,%ymm1,%ymm2
0x0000000100001742 <+18>: vpmovmskb %ymm2,%eax
0x0000000100001746 <+22>: vpslld $0x1,%ymm2,%ymm2
0x000000010000174b <+27>: vpmovmskb %ymm2,%ecx
0x000000010000174f <+31>: vxorps %ymm3,%ymm3,%ymm3
0x0000000100001753 <+35>: vmovd %ecx,%xmm4
0x0000000100001757 <+39>: vpinsrd $0x1,%eax,%xmm4,%xmm4
0x000000010000175d <+45>: vpmovmskb %ymm1,%eax
0x0000000100001761 <+49>: vpinsrd $0x2,%eax,%xmm4,%xmm1
0x0000000100001767 <+55>: vpslld $0x1,%ymm2,%ymm2
0x000000010000176c <+60>: vpslld $0x1,%ymm2,%ymm4
0x0000000100001771 <+65>: vpslld $0x1,%ymm4,%ymm5
0x0000000100001776 <+70>: vpmovmskb %ymm0,%eax
0x000000010000177a <+74>: vpinsrd $0x3,%eax,%xmm1,%xmm0
0x0000000100001780 <+80>: vpmovmskb %ymm5,%eax
0x0000000100001784 <+84>: vpslld $0x1,%ymm5,%ymm1
0x0000000100001789 <+89>: vpmovmskb %ymm1,%ecx
0x000000010000178d <+93>: vmovd %ecx,%xmm1
0x0000000100001791 <+97>: vpinsrd $0x1,%eax,%xmm1,%xmm1
0x0000000100001797 <+103>: vpmovmskb %ymm4,%eax
0x000000010000179b <+107>: vpinsrd $0x2,%eax,%xmm1,%xmm1
0x00000001000017a1 <+113>: vpmovmskb %ymm2,%eax
0x00000001000017a5 <+117>: vpinsrd $0x3,%eax,%xmm1,%xmm1
0x00000001000017ab <+123>: vinserti128 $0x1,%xmm0,%ymm1,%ymm0
0x00000001000017b1 <+129>: vpslld $0x10,%ymm0,%ymm0
0x00000001000017b6 <+134>: vblendvps %ymm0,0x1e0(%rdi),%ymm3,%ymm1
0x00000001000017c0 <+144>: vpslld $0x1,%ymm0,%ymm0
0x00000001000017c5 <+149>: vblendvps %ymm0,0x1c0(%rdi),%ymm3,%ymm2
0x00000001000017cf <+159>: vaddps %ymm2,%ymm1,%ymm1
0x00000001000017d3 <+163>: vpslld $0x1,%ymm0,%ymm0
0x00000001000017d8 <+168>: vblendvps %ymm0,0x1a0(%rdi),%ymm3,%ymm2
0x00000001000017e2 <+178>: vaddps %ymm2,%ymm1,%ymm1
0x00000001000017e6 <+182>: vpslld $0x1,%ymm0,%ymm0
0x00000001000017eb <+187>: vblendvps %ymm0,0x180(%rdi),%ymm3,%ymm2
0x00000001000017f5 <+197>: vaddps %ymm2,%ymm1,%ymm1
0x00000001000017f9 <+201>: vpslld $0x1,%ymm0,%ymm0
0x00000001000017fe <+206>: vblendvps %ymm0,0x160(%rdi),%ymm3,%ymm2
0x0000000100001808 <+216>: vaddps %ymm2,%ymm1,%ymm1
0x000000010000180c <+220>: vpslld $0x1,%ymm0,%ymm0
0x0000000100001811 <+225>: vblendvps %ymm0,0x140(%rdi),%ymm3,%ymm2
0x000000010000181b <+235>: vaddps %ymm2,%ymm1,%ymm1
0x000000010000181f <+239>: vpslld $0x1,%ymm0,%ymm0
0x0000000100001824 <+244>: vblendvps %ymm0,0x120(%rdi),%ymm3,%ymm2
0x000000010000182e <+254>: vaddps %ymm2,%ymm1,%ymm1
0x0000000100001832 <+258>: vpslld $0x1,%ymm0,%ymm0
0x0000000100001837 <+263>: vblendvps %ymm0,0x100(%rdi),%ymm3,%ymm2
0x0000000100001841 <+273>: vaddps %ymm2,%ymm1,%ymm1
0x0000000100001845 <+277>: vpslld $0x1,%ymm0,%ymm0
0x000000010000184a <+282>: vblendvps %ymm0,0xe0(%rdi),%ymm3,%ymm2
0x0000000100001854 <+292>: vaddps %ymm2,%ymm1,%ymm1
0x0000000100001858 <+296>: vpslld $0x1,%ymm0,%ymm0
0x000000010000185d <+301>: vblendvps %ymm0,0xc0(%rdi),%ymm3,%ymm2
0x0000000100001867 <+311>: vaddps %ymm2,%ymm1,%ymm1
0x000000010000186b <+315>: vpslld $0x1,%ymm0,%ymm0
0x0000000100001870 <+320>: vblendvps %ymm0,0xa0(%rdi),%ymm3,%ymm2
0x000000010000187a <+330>: vaddps %ymm2,%ymm1,%ymm1
0x000000010000187e <+334>: vpslld $0x1,%ymm0,%ymm0
0x0000000100001883 <+339>: vblendvps %ymm0,0x80(%rdi),%ymm3,%ymm2
0x000000010000188d <+349>: vaddps %ymm2,%ymm1,%ymm1
0x0000000100001891 <+353>: vpslld $0x1,%ymm0,%ymm0
0x0000000100001896 <+358>: vblendvps %ymm0,0x60(%rdi),%ymm3,%ymm2
0x000000010000189d <+365>: vaddps %ymm2,%ymm1,%ymm1
0x00000001000018a1 <+369>: vpslld $0x1,%ymm0,%ymm0
0x00000001000018a6 <+374>: vblendvps %ymm0,0x40(%rdi),%ymm3,%ymm2
0x00000001000018ad <+381>: vaddps %ymm2,%ymm1,%ymm1
0x00000001000018b1 <+385>: vpslld $0x1,%ymm0,%ymm0
0x00000001000018b6 <+390>: vblendvps %ymm0,0x20(%rdi),%ymm3,%ymm2
0x00000001000018bd <+397>: vaddps %ymm2,%ymm1,%ymm1
0x00000001000018c1 <+401>: vpslld $0x1,%ymm0,%ymm0
0x00000001000018c6 <+406>: vblendvps %ymm0,(%rdi),%ymm3,%ymm0
0x00000001000018cc <+412>: vaddps %ymm0,%ymm1,%ymm0
0x00000001000018d0 <+416>: vpshufd $0x1,%xmm0,%xmm1
0x00000001000018d5 <+421>: vaddss %xmm1,%xmm0,%xmm1
0x00000001000018d9 <+425>: vmovhlps %xmm0,%xmm0,%xmm2
0x00000001000018dd <+429>: vaddss %xmm1,%xmm2,%xmm1
0x00000001000018e1 <+433>: vpshufd $0x3,%xmm0,%xmm2
0x00000001000018e6 <+438>: vaddss %xmm1,%xmm2,%xmm1
0x00000001000018ea <+442>: vextracti128 $0x1,%ymm0,%xmm0
0x00000001000018f0 <+448>: vaddss %xmm1,%xmm0,%xmm1
0x00000001000018f4 <+452>: vpshufd $0x1,%xmm0,%xmm2
0x00000001000018f9 <+457>: vaddss %xmm1,%xmm2,%xmm1
0x00000001000018fd <+461>: vpshufd $0x3,%xmm0,%xmm2
0x0000000100001902 <+466>: vmovhlps %xmm0,%xmm0,%xmm0
0x0000000100001906 <+470>: vaddss %xmm1,%xmm0,%xmm0
0x000000010000190a <+474>: vaddss %xmm0,%xmm2,%xmm0
0x000000010000190e <+478>: pop %rbp
0x000000010000190f <+479>: vzeroupper
0x0000000100001912 <+482>: retq
End of assembler dump.正如我们所看到的，内核现在是3条指令(vblendvps，vaddps，vpslld)。
Stack Overflow用户
发布于 2013-04-19 21:13:45
如果允许在数据掩码( data float data[128] )中进行略微不同的排列，或者在__m128 mask;的位掩码中进行相应的排列，则可以略微改进克里斯·多德( Chris Dodd )提出的算法。(不包括置换掩码所需的时间，此实现(+开销)大约快25% )。当然，这只是我在评论中提供的想法的快速草稿。
union {
unsigned int i[4];
float f[4];
__m128 xmm;
} mask = { 0xFF00FF00, 0xF0F0F0F0, 0xCCCCCCCC, 0xAAAAAAAA };
float dot2(__m128 *a, __m128 mask);
// 20M times 1.161s
float dotref(__m128 *a, unsigned int *mask) // 20M times 8.174s
{
float z=0.0f;
int i;
for (i=0;i<32;i++) {
if (mask[0] & (0x80000000U >> i)) z+= a[i][0];
if (mask[1] & (0x80000000U >> i)) z+= a[i][1];
if (mask[2] & (0x80000000U >> i)) z+= a[i][2];
if (mask[3] & (0x80000000U >> i)) z+= a[i][3];
}
return z;
}相应的汇编程序实现将是：
dot2:
// warm up stage: fill in initial data and
// set up registers
pxor %xmm1, %xmm1 ;; // clear partial sum1
pxor %xmm2, %xmm2 ;; // clear partial sum2
movaps (%rdi), %xmm3 ;; // register warm up stage1
movaps 16(%rdi), %xmm4 ;; // next 4 values
pxor %xmm5, %xmm5
pxor %xmm6, %xmm6
lea 32(%rdi), %rdi
movl $16, %ecx ;; // process 2x4 items per iteration (total=128)
a: ;; // inner loop -- 2 independent data paths
blendvps %xmm3, %xmm5
pslld $1, %xmm0
movaps (%rdi), %xmm3
blendvps %xmm4, %xmm6
pslld $1, %xmm0
movaps 16(%rdi), %xmm4
addps %xmm5, %xmm1
pxor %xmm5, %xmm5
addps %xmm6, %xmm2
pxor %xmm6, %xmm6
lea 32(%rdi), %rdi
loop a
;; // cool down stage: gather results (xmm0 = xmm1+xmm2)
;; // in beautiful world this stage is interleaved
;; // with the warm up stage of the next block
addps %xmm2, %xmm1
movaps %xmm1, %xmm2
movaps %xmm1, %xmm0
shufps $85, %xmm1, %xmm2
addss %xmm2, %xmm0
movaps %xmm1, %xmm2
unpckhps %xmm1, %xmm2
shufps $255, %xmm1, %xmm1
addss %xmm2, %xmm0
addss %xmm1, %xmm0
ret页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接：
https://stackoverflow.com/questions/16051365
复制相关文章
相似问题