Прийшов час осідлати справжнього Буцефала🏇🏻Приборкай норовливого коня разом з Newxel🏇🏻Умови на сайті
×Закрыть

Покритикуйте реализацию на SSE2

double _sse2_imabsdiffmeanROI_LC_8UC3(const uint8_t * const src0, const uint8_t * const src1,
											 const size_t sz_r0, const size_t sz_c0, const size_t sz_r1, const size_t sz_c1,
											 const size_t sz_r0b, const size_t sz_r0e, const size_t sz_c0b, const size_t sz_c0e,
											 const size_t sz_r1b, const size_t sz_r1e, const size_t sz_c1b, const size_t sz_c1e)
{
	uint32_t res = 0;

	const size_t sz_r0t = sz_r0e-sz_r0b;
	const size_t sz_r1t = sz_r1e-sz_r1b;
	const size_t sz_c0t = sz_c0e-sz_c0b;
	const size_t sz_c1t = sz_c1e-sz_c1b;
	if ((sz_c1t != sz_c0t) || (sz_r1t != sz_r0t))
		return std::numeric_limits<double>::infinity();
	const size_t sz_l0 = sz_r0*sz_c0;
	const size_t sz_l1 = sz_r1*sz_c1;

	const size_t rem = sz_r0t/16*16;
	for(size_t i=0;i<3; ++i)
	{
		for(size_t c=0; c<sz_c0t; ++c)
		{
			const size_t i0 = ((sz_c0b+c)*sz_r0+sz_r0b)+i*sz_l0;
			const size_t i1 = ((sz_c1b+c)*sz_r1+sz_r1b)+i*sz_l1;
			__m128i vsum = _mm_set1_epi32(0);
			for(size_t r=0; r<rem; r+=16)
			{
				const __m128i v0 = _mm_lddqu_si128((__m128i*)(src0+i0+r));
				const __m128i v1 = _mm_lddqu_si128((__m128i*)(src1+i1+r));
				const __m128i v2 = _mm_sad_epu8(v0, v1);
				vsum = _mm_add_epi32(vsum, v2);
			}
			vsum = _mm_add_epi32(vsum, _mm_srli_si128(vsum, 8));
			vsum = _mm_add_epi32(vsum, _mm_srli_si128(vsum, 4));
			res += _mm_cvtsi128_si32(vsum);
			for(size_t r=rem; r<sz_r0t; r++)
			{
				if (src0[i0+r] > src1[i1+r])
					res += src0[i0+r]-src1[i1+r];
				else
					res += src1[i1+r]-src0[i0+r];
			}
		}
	}
	_mm_sfence();
	return (double)res/(sz_c0t*sz_r0t*3);
}

Это функция вычисляет mean(imabsdiff) двух одинаковых прямоугольников на двух картинках в формате матлаба. Цвета по отдельным слоям разложены, в слое по столбцам матрица расположена.
Получилось 80 mks на прямоугольниках 300×600 на АМД FX-8320.

И да могу сказать, зачем она мне понадобилась. В IPP и подобном есть уже подобное. Но ни в одной либе нет функции, что будет сразу считать норму для 3-х мерной матрицы (или я не нашел) с учетом ROI. Есть отдельно inabsdiff и mean. А это запись результирующего массива в память и его вычитывание — это в 3 раза замедлит эту операцию.
Это к тому, что тут недавно вспоминали OpenCV и что нынче нафиг уже не нужно ручками оптимизировать.

LinkedIn
Допустимые теги: blockquote, a, pre, code, ul, ol, li, b, i, del.
Ctrl + Enter
Допустимые теги: blockquote, a, pre, code, ul, ol, li, b, i, del.
Ctrl + Enter

Виктор, может лучше сразу Майку Горчаку в личку писать?

нет, я тоже послушать хочу

Так напрягают програмерские топики на ДОУ? Выбиваются из типичных? Как вайтивайти?

Просто в каждой теме это дискуссия ваша и Майка)

Так присоединяйся. Я с удовольствием прочитаю твои высказывания по данной теме.

Не надо %), у меня личка тут на уровне мусорного ящика, куда приходят подписки со всех мейл-листов (около 1000 писем в день), там потерять письмо проще простого.

Глупый вопрос: а сам процессор, без помощи других компонентов, умеет какие-то SIMD-команды выполнить? Потому как 300×600 это таки не самая маленькая цифра.

Ты вопрос-то корректно сформулируй, а так непонятно, о чем ты спросил.

