Mat*_* K. 1 c bitmap intrinsics avx avx2
我有 2 个位图。我想以 80:20 的比例混合它们，所以我只需将像素值乘以 0,8 和 0,2。用 C 编写的代码（作为 for 循环）运行良好，但使用 AVX2 指令会导致输出图像错误。
#include <stdio.h>
#include <stdlib.h>
#include <immintrin.h>
#define ARRSIZE 5992826
void main(void){
FILE *bmp = fopen("filepath1", "rb"),
*bmpp = fopen("filepath2", "rb"),
*write = fopen("output", "wb");
unsigned char *a = aligned_alloc(32, ARRSIZE),
*b = aligned_alloc(32, ARRSIZE),
*c = aligned_alloc(32, ARRSIZE);
fread(c, 1, 122, bmp);
rewind(bmp);
fread(a, 1, ARRSIZE, bmp);
fread(b, 1, ARRSIZE, bmpp);
__m256i mm_a, mm_b;
__m256d mm_two = _mm256_set1_pd(2),
mm_eight = _mm256_set1_pd(8);
__m256d mm_c, mm_d,
mm_ten = _mm256_set1_pd(10.0);
int i = 122;
for(; i < ARRSIZE; i+=32){
// c[i] = ((a[i] * 0.8) + (b[i] * 0.2));
mm_a = _mm256_loadu_si256((__m256i *)&(a[i]));
mm_b = _mm256_loadu_si256((__m256i *)&(b[i]));
mm_c = _mm256_div_pd((__m256d)mm_a, mm_ten);
mm_d = _mm256_div_pd((__m256d)mm_b, mm_ten);
mm_a = (__m256i)_mm256_floor_pd(_mm256_mul_pd(mm_c, mm_eight));
mm_b = (__m256i)_mm256_floor_pd(_mm256_mul_pd(mm_d, mm_two));
mm_a = _mm256_add_epi8(mm_a, mm_b);
_mm256_storeu_si256((__m256i *)&(c[i]), mm_a);
}
fwrite(c, 1, ARRSIZE, write);
fclose(bmp);
fclose(bmpp);
fclose(write);
free(a);
free(b);
free(c);
}
您的代码存在一个问题，即向量类型之间的转换不是保值转换，而是重新解释。所以(__m256d)mm_a实际上意味着“将这 32 个字节解释为 4 个双精度数”。这可能没问题，但如果数据是 RGB888 打包的，那么将其重新解释为双打就不好了。
可以使用适当的转换，但为此使用浮点算术（尤其是双精度）就太过分了。使用较小的类型可以使更多的类型适合向量，因此通常速度更快，因为可以使用指令处理更多项目。
此外，122 字节标头不应放入对齐数组中，它的存在会立即使实际像素数据的位置不对齐。它可以单独写入输出文件。
例如，一种策略是扩大到 16 位，使用_mm256_mulhi_epu16缩放大约 80% 和大约 20%，用 相加_mm256_add_epi16，然后再次缩小到 8 位。对于 256 位向量，解包到 16 位然后打包回 8 位的工作方式有点奇怪，可以将其视为并排操作的 2 倍 128 位操作。为了防止过早截断，8位源数据可以通过左移8位进行解包，将数据字节放在相应字的高字节中。这样，乘法高位将创建 16 位中间结果，而不是立即将它们截断为 8 位，这样我们可以在进行加法后进行舍入，这是更合适的（这确实需要额外的移位，并且可以选择添加）。例如这样（未测试）：
const uint16_t scale_a = uint16_t(0x10000 * 0.8);
const uint16_t scale_b = uint16_t(0x10000 - scale_a);
__m256i roundoffset = _mm256_set1_epi16(0x80);
__m256i zero = _mm256_setzero_si256();
for(int i = 0; i < ARRSIZE; i += 32) {
// c[i] = ((a[i] * 0.8) + (b[i] * 0.2));
// c[i] = ((a[i] << 8) * scale_a) + ((b[i] << 8) * scale_b) >> 7;
__m256i raw_a = _mm256_loadu_si256((__m256i *)&(a[i]));
__m256i raw_b = _mm256_loadu_si256((__m256i *)&(b[i]));
__m256i data_al = _mm256_unpacklo_epi8(zero, raw_a);
__m256i data_bl = _mm256_unpacklo_epi8(zero, raw_b);
__m256i data_ah = _mm256_unpackhi_epi8(zero, raw_a);
__m256i data_bh = _mm256_unpackhi_epi8(zero, raw_b);
__m256i scaled_al = _mm256_mulhi_epu16(data_al, _mm256_set1_epi16(scale_a));
__m256i scaled_bl = _mm256_mulhi_epu16(data_bl, _mm256_set1_epi16(scale_b));
__m256i scaled_ah = _mm256_mulhi_epu16(data_ah, _mm256_set1_epi16(scale_a));
__m256i scaled_bh = _mm256_mulhi_epu16(data_bh, _mm256_set1_epi16(scale_b));
__m256i suml = _mm256_add_epi16(scaled_al, scaled_bl);
__m256i sumh = _mm256_add_epi16(scaled_ah, scaled_bh);
__m256i roundedl = _mm256_srli_epi16(_mm256_add_epi16(suml, roundoffset), 8);
__m256i roundedh = _mm256_srli_epi16(_mm256_add_epi16(sumh, roundoffset), 8);
__m256i packed = _mm256_packus_epi16(roundedl, roundedh);
_mm256_storeu_si256((__m256i *)&(c[i]), packed);
}
其中有相当多的洗牌操作，这将吞吐量限制为每 5 个周期迭代一次（在没有其他限制器的情况下），每个周期大约为 1 个像素（作为输出）。
另一种策略是使用_mm256_maddubs_epi16，其混合因子的近似精度较低。它将第二个操作数视为有符号字节并执行有符号饱和，因此这次仅适合比例的 7 位近似值。由于它对 8 位数据进行操作，因此解包较少，但由于需要交错来自两个图像的数据，因此仍然存在一些解包。也许像这样（也没有测试过）：
const uint8_t scale_a = uint8_t(0x80 * 0.8);
const uint8_t scale_b = uint8_t(0x80 - scale_a);
__m256i scale = _mm256_set1_epi16((scale_b << 8) | scale_a);
__m256i roundoffset = _mm256_set1_epi16(0x80);
//__m256i scale = _mm256_set1_epi16();
for(int i = 0; i < ARRSIZE; i += 32) {
// c[i] = ((a[i] * 0.8) + (b[i] * 0.2));
// c[i] = (a[i] * scale_a) + (b[i] * scale_b) >> 7;
__m256i raw_a = _mm256_loadu_si256((__m256i *)&(a[i]));
__m256i raw_b = _mm256_loadu_si256((__m256i *)&(b[i]));
__m256i data_l = _mm256_unpacklo_epi8(raw_a, raw_b);
__m256i data_h = _mm256_unpackhi_epi8(raw_a, raw_b);
__m256i blended_l = _mm256_maddubs_epi16(data_l, scale);
__m256i blended_h = _mm256_maddubs_epi16(data_h, scale);
__m256i roundedl = _mm256_srli_epi16(_mm256_add_epi16(blended_l, roundoffset), 7);
__m256i roundedh = _mm256_srli_epi16(_mm256_add_epi16(blended_h, roundoffset), 7);
__m256i packed = _mm256_packus_epi16(roundedl, roundedh);
_mm256_storeu_si256((__m256i *)&(c[i]), packed);
}
只需 3 次洗牌，吞吐量也许可以达到每 3 个周期 1 次迭代，即每个周期几乎 1.8 个像素。
希望有更好的方法来做到这一点。这两种方法都没有接近最大化乘法，这似乎应该是目标。但我不知道如何到达那里。
另一种策略是使用多轮平均来接近所需的比率，但接近并不那么接近。也许是这样的（未测试）：
for(int i = 0; i < ARRSIZE; i += 32) {
// c[i] = round_somehow((a[i] * 0.8125) + (b[i] * 0.1875));
__m256i raw_a = _mm256_loadu_si256((__m256i *)&(a[i]));
__m256i raw_b = _mm256_loadu_si256((__m256i *)&(b[i]));
__m256i mixed_8_8 = _mm256_avg_epu8(raw_a, raw_b);
__m256i mixed_12_4 = _mm256_avg_epu8(raw_a, mixed_8_8);
__m256i mixed_14_2 = _mm256_avg_epu8(raw_a, mixed_12_4);
__m256i mixed_13_3 = _mm256_avg_epu8(mixed_12_4, mixed_14_2);
_mm256_storeu_si256((__m256i *)&(c[i]), mixed_13_3);
}
但_mm256_avg_epu8归根结底，也许堆叠这么多次是不好的。没有“avg round down”指令，但是avg_down(a, b) == ~avg_up(~a, ~b)。这不会导致互补性的巨大混乱，因为它们中的大多数相互抵消。如果仍有舍入，则将其留到最后一次操作是有意义的。不过，总是向下舍入可以节省异或运算。也许是这样的（未测试）
__m256i ones = _mm256_set1_epi8(-1);
for(int i = 0; i < ARRSIZE; i += 32) {
// c[i] = round_somehow((a[i] * 0.8125) + (b[i] * 0.1875));
__m256i raw_a = _mm256_loadu_si256((__m256i *)&(a[i]));
__m256i raw_b = _mm256_loadu_si256((__m256i *)&(b[i]));
__m256i inv_a = _mm256_xor_si256(ones, raw_a);
__m256i inv_b = _mm256_xor_si256(ones, raw_b);
__m256i mixed_8_8 = _mm256_avg_epu8(inv_a, inv_b);
__m256i mixed_12_4 = _mm256_avg_epu8(inv_a, mixed_8_8);
__m256i mixed_14_2 = _mm256_avg_epu8(inv_a, mixed_12_4);
__m256i mixed_13_3 = _mm256_avg_epu8(_mm256_xor_si256(mixed_12_4, ones),
_mm256_xor_si256(mixed_14_2, ones));
_mm256_storeu_si256((__m256i *)&(c[i]), mixed_13_3);
}