Внедрение Lanczos SSE/AVX
Кто-нибудь есть какие-либо советы о том, как реализовать алгоритм пересчёта (увеличения и уменьшения) изображений Lanczos с использованием SSE/AVX (встроенные функции или сборка)?
Я посмотрел на некоторые реализации C, но есть много разветвлений, и я не знаю, как это реализовать разумно, используя SSE/AVX.
Пример - нормализованный кардинальный грех:
// C implementation
if (!x)
return sin(x*M_PI)/(x*M_PI);
else
return 1;
// AVX implementation
PXOR ymm0, ymm0
MOVAPD ymm1, [x] // x - array of double
CMPPD ymm0, ymm1, 0 // if (!x)
// what now?
MOVAPD ymm3, [pi] // pi - array of double = M_PI - is there better way?
PMULPD ymm1, ymm3 // ymm1 = x*pi
(SINPD ymm2, ymm1) // found intrinsic _mm256_sin_pd - Intel math library, intrinsic functions are OK with me
DIVPD ymm2, ymm1 // result in ymm2
Как я должен вернуть 1 для значений х == 0? По этим индексам у меня 11...11 (true) после CMPPD.
Кроме того, я делаю это для 8-битного изображения в оттенках серого, так что один пиксель только (0..255). Какое влияние на качество будет иметь использование float вместо double? И можно ли будет работать с u_int8 все время и вообще не конвертировать в реальные числа (ошибка будет, вероятно, существенной)?
1 ответ
Если вы еще не знаете asm или SSE/AVX, изучение одного за раз может быть проще. Написание векторных алгоритмов со встроенными функциями C/C++ даст вам более переносимую реализацию, чем непосредственное использование asm. (Скомпилируйте для 32 против 64 бит и Windows или всего остального, вместо того, чтобы нуждаться в 2 или 4 различных версиях asm (или #ifdef-эквивалентных макросах в вашем asm).
Просмотр выходных данных компилятора при написании C может быть полезен, чтобы убедиться, что загрузки / хранилища происходят так, как вы ожидаете, и компилятор не делает глупостей с раздутым кодом из-за допущений по псевдонимам / выравниванию (отсутствию) или сохранение / генерация констант.
Отладка векторного кода достаточно сложна (так как существует гораздо больше состояний, которые нужно отслеживать, и вы должны мысленно следить за вещами через случайные последовательности).
Я бы начал с поиска некоторых частей C, которые могли бы векторизоваться, если компилятор уже не выполняет автоматическую векторизацию, и использования встроенных функций в C. Затем, когда это сработало, я мог бы взять выходные данные компилятора и настроить его вручную. в местах, где компилятор не создавал оптимальный код. (См. http://agner.org/optimize/)
Что касается обработки данных изображения с плавающей точкой по сравнению с int, хорошо, если вам удастся использовать 16-битную фиксированную точку, это будет быстрее (если не требуется больше инструкций). Смотрите мой ответ на другой вопрос обработки изображений об использовании чисел с плавающей запятой и с фиксированной запятой.
Единственная математическая инструкция в SSE (кроме базового add/sub/mul/div) sqrt, Trig / log / exp - все функции библиотеки. Обратите внимание, что во внутреннем руководстве Intel поле "инструкция" не заполнено, что означает, что оно соответствует нескольким инструкциям. Только компилятор Intel даже обеспечивает эти составные встроенные функции.
В любом случае, вам нужно найти sin Реализация, чтобы встроить, или сохранить некоторые регистры и сделать вызов функции. В зависимости от ABI (окна или все остальное), некоторые или все регистры xmm могут быть закрыты функциями. Используя определенный sin реализация позволит вам узнать, какие регистры ему нужны, и только пролить их. (Поскольку вы программируете в asm, вы можете создавать функции, которые знают больше друг о друге, чем если бы они просто следовали ABI.)
Если вам нужно рассчитатьsin(x*PI)Вы можете сделать заказ sin функция, которая делает это, избавляя от необходимости предварительного умножения на PI. Так как идеальная реализация sin выбирая какой алгоритм использовать на основе диапазона входных данных, вы, вероятно, не сможете получить векторизованную реализацию, которая точна до последнего бита мантиссы. К счастью, вам, вероятно, это не нужно, так что запустите реализацию SSE sin(x).
Условные выражения в векторах SIMD обрабатываются сравнениями, образующими вектор элементов со всеми нулями или со всеми единицами. Затем вы можете И или ИЛИ те на другие векторы. Это хорошо работает для таких вещей, как добавление, где значение идентичности 0, (x + 0 = x, так что вы можете отфильтровать некоторые элементы из вектора перед добавлением вектора в аккумулятор). Если вам нужно выбрать один из двух исходных элементов на основе вектора 0 / -1, вы можете делать операции И / ИЛИ вместе или выполнять ту же работу быстрее с blendvps (скаляр с переменным смешиванием, а не постоянное смешивание).
Эта идея ломается немного, если вы хотите избежать вычисления медленного деления на ноль, в первую очередь, вместо того, чтобы просто делать вычисления для всего, а затем маскировать / смешивать. Так как вы хотите, чтобы результат вышел на 1 когда x == 0.0Лучше всего установить нулевые элементы x в FLT_MIN * 16 или что-то еще до вычисления любого из sin(x*PI)/(x*PI), Таким образом, вы избегаете деления на ноль, и результат деления получается очень близким к 1. Если вам нужно, чтобы оно получилось ровно 1,0f (и нет значения x что делает sin(x*PI) == x*PI с вашим sin реализации), то вам нужно смешать дважды: один раз в числителе и один раз в знаменателе. (Чтобы установить для них обоих одинаковое ненулевое значение).
vxorps xmm15, xmm15, xmm15 ; if you can spare a reg to hold a zero constant
; inside your loop: xmm0 holds { x3, x2, x1, x0 }.
vcmpeqps xmm1, xmm0, xmm15 ;; mnemonic for vcmpps xmm1, xmm0, xmm15, 0.
;; Different predicates are an immediate operand, not different opcodes
vblendvps xmm2, xmm0, [memory_holding_vector_of_float_min], xmm1 ; Or cache it in a reg if you have one to spare
; blendv takes elements from the 2nd source operand when the selector (xmm1) has a 1-bit in the MSB (sign bit)
; xmm2 = (x==0.0f) ? FLT_MIN : x
; xmm1 holds { sin(x3*pi), sin(x2*pi)... }
Обратите внимание, что cmpps имеет более широкий выбор предикатов в версии AVX с кодированием VEX, чем в версии SSE.