Та да, уже глянул. В формальной модели этого нет, вывели в спец.наборы инструкций со всеми вытекающими. То есть делать надо только под SSE2.

Пение, ты потрясаешь меня. Я аж улыбнулся от этого твоего поста, спасибо, насмешил.

Напиши нормальный код с выводом конвейра в SSE2, тогда и посмеёшься.
До тех пор это что угодно, только не highload-оптимизация.

Ну ты и пи...бол. Покажи, как надо, если выше всё неправильно.

ну как же
import java.ceeewl.cpu.extended.instructions.sse2.*

Совсем непонятно, зачем тебе столько параметров, неужто типизированные данные типа struct или typdef вызывают какие-то проблемы? ОЧЕНЬ легко ошибиться при вызове функции, вводя кучу параметров. И читать код тяжело.

Если это функция, что мешает передать все параметры по ссылке вместо const? Вроде как разыменование ссылок операция не самая долгая. Могу быть не прав.

Объясни цель операции
const size_t rem = sz_r0t/16*16
по порядку действий ты сначала на 16 поделишь (с плавающей точкой), потом конвертнёшь в целочисленное (твоё счастье что на 16 делишь, а то б ещё на точности единицу мог потерять), и потом умножаешь снова на 16. Даже если там всё правильно, ОЧЕНЬ сложно потом это узнать, что ты планировал. Если это какой-то финт ушами, ПИШИ КОММЕНТ.

Прошлое замечание в силе: убери нахрен цикл от 0 до 3. Лучше впиши по 4 одинаковых строки кода на всё, где переменная i используется. Так ты подскажешь процессору что операция подлежит конвейризации. А читателю станет понятно что ты творишь. Вплоть до того, что впиши в одну строку 4 команды, так очень легко читается.

Есть ли в SSE2 понятие регистра? То есть можно ли с ходу грузануть основные коэфициенты в регистры и не подгружать каждый раз из оперативы?

Проверка if ((sz_c1t != sz_c0t) || (sz_r1t != sz_r0t)) — зачем она тебе? Это действительно необходимо, есть риск выхода за пределы массивов? Да он у тебя и так кажется есть. Считается серьёзной системной угрозой, позволяющей при помощи данных перехватить управление системой. Займись этой темой отдельно.

Не вижу выгоды в применении SSE2 исключительно к сложению. Сложение — одна из самых быстрых операций. Даже для плавающей точки. Но с плавающей точкой ты можешь снизить точность и уйти от этого 8-байтного монстра double.

Лично я бы нарисовал это процедурой, где все данные, включая возвращаемое значение, передавались по ссылкам (то есть void на возврат). А потом поставил перед этой функцией директиву inline. То есть чтобы её код встраивался компилятором прямо в поток команд. Почему так: это достаточно мизерный кусок кода с точки зрения компьютера. Но достаточно обособленный с точки зрения человека. И компилятор сделает грязную работу за тебя, ты увидишь лишь скорость на выходе.

Если это какой-то финт ушами, ПИШИ КОММЕНТ.

const size_t rem = sz_r0t/16*16;

Этому финту ушами скоро будет лет 80. То дед просто математик и пишет, как умеет, а те, кто проходил курс дискретки пишут так: sz_r0t & ~0xFULL; Хотя компилятору всё равно и в конечном итоге он сделает так как я написал, т.к. он уже давно распознаёт паттерны: /16*16, >>4<<4, & ~0xF. Цель — потерять младшие четыре бита числа.

/int*int мне всегда было понятнее читать, чем игры с битами. А процу всё едино уже давно, да и компилятор с таким справиться. А про этот трюк я еще в школе узнал, когда для EC 1022 писал (или какая-то подобная, уже и не вспомню).

Ответил комментом выше. Ты не прав по производительности. Умножение — достаточно медленная операция, деление — экстремально медленная. Тебя спасает только то, что это одна строчка, а не в цикле.

Его спасает то, что умножение на степень двойки любой нормальный компилятор заменяет на сдвиг (а два сдвига — на маскирование).

Ну извини, сейчас на дворе 2018 год, а не 1985.

А чего извиняться-то? ;)

Про спасает ;).
Это в 1985 было важно, а за более чем 30 лет много чего изменилось. Но Пение похоже был в глубокой заморозке.
Сейчас же и такты не столько важны, сколько работа с кешем. В те года память и проц на одной скорости работали, а сейчас разница в 10-ки раз и стадо кешей разных уровней и конвейеры.

