Проблема с кодом, который у вас был, заключается в том, что приведение между векторными типами не является сохраняющим значение преобразованием, а представляет собой реинтерпретацию. Таким образом, (__m256d)mm_a на самом деле означает «взять эти 32 байта и интерпретировать их как 4 двойных». Это может быть хорошо, но если данные упакованы в RGB888, то их повторная интерпретация как double не годится.
Можно использовать правильные преобразования, но использование арифметики с плавающей запятой (особенно двойной точности) для этого является излишним. Использование меньших типов приводит к тому, что большее количество из них помещается в векторе, что обычно происходит быстрее, поскольку с инструкцией можно обрабатывать больше элементов.
Кроме того, 122-байтовый заголовок не следует помещать в выровненные массивы, его присутствие в нем сразу же выравнивает положение фактических данных пикселя. Он может быть записан в выходной файл отдельно.
Например, одна из стратегий для этого заключается в расширении до 16 бит, используйте _mm256_mulhi_epu16 для масштабирования приблизительно на 80% и приблизительно на 20%, добавьте их с помощью _mm256_add_epi16, затем снова сузьте до 8 бит. Распаковка в 16-битную и последующая упаковка в 8-битную версию работает немного странно с 256-битными векторами, представьте, что это вдвое больше, чем 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 шаффлами, возможно, пропускная способность может достигать 1 итерации за 3 цикла, что будет почти 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). Это не приводит к огромному беспорядку дополнений, потому что большинство из них отменяет друг друга. Если все еще есть округление, имеет смысл оставить это для последней операции. Всегда округление сохраняет 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);
}