我有两个位图。我想把它们混合成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,那么将其重新解释为doubles是不好的。
可以使用适当的转换,但使用浮点运算(尤其是双精度(是过分的。使用较小的类型可以使更多的类型适合向量,这样通常会更快,因为一条指令可以处理更多的项目。
此外,122字节的头部不应被放入对齐的阵列中,其在那里的存在立即与实际像素数据的位置不对齐。它可以单独写入输出文件。
例如,一种策略是扩展到16位,使用_mm256_mulhi_epu16
缩放大约80%和大约20%,将它们与_mm256_add_epi16
相加,然后再次缩小到8位。解包到16位,然后打包回8位,对于256位矢量来说,工作有点奇怪,可以把它看作是128位运算的两倍。为了防止过早截断,可以对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_down(a, b) == ~avg_up(~a, ~b)
。这并不会导致大量的补语,因为它们大多相互抵消。如果仍有四舍五入,那么将其留给最后一次操作是有意义的。不过,总是四舍五入可以节省XOR。也许像这样(未测试(
__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);
}