Кстати я видел, что gcc понимает умножение от 1 до 10 таким же образом:
3: shl + add
5: shl + shl + add
6: shl + add + shl
etc.

Да, есть такое. И не только до 10. Вот берём такой код (кусок base36 encoding, но я специально заменил 36 на дефайн):

static const char digits[] = "0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ";
void encode(char *tbuf, unsigned long val) {
  // Buffer shall be reversed after setting
  if (val == 0) {
    *tbuf++ = '0';
    *tbuf++ = '\0';
    return;
  }
  while (val > 0) {
    *tbuf++ = digits[val % DDD];
    val = val / DDD;
  }
  *tbuf++ = '\0';
}

и смотрим, во что скомпилировалось при -D DDD=33:

encode:
        movabsq $1117984489315730401, %rcx ; константа для обр. умножения
        testq   %rsi, %rsi
        jne     .L4
        movb    $48, (%rdi)
        movb    $0, 1(%rdi)
        ret
.L4:
        addq    $1, %rdi
        movq    %rsi, %rax
        mulq    %rcx ; умножили на ceil(2/33 * 2**64)
        shrq    %rdx ; финальное деление на 2 для получения частного
        ; теперь надо получить остаток вычитанием
        movq    %rdx, %rax
        salq    $5, %rax ; val*32
        addq    %rdx, %rax ; val+val*32 == val*33
        subq    %rax, %rsi ; остаток от деления
        movzbl  digits(%rsi), %eax
        movb    %al, -1(%rdi)
        movq    %rdx, %rsi
        testq   %rdx, %rdx
        jne     .L4
        ret

Я проверил начиная от 32 числа вверх — первое, где он решил, что явный imul с константой был выгоднее, было 46, затем 58, 70, 78 (просматривается явный шаблон). До этого старался всё сделать через lea.

В ARM/64, я вижу, используется такая же хохма:

        mov     x5, 33761
        movk    x5, 0x3e0f, lsl 16
        movk    x5, 0xe0f8, lsl 32
        movk    x5, 0xf83, lsl 48
.L4:
        umulh   x2, x1, x5
        lsr     x2, x2, 1
        add     x3, x2, x2, lsl 5 ; <-- тут
        sub     x1, x1, x3

а тут на тех же 46, 58... он всё равно начал колдовать сдвиги и вычитания.

Кто знает подноготную, за этот финт могут и яйки усечь на 4 бита.
Умножение — операция не самая быстрая. Если вы знаете как оно происходит, то хорошо представляете, что в процессоре НЕТ таблицы умножения размерностью 18трлн * 18трлн, то есть оно НЕ МОЖЕТ выполниться за 1 такт. С делением всё ещё плачевнее. Так что этот финт ушами оставьте ЖабаСкриптерам, считающим себя программистами. Я уже молчу о том, что нестабильность такта умножения и деления вызывает определённые тормоза в предсказании конвейра. А ещё — требует регистров двойной точности, либо обслужить прерывания при переполнении — это ещё мерзопакостней с точки зрения оптимизации конвейра.

Сравните с битовой операцией, которая гарантированно выполняется в 1 такт, не имеет внутренних прерываний как таковых. Потому всё что можно делать побитово — делается побитово. И всё что можно делать БЕЗ ветвлений алгоритма — делается без них. Именно для того, чтобы не прерывать ламинарный поток, чтобы не зависеть от возможностей оптимизатора процессора (а они весьма скудные), а в идеале — компилировать под специфический процессор или семейство, просто меняя сам код для разных процов, если мы говорим об экстремальной производительности.

Ты написал какую-то дикую смесь кэпства и откровенно сомнительных утверждений.

Если вы знаете как оно происходит, то хорошо представляете, что в процессоре НЕТ таблицы умножения размерностью 18трлн * 18трлн, то есть оно НЕ МОЖЕТ выполниться за 1 такт.

