Вычисление SAD для 128 элементов с учетом двух массивов uint8_t

У меня есть два массива uint8_t, каждый из которых имеет 64 элемента. «Лучший» способ, который я придумал, чтобы вычислить SAD для всех из них, состоит в том, чтобы загрузить 4x 16 элементов, поместить их в два регистра m128i, а затем поместить их оба в регистр m256. Это делается для обоих массивов uint8_t, например:

__m128i a1, a2, b1, b2, s1, s2;
__m256i u, v, c;

// 128 bit of data x 2
a1 = _mm_set_epi64(*(__m64*)block1, *((__m64*)(block1 + stride)));
block1 += stride + stride;
a2 = _mm_set_epi64(*(__m64*)block1, *((__m64*)(block1 + stride)));

// the upper 128 bits of the result are undefined
u = _mm256_castsi128_si256(a1);
// Copy a to dst, then insert 128 bits from b into dst at the location specified by imm.
u = _mm256_insertf128_si256(u, a2, 0x1);

b1 = _mm_set_epi64(*(__m64*)block2, *((__m64*)(block2 + stride)));
block2 += stride + stride;
b2 = _mm_set_epi64(*(__m64*)block2, *((__m64*)(block2 + stride)));

// the upper 128 bits of the result are undefined
v = _mm256_castsi128_si256(b1);
// Copy a to dst, then insert 128 bits from b into dst at the location specified by imm.
v = _mm256_insertf128_si256(v, b2, 0x1);

Теперь у меня есть два регистра m256, u и v, и я могу вычислить SAD:

c = _mm256_sad_epu8(u, v);

Однако, вероятно, из-за позднего рабочего дня я не смог придумать лучшего способа получения результата... Вот что у меня есть сейчас:

s1 = _mm256_extractf128_si256(c, 0x0);
s2 = _mm256_extractf128_si256(c, 0x1);

int p, q;
p = _mm_extract_epi32(s1, 0x0);
q = _mm_extract_epi32(s1, 0x2);
*result += p + q;

p = _mm_extract_epi32(s2, 0x0);
q = _mm_extract_epi32(s2, 0x2);
*result += p + q;

result — это int, если неясно.

Это генерирует довольно много инструкций. На мой взгляд, это единственный способ загрузить все нужные мне элементы. Однако, вероятно, это не лучший способ получить результат из регистра m256i c.

Что ты говоришь? Можете ли вы помочь мне сделать это более оптимизированным способом?

В совокупности функция выглядит примерно так:

void foobar(uint8_t *block1, uint8_t *block2, int stride, int *result)
{
  *result = 0;
  int i;
  __m128i a1, a2, b1, b2, s1, s2;
  __m256i u, v, c;

  for (i = 0; i < 2; ++i) {
    // loading of uints
    // calculating SAD, and getting result

    block1 += stride; block2 += stride;
    block1 += stride; block2 += stride;
  }
}

Из-за особенностей организации uints я могу загрузить только восемь за раз, прежде чем мне придется увеличивать адрес на stride. Загрузка, т.е. шестнадцати сразу, даст плохой результат.