Реально следует считать, что там, где сложение выполняется за 1 такт, умножение выполняется за 3 такта. Это на современных процессорах, где умножение даже 64*64 делается квадратной AND-сеткой аргументов и затем двоичным деревом сумматоров. Сумматоры сами параллельные и работают параллельно, и всё это вместе занимает... смотрим документ Intel® 64 and IA-32 Architectures Optimization Reference Manual (order number 248966), версия 33, таблица С-17:
где ADC/SBB — 1 такт (Skylake) или 2 такта (Broadwell) — умножение «ровное» (64*64->64, IMUL трёхаргументный) — 3 такта, «косое» (64*64->128, двухаргументное) — 4-5 тактов.
Теперь плавающее: ADDSS — 3-4 такта, MULSS — тоже 3-4 — сюрприз! (Кто знает внутренности — ничего удивительного — основное в этих тактах съедают распаковка, запаковка, jamming shift, округление и т.д.)
ADDSD — 3-4 такта, MULSD — 3-5 тактов.
Разница, считаем, вообще отсутствует.
Деление там же: DIVSS — 11 тактов, DIVSD — 14 тактов (внутренних данных нет, но это явно не SRT, это скорее Goldschmidt; про последнего в AMD есть статья). DIV целое — вот тут, да, цифры страшные (до 85). Но оно никак не относится ни к делению на 16, ни к делению на заранее известную константу в любом нормальном компиляторе, и вообще говорят, что Intel забил на его оптимизацию.

Я уже молчу о том, что нестабильность такта умножения и деления вызывает определённые тормоза в предсказании конвейра.

С этим вопросом справились ещё во времена ранних Core.

А ещё — требует регистров двойной точности

С кого это оно требует?
Плавающее умножение, деление не увеличивает размерность.
Целочисленное умножение x86, ARM имеет варианты равной ширины множителей и произведения.
Деление у x86, да, только «косое» (2N/N -> N,N). Но и то есть оптимизация на типичный случай, когда старшая часть является расширением младшей.
Для сравнения, ARM отказался от «косого» деления — видимо, когда стало понятно, что всё равно будут работать через обратное умножение. У x86 же остались раритетные подходы времён S/360.

Но это всё имеет отношение к теме, если считать, что там будет реальное умножение и деление. Потому что его не будет — будут те самые битовые операции.

Дальнейшее кэпство пропускаю, но возражать им Горчаку никакого смысла не было.

(с плавающей точкой)

я что то забыл? с каких пор при делении двух целых чисел получается с плавающей точкой?

с каких пор при делении двух целых чисел получается с плавающей точкой?

Это в JS, Python3, Pascal и аналогичных местах :)
Хотя в Python3 есть //, а в куче других языков есть div.

((sz_c1t != sz_c0t) || (sz_r1t != sz_r0t))

Уговорил, объясню. Входные данный могут оказаться и с разными ROI.

Есть ли в SSE2 понятие регистра?

Интересный вопрос. Начни отсюда ru.wikipedia.org/wiki/SSE2.

Не вижу выгоды в применении SSE2 исключительно к сложению.

А оно есть, выйгрыш в скорости до 16 раз для uint8_t.

Выгода SSE-инстркций прежде всего в конвейрах. Копай в эту сторону. Твоя операция явно поддаётся как последовательности с несвязанными результатами, так и SIMD.

«до 16 раз» — расскажешь на реальном измерении. Не забывай, что это работа со специальными регистрами, и туда надо переместить, и оттуда достать.

Пение, пиши лучше во флудильных топиках, там у тебя лучше получается. А тут ты выглядишь ну очень никак, как программист.

Если честно, то я вообще не понимаю, чего ты тратишь время на SSE и интринсики:

Сделай просто два цикла, как у тебя сейчас, один выравненный и один цикл для остатка, только не вычисляй abs() руками, а замени на билтин:

for(int r=rem; r<sz_r0t; r++)
{
    res += __builtin_abs(src0[i0+r]-src1[i1+r]);
}

И посмотри на результирующий код, результат тебя приятно удивит — так ты на SSE в жизни не напишешь. Оно даже сделало подобие duff’s device.

root@mikenfs:~# gcc —version
gcc (Ubuntu 5.4.0-6ubuntu1~16.04.9) 5.4.0 20160609

Вот тут, как всегда:
gcc.gnu.org/...​s/gcc/Other-Builtins.html

Там список билтин функций, к каждой можешь добавлять __builtin_ префикс, если хочешь её использовать, иначе по усмотрению твоей libc или компилятора.

double _sse2_imabsdiffmeanROI_LC_8UC3_1(const uint8_t * const src0, const uint8_t * const src1,
											 const size_t sz_r0, const size_t sz_c0, const size_t sz_r1, const size_t sz_c1,
											 const size_t sz_r0b, const size_t sz_r0e, const size_t sz_c0b, const size_t sz_c0e,
											 const size_t sz_r1b, const size_t sz_r1e, const size_t sz_c1b, const size_t sz_c1e)
{
	uint32_t res = 0;

	const size_t sz_r0t = sz_r0e-sz_r0b;
	const size_t sz_r1t = sz_r1e-sz_r1b;
	const size_t sz_c0t = sz_c0e-sz_c0b;
	const size_t sz_c1t = sz_c1e-sz_c1b;
	if ((sz_c1t != sz_c0t) || (sz_r1t != sz_r0t))
		return std::numeric_limits<double>::infinity();
	const size_t sz_l0 = sz_r0*sz_c0;
	const size_t sz_l1 = sz_r1*sz_c1;

	const size_t rem = sz_r0t/16*16;
	for(size_t i=0;i<3; ++i)
	{
		for(size_t c=0; c<sz_c0t; ++c)
		{
			const size_t i0 = ((sz_c0b+c)*sz_r0+sz_r0b)+i*sz_l0;
			const size_t i1 = ((sz_c1b+c)*sz_r1+sz_r1b)+i*sz_l1;
			for(size_t r=0; r<rem; r++)
			{
				res += __builtin_abs(src0[i0+r]-src1[i1+r]);
			}
			for(size_t r=rem; r<sz_r0t; r++)
			{
				res += __builtin_abs(src0[i0+r]-src1[i1+r]);
			}
		}
	}
	return (double)res/(sz_c0t*sz_r0t*3);
}
130 mks против 80 mks вверху GCC 5.5.

Я тебя когда-то рассказывал, что основной момент оптимизации для tree-vectorize — это конечный размер цикла, т.е. тебе надо сделать наоборот, чем ты делаешь сейчас. Делишь на 256 и крутишь два цикла от i=0 ; i «<» size/256; i++, внутри его от 0 до 255, затем обрабатываешь остатки, как у тебя. Цикл от 0 до 255 будет самым оптимизированным компилятором.

чего ты тратишь время на SSE и интринсики

я бы даже больше спросил, зачем матлаб и все остальное
если это продакшен то нужно смотреть как бы фреймы с камеры с минимумом копирований загонять на GPU/специальный IP и там их процесить.
если это прототип то слижком много усилий как для прототипа

А нужен именно CPU. И это не прототип. Но собственно алгоритм придумывать мне сильно проще на матлабе. Но начали напрягать меня расчеты по базе. Вот пока считает, я с этим ковыряюсь. Но дальше часть из этого идет на прод (а там теслы ставить не хотят — таков заказчик сейчас).
У меня сейчас самый тормозной момент — это прочитать быстро jpg в матлаб. Получается пока на чтение одного файла 1920×1080×3 около 40 fps c turbojpeg и чуть больше половины времени из этого занимает копирование в формат матлаба.
А то, что я здесь доставал с SSE2 я на этом уже капитально ускорился.

и что на таргет платформе будет тоже x86 c SSE2?

_mm_sfence();

это зачем?

((sz_c0b+c)*sz_r0+sz_r0b)+i*sz_l0

во первых скобочки лишние во вторых sz_r0b+i*sz_l0 не зависит от c и можно вынести за цикл

Скобочки для меня.
Как бы не очень охота плодить лишнюю сущность (так понятнее, что за переменная и что в ней), а на быстродействии не сказывается.

Пробовал с остатком играться и заменять if — по быстродействию одинаково, а так понятнее (мне кажется).

_mm_sfence();

Это где-то увидел, типа выгружает всё из сопроцессора, если не выгрузилось. Но не уверен, что в ней есть смысл.

типа выгружает всё из сопроцессора

на самом деле оно делает так что все store операции до sfence становятся глобально видимыми раньше чем store операции после, все это важно только если у вас несколько потоков

все это важно только если у вас несколько потоков

На самом деле, когда много потоков — это не важно, в SMP все кеши когерентны между процессорами и ядрами. Это важно только если ты передаёшь данные из одного домена в другой, например, между CPU и GPU.

это не про кеш это про memory ordering, этот fence заставляет CPU делать store операции в определенном порядке

Смысла скорее всего нет, но запросто окажется за этой фенечкой хардовое прерывание, которое само по себе может десятки милисекунд отожрать.

Лучше почитай о работе сопроцессора, что именно из созданного должно подлежать деструкции после создания там. А именно, какие ресурсы являются монопольными. Не похоже чтобы ты хоть один такой использовал. Хотя бы потому, что если бы было так, ты бы при многократном вызове функции словил исключение.

Подписаться на комментарии