person DavidS    schedule 04.09.2014    source источник
comment
Можете ли вы объяснить, почему вы не можете загрузить более 8 за раз? Можете ли вы описать словами или диаграммой ASCII структуру памяти 64 пар 8-битных целых чисел без знака и какие пары не учитываются?   -  person Iwillnotexist Idonotexist    schedule 05.09.2014
comment
Посмотрите на цикл for в последнем примере кода. Обратите внимание, как я должен увеличить адреса памяти block1 и block2 на stride. Элементы, которые я хочу, следуют этому уравнению: block1[stride * y + x], где y = [0, > и x = [0, 7]   -  person DavidS    schedule 05.09.2014
comment
Что такое stride? Я так понимаю не 8? Если бы это было так, ваш код выше можно было бы заменить только одной невыровненной загрузкой __m256. Если это 16, я бы все равно использовал нагрузку __m256. Если это не что-то из этого, тогда я бы даже не стал возиться с тасованием, поскольку каждый цикл, который вы тратите на перетасовку 8-байтовых групп, — это цикл, который вы не тратите на MOVDQU, PSADBW и PADDQ.   -  person Iwillnotexist Idonotexist    schedule 05.09.2014
comment
Чередуются ли block1 и block2?   -  person Iwillnotexist Idonotexist    schedule 05.09.2014
comment
Это часть алгоритма оценки движения, а шаг — это, по сути, ширина видео. Так что значение шага может быть любым.   -  person DavidS    schedule 05.09.2014
comment
Это алгоритм ME! О, тогда вы можете использовать VMPSADBW? Эта инструкция - золото для МЕНЯ.   -  person Iwillnotexist Idonotexist    schedule 05.09.2014
comment
Я посмотрел на это, но это потребует, чтобы я вычислил SAD четырех 8-битных блоков, верно? Не могли бы вы привести краткий пример использования? Благодарю вас!   -  person DavidS    schedule 05.09.2014
comment
Конечно, могу я просто подтвердить, что вы SAD-блок 8x8 8-бит без знака против нескольких соседних блоков 8x8? Кроме того, информационный документ Intel по этой теме: здесь   -  person Iwillnotexist Idonotexist    schedule 05.09.2014
comment
да, это идея. За исключением параметра stride, это правильно. Я загружаю восемь и восемь байтов за раз, затем увеличиваю stride и загружаю снова. Одинаковая структура для block1 и block2   -  person DavidS    schedule 05.09.2014
comment
И с каким районом, если хотите, сравнивается блок 8x8? Это, скажем, ± 4 по горизонтали и вертикали относительно центральной точки? Поскольку MPSADBW выполняет четырехэлементный 8-битный SAD без знака в восьми последовательных пространственных местоположениях; Так что, если вам это не удобно, мы должны отказаться от этой идеи.   -  person Iwillnotexist Idonotexist    schedule 06.09.2014
comment
Извините, если я не понимаю вас, но значение SAD block1 и block2 сравнивается с предыдущим значением SAD. После каждого расчета мы переходим к следующей (x, y) точке. Предполагая, что мы находимся на одной и той же оси y, адрес block2 просто увеличивается на 1 каждый раз. Не знаю, ответил ли я вам вообще .. Я не разбираюсь во всей этой оценке движения, но я, вероятно, смогу узнать, можете ли вы уточнить, что вам нужно ответить? Снова простите.   -  person DavidS    schedule 06.09.2014
comment
Я имею в виду, что когда вы переходите к следующей точке (x, y), сколько приращений x вы делаете, прежде чем увеличивать y?   -  person Iwillnotexist Idonotexist    schedule 06.09.2014
comment
Хех, это зависит, правда.. Если вы посмотрите здесь: pastebin.com/XaCFG5RR, вы увидите, что x начинается с left и продолжается до right. Эти значения... ну, не всегда одинаковы. Зависит от видео.   -  person DavidS    schedule 06.09.2014
comment
Почти готово; Я точно после значения range. Что можно сказать о его ценности для общего случая, если вы натыкаетесь на границу? Это 1? 8? 16?   -  person Iwillnotexist Idonotexist    schedule 06.09.2014
comment
У вас была возможность посмотреть на мой пример использования VMPSADBW?   -  person Iwillnotexist Idonotexist    schedule 08.09.2014


Ответы (3)


Что касается получения суммы абсолютных разностей из двух байтовых массивов, то я бы сделал это с SSE:

__m128i sum1 = _mm_sad_epu8(u,v);
__m128i sum2 = _mm_shuffle_epi32(sum1,2);
__m128i sum3 = _mm_add_epi16(sum1,sum2);
int8_t  sum4 = (int8_t)_mm_cvtsi128_si32(sum3);

Я не могу проверить это на AVX2 прямо сейчас, но вот непроверенный код, который я бы попробовал в первую очередь.

__m256i sum1 = _mm256_sad_epu8(u,v);
__m256i sum2 = _mm256_shuffle_epi32(sum1,2);
__m256i sum3 = _mm256_add_epi16(sum1,sum2);  
__m128i sum4 = _mm_add_epi16(_mm256_castsi256_si128(sum3),
_mm256_extracti128_si256(sum3,1));
int8_t  sum5 = (int8_t)_mm_cvtsi128_si32(sum4);

Я могу проверить это позже.

person Z boson    schedule 05.09.2014

Учитывая приведенные выше комментарии, я набросал рабочий пример использования инструкции VMPSADBW для оценки движения 8x8. Я несколько разочарован тем, что для этого генерирует GCC-4.8.1, но это очень хорошее начало. Он включает в себя два теста для проверки функциональности, а также для демонстрации использования моей новой функции sad_block_8x8_range().

Внутренний цикл вычисляет SAD блока 8x8 по сравнению с 8 перекрывающимися блоками в исходном изображении, используя 8 загрузок, 8 VMPSADBWs, 7 вертикальных добавлений, перемешивание и уменьшение. После маскирования недопустимых SAD с помощью | 0xFFFFU младший SAD и его индекс немедленно предоставляются благодаря PHMINPOSUW, инструкции, которая обеспечивает минимум и индекс самого младшего из восьми 16-битных значений без знака в регистре xmm, и если этот SAD даже ниже текущего лучшего, он сохраняется вместе с указанным индексом.

/* Includes */
#include <stdint.h>
#include <string.h>
#include <stdio.h>
#include <immintrin.h>




/* Typedefs */
typedef uint8_t  u8;
typedef uint16_t u16;
typedef uint32_t u32;
typedef uint64_t u64;






/* Functions */

/**
 * 
 * 
 * @param [in]  orig   A pointer into the image within which to run ME. Points
 *                     to a base offset from which a window of maxDx pixels to
 *                     the right and maxDy pixels down is explored to find the
 *                     lowest SAD.
 * @param [in]  os     The span of the original image.
 * @param [in]  ref    A pointer to the 8x8 reference block being SAD-ed for in
 *                     the original image.
 * @param [in]  rs     The span of the 8x8 reference block.
 * @param [in]  maxDx  The width of the search window for ME. Cannot be 0.
 * @param [in]  maxDy  The height of the search window for ME. Cannot be 0.
 * @param [out] sadOut The lowest SAD found.
 * @param [out] dxOut  The corresponding best vector found, x-coordinate.
 * @param [out] dyOut  The corresponding best vector found, y-coordinate.
 */

void sad_block_8x8_range(const u8* orig,
                         unsigned  os,
                         const u8* ref,
                         unsigned  rs,
                         unsigned  maxDx,
                         unsigned  maxDy,
                         unsigned* sadOut,
                         unsigned* dxOut,
                         unsigned* dyOut){
    __m128  tmp01f, tmp23f, tmp45f, tmp67f;
    __m128i tmp01,  tmp23,  tmp45,  tmp67;
    __m256i r01,    r23,    r45,    r67;
    __m256i o0, o1, o2, o3, o4, o5, o6, o7;
    const u8* refTmp;
    const u8* origTmp;
    int i;

    unsigned tmpDx, dx, dy, sad;
    unsigned minDx = 0, minDy = 0, minSAD = 0xFFFF;
    const static u16 MASKTBLw[] = {
        0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF,
        0x0000, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF,
        0x0000, 0x0000, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF,
        0x0000, 0x0000, 0x0000, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF,
        0x0000, 0x0000, 0x0000, 0x0000, 0xFFFF, 0xFFFF, 0xFFFF, 0xFFFF,
        0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xFFFF, 0xFFFF, 0xFFFF,
        0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xFFFF, 0xFFFF,
        0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0x0000, 0xFFFF,
    };
    const static __m128i* MASKTBL = (const __m128i*)MASKTBLw;

    /* Load the eight rows of 8 bytes of the reference block. */
    refTmp = ref;
    tmp01f = _mm_loadl_pi(tmp01f, (const __m64*)(refTmp));refTmp+=rs;/* tmp_a = [ x x x x x x x x 7 6 5 4 3 2 1 0 ] */
    tmp01f = _mm_loadh_pi(tmp01f, (const __m64*)(refTmp));refTmp+=rs;/* tmp_a = [ f e d c b a 9 8 7 6 5 4 3 2 1 0 ] */
    tmp23f = _mm_loadl_pi(tmp23f, (const __m64*)(refTmp));refTmp+=rs;
    tmp23f = _mm_loadh_pi(tmp23f, (const __m64*)(refTmp));refTmp+=rs;
    tmp45f = _mm_loadl_pi(tmp45f, (const __m64*)(refTmp));refTmp+=rs;
    tmp45f = _mm_loadh_pi(tmp45f, (const __m64*)(refTmp));refTmp+=rs;
    tmp67f = _mm_loadl_pi(tmp67f, (const __m64*)(refTmp));refTmp+=rs;
    tmp67f = _mm_loadh_pi(tmp67f, (const __m64*)(refTmp));
    tmp01  = _mm_castps_si128(tmp01f);/* A cast is needed to integer. */
    tmp23  = _mm_castps_si128(tmp23f);
    tmp45  = _mm_castps_si128(tmp45f);
    tmp67  = _mm_castps_si128(tmp67f);

    /**
     * Combine them into 4 ymm registers each holding two rows in duplicate;
     * One in high half and once in low half.
     */

    r01  = _mm256_inserti128_si256(_mm256_castsi128_si256(tmp01), tmp01, 1);/* r_ab = [ f e d c b a 9 8 7 6 5 4 3 2 1 0 f e d c b a 9 8 7 6 5 4 3 2 1 0 ] */
    r23  = _mm256_inserti128_si256(_mm256_castsi128_si256(tmp23), tmp23, 1);
    r45  = _mm256_inserti128_si256(_mm256_castsi128_si256(tmp45), tmp45, 1);
    r67  = _mm256_inserti128_si256(_mm256_castsi128_si256(tmp67), tmp67, 1);

    /* Iterate over x and y of search space. */
    for(dy=0;dy<maxDy;dy++){
        for(dx=0;dx<maxDx;dx+=8){
            /* Broadcast 16-byte rows to both halves of ymm register */
            origTmp = orig + dy*os + dx;
            o0 = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(origTmp)));origTmp += os;
            o1 = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(origTmp)));origTmp += os;
            o2 = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(origTmp)));origTmp += os;
            o3 = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(origTmp)));origTmp += os;
            o4 = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(origTmp)));origTmp += os;
            o5 = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(origTmp)));origTmp += os;
            o6 = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(origTmp)));origTmp += os;

/**
 * Define to 0 if the image can be guaranteed to always have 8 extra allocated
 * bytes beyond its nominal end.
 */
#define NO_OVERALLOCATION 1
            if(NO_OVERALLOCATION && maxDx-dx < 9){/* i.e., maxDx+7-dx < 16, the load size. */
                /**
                 * Special-case code for last row.
                 *      maxDx+7-dx   is the number of bytes to be loaded.
                 *      maxDx-dx     is the number of valid elements.
                 */
#if 1
                __m128i dealigned = _mm_loadu_si128((const __m128i*)(origTmp+maxDx-dx-9));
                __m128i shufmsk   = _mm_add_epi8(_mm_set_epi8(15, 14, 13, 12, 11, 10, 9, 8, 7, 6, 5, 4, 3, 2, 1, 0),
                                                 _mm_set1_epi8(7));
                shufmsk = _mm_add_epi8(shufmsk, _mm_set1_epi8(maxDx-dx));
                o7 = _mm256_broadcastsi128_si256(_mm_shuffle_epi8(dealigned, shufmsk));
#else
                u8 tmpArr[16] = {0};
                for(i=0;i < maxDx+7-dx;i++){
                    tmpArr[i] = (orig+(dy+7)*os+dx)[i];
                }
                o7 = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(tmpArr)));
#endif
            }else{
                o7 = _mm256_broadcastsi128_si256(_mm_loadu_si128((const __m128i*)(origTmp)));
            }


            /**
             * ACTUAL ACTION.
             * 
             * The upper 128-bit lane calculates the SAD for the right 4 bytes
             * of each row of the reference block, while the lower 128-bit lane
             * does similarly for the left 4 bytes of each row of the reference
             * block.
             * 
             * Once the individual SADs for each 4-byte half of each row are
             * obtained against eight consecutive neighbours, add the sixteen
             * 4-byte row halves to get the SADs for the full 8x8 blocks.
             * 
             * After masking for invalid entries, find the minimum SAD and its
             * index using PHMINPOSUW.
             * 
             * Compare the old to the new SAD and if it is a record-setter, save
             * it.
             */

            /* MPSADBWs */
            __m256i s0 = _mm256_mpsadbw_epu8(o0, r01, 1<<5 | 1<<3 | 0<<2 | 0<<0);
            __m256i s1 = _mm256_mpsadbw_epu8(o1, r01, 1<<5 | 3<<3 | 0<<2 | 2<<0);
            __m256i s2 = _mm256_mpsadbw_epu8(o2, r23, 1<<5 | 1<<3 | 0<<2 | 0<<0);
            __m256i s3 = _mm256_mpsadbw_epu8(o3, r23, 1<<5 | 3<<3 | 0<<2 | 2<<0);
            __m256i s4 = _mm256_mpsadbw_epu8(o4, r45, 1<<5 | 1<<3 | 0<<2 | 0<<0);
            __m256i s5 = _mm256_mpsadbw_epu8(o5, r45, 1<<5 | 3<<3 | 0<<2 | 2<<0);
            __m256i s6 = _mm256_mpsadbw_epu8(o6, r67, 1<<5 | 1<<3 | 0<<2 | 0<<0);
            __m256i s7 = _mm256_mpsadbw_epu8(o7, r67, 1<<5 | 3<<3 | 0<<2 | 2<<0);

            /* Accumulate half-row results together into half-block results */
            s0 = _mm256_add_epi16(s0, s1);
            s0 = _mm256_add_epi16(s0, s2);
            s0 = _mm256_add_epi16(s0, s3);
            s0 = _mm256_add_epi16(s0, s4);
            s0 = _mm256_add_epi16(s0, s5);
            s0 = _mm256_add_epi16(s0, s6);
            s0 = _mm256_add_epi16(s0, s7);

            /* Accumulate half-block results into block results */
            __m128i t0 = _mm256_extracti128_si256(s0, 0);
            __m128i t1 = _mm256_extracti128_si256(s0, 1);
            __m128i t  = _mm_add_epi16(t0, t1);

            /* Find horizontal minimum using PHMINPOSUW */
            __m128i hm = maxDx-dx < 8 ? MASKTBL[maxDx-dx] : _mm_setzero_si128();
            __m128i h  = _mm_minpos_epu16(_mm_or_si128(t, hm));
            sad   =      (u16)_mm_extract_epi16(h, 0);
            tmpDx = dx + (u16)_mm_extract_epi16(h, 1);

            /* Save the result if it is the best so far. */
            if(sad < minSAD){
                minDx  = tmpDx;
                minDy  = dy;
                minSAD = sad;
            }
        }
    }

    sadOut && (*sadOut = minSAD);
    dxOut  && (*dxOut  = minDx);
    dyOut  && (*dyOut  = minDy);
}

/**
 * MAIN.
 * 
 * Runs two testcases.
 */

int main(){
    const u8 ref[] = {
        0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07,
        0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f,
        0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17,
        0x18, 0x19, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f,
        0x20, 0x21, 0x22, 0x23, 0x24, 0x25, 0x26, 0x27,
        0x28, 0x29, 0x2a, 0x2b, 0x2c, 0x2d, 0x2e, 0x2f,
        0x30, 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37,
        0x38, 0x39, 0x3a, 0x3b, 0x3c, 0x3d, 0x3e, 0x3f,
    };
    const u8 img0[] = {
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x18, 0x19, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x20, 0x21, 0x22, 0x23, 0x24, 0x25, 0x26, 0x27, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x28, 0x29, 0x2a, 0x2b, 0x2c, 0x2d, 0x2e, 0x2f, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x30, 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x38, 0x39, 0x3a, 0x3b, 0x3c, 0x3d, 0x3e, 0x3f, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
    };
    const u8 img1[] = {
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x00, 0x01, 0x02, 0x03, 0x04, 0x05, 0x06, 0x07,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x08, 0x09, 0x0a, 0x0b, 0x0c, 0x0d, 0x0e, 0x0f,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x10, 0x11, 0x12, 0x13, 0x14, 0x15, 0x16, 0x17,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x18, 0x19, 0x1a, 0x1b, 0x1c, 0x1d, 0x1e, 0x1f,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x20, 0x21, 0x22, 0x23, 0x24, 0x25, 0x26, 0x27,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x28, 0x29, 0x2a, 0x2b, 0x2c, 0x2d, 0x2e, 0x2f,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x30, 0x31, 0x32, 0x33, 0x34, 0x35, 0x36, 0x37,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0x38, 0x39, 0x3a, 0x3b, 0x3c, 0x3d, 0x3e, 0x40,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF,
        0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF, 0xFF
    };
    unsigned sad, dx, dy;


    sad_block_8x8_range(img0, 16, ref, 8, 9, 9, &sad, &dx, &dy);
    if(sad == 0 && dx == 7 && dy == 3){
        printf("Test 1 PASSED!\n");
    }else{
        printf("Test 1 FAILED! (SAD = %u, MV=(%u, %u))\n", sad, dx, dy);
    }


    sad_block_8x8_range(img1, 16, ref, 8, 9, 9, &sad, &dx, &dy);
    if(sad == 1 && dx == 8 && dy == 4){
        printf("Test 2 PASSED!\n");
    }else{
        printf("Test 2 FAILED! (SAD = %u, MV=(%u, %u))\n", sad, dx, dy);
    }

    return 0;
}
person Iwillnotexist Idonotexist    schedule 06.09.2014

Чтобы уточнить, я предположил, что вы хотите оптимизировать следующий код:

int sad_2x64_normal(uint8_t *ptr0, uint8_t *ptr1)
{
  int  sum = 0;
  int  v0,v1;  
  for(int i = 0; i < 64; ++i) {
    v0 = static_cast<int>(*ptr0++);
    v1 = static_cast<int>(*ptr1++);
    sum += abs(v0-v1);
  }
  return sum;
}

Вопрос помечен AVX2, мое короткое решение

int sad_2x64_avx2(uint8_t *ptr0, uint8_t *ptr1)
{
  register __m256i  r0;
  register __m256i  r1;
  register __m256i  r2;
  register __m256i  r3;

  r0 = _mm256_load_si256(reinterpret_cast<__m256i const *>(ptr0)); // load 32 bytes (aligned)
  r1 = _mm256_load_si256(reinterpret_cast<__m256i const *>(ptr1)); // load 32 bytes (aligned)

  r2 = _mm256_sad_epu8(r0, r1);    // 4 unsigned 64 bit value

  r0 = _mm256_load_si256(reinterpret_cast<__m256i const *>(ptr0+32));
  r1 = _mm256_load_si256(reinterpret_cast<__m256i const *>(ptr1+32));

  r3 = _mm256_sad_epu8(r0, r1);   

  r2 = _mm256_add_epi16(r2, r3);
  r2 = _mm256_shuffle_epi32(r2, 0xE8);    
  r2 = _mm256_hadd_epi32(r2, r2);
  r2 = _mm256_permute4x64_epi64(r2, 0xE8);
  r2 = _mm256_shuffle_epi32(r2, 0xE8);   
  r2 = _mm256_hadd_epi32(r2, r2);    
  return _mm_extract_epi16(_mm256_castsi256_si128(r2), 0);
}

То же, но со всеми комментариями

#include "emmintrin.h"
#include "immintrin.h"

// 
// returns \sum_0^{63} abs(ptr0[i]-ptr1[i])
// assume ptr0 and ptr1 are 32 byte aligned
//
int sad_2x64_avx2(uint8_t *ptr0, uint8_t *ptr1)
{
  register __m256i  r0;
  register __m256i  r1;
  register __m256i  r2;
  register __m256i  r3;

  // 1st 32 bytes

  r0 = _mm256_load_si256(reinterpret_cast<__m256i const *>(ptr0));    // load 32 bytes (aligned)
  r1 = _mm256_load_si256(reinterpret_cast<__m256i const *>(ptr1));    // load 32 bytes (aligned)

  r2 = _mm256_sad_epu8(r0, r1);    // results stored as 4x64

  // 2nd 32 bytes

  r0 = _mm256_load_si256(reinterpret_cast<__m256i const *>(ptr0+32));    // load 32 bytes (aligned)
  r1 = _mm256_load_si256(reinterpret_cast<__m256i const *>(ptr1+32));    // load 32 bytes (aligned)

  r3 = _mm256_sad_epu8(r0, r1);


  // after sad_epu8
  //
  //       255                             127             63      31
  //      |                               |               |       |
  // r2 = | 0 | 0 | 0 | a | 0 | 0 | 0 | b | 0 | 0 | 0 | c | 0 | 0 | 0 | d |
  // r3 = | 0 | 0 | 0 | e | 0 | 0 | 0 | f | 0 | 0 | 0 | g | 0 | 0 | 0 | h | 

  r2 = _mm256_add_epi16(r2, r3); 

  // after add_epi16
  //
  //       255                             127             63      31
  //      |                               |               |       |
  // r2 = | 0 | 0 | 0 | i | 0 | 0 | 0 | j | 0 | 0 | 0 | k | 0 | 0 | 0 | l |

  r2 = _mm256_shuffle_epi32(r2, 0xE8); // binary 11_10_10_00

  // after shuffle
  //
  //       255                             127             63      31
  //      |                               |               |       |
  // r2 = | 0 | 0 | 0 | i | 0 | i | 0 | j | 0 | 0 | 0 | k | 0 | k | 0 | l | 

  r2 = _mm256_hadd_epi32(r2, r2);  

  // after hadd
  //
  //       255                             127             63      31
  //      |                               |               |       |
  // r2 = |     i |   i+j |     i |  i+j  |     k |   k+l |     k |   k+l | 

  r2 = _mm256_permute4x64_epi64(r2, 0xE8); // 11_10_10_00
  r2 = _mm256_shuffle_epi32(r2, 0xE8); // binary 11_10_10_00

  // after permute and shuffle
  //
  //       255                             127             63      31
  //      |                               |               |       |
  // r2 = |     i |   i+j |   i+j |  i+j  |     i |   i+j |   i+j |   k+l | 

  r2 = _mm256_hadd_epi32(r2, r2);

  // after hadd and shuffle
  //
  //       255                             127             63      31
  //      |                               |               |       |
  // r2 = | ....  | ....  | ....  | ....  |  .... | ....  | ....  |i+j+k+l| 


  return _mm_extract_epi16(_mm256_castsi256_si128(r2), 0);
}

Надеюсь, поможет!

На моей машине (Haswell Core i7 4900MQ) с g++ 4.8.2

g++ -march=core-avx2 -Wall -Wextra -std=c++11 -O3 ...

Я наблюдал увеличение скорости x58 между обычной версией и версией avx2!

person user3636086    schedule 24.09.2014