Улучшение рекурсивного преобразования Хадамара
У меня есть следующий код для вычисления преобразования Адамара. Прямо сейчас функция хадамара является узким местом моей программы. Видите ли вы какой-нибудь потенциал для его ускорения? Может по инструкции AVX2? Типичные размеры ввода составляют около 512 или 1024.
С уважением, Том
#include <stdio.h>
void hadamard(double *p, size_t len) {
double tmp = 0.0;
if (len == 2) {
tmp = p[0];
p[0] = tmp + p[1];
p[1] = tmp - p[1];
} else {
hadamard(p, len/2);
hadamard(p+len/2, len/2);
for(int i = 0; i < len/2; i++) {
tmp = p[i];
p[i] = tmp + p[i+len/2];
p[i+len/2] = tmp - p[i+len/2];
}
}
}
int main(int argc, char* argv[]) {
double a[] = {1.0, 2.0, 3.0, 4.0};
hadamard(a, 4);
}
Вы определенно можете избавиться от некоторого кода. Это не улучшит производительность, но сделает его чище. Например, len/2 можно рассчитать только один раз, а не для всей функции.
Думаю, это не по теме ТАК. Насколько я понимаю, вы можете получить обзор кода и предложения по улучшениям на Stack Exchange.
Ссылка: Преобразование Адамара
Это часть кода, в котором мне нужно свернуть векторы и заменить свертку преобразованием Адамара, чтобы уменьшить сложность с O (n ^ 2) до O (N * log (N)). К сожалению, его вызывают так часто, что он все еще остается узким местом. Тем временем я заменил тогда вызовы len / 2, но думаю, что компилятор уже сделал это за меня;)
@ tomseidel1 Вы пересчитываете одни и те же матрицы снова и снова? Если так, то простое вычисление матрицы один раз должно творить чудеса.
Вы можете видеть, что я делаю выше :) У меня рекурсивная реализация, и я вообще не вычисляю никаких матриц, а использую структуру бабочки.
@templatetypedef Явное вычисление матрицы не имело бы смысла, поскольку произведение матрица-вектор - это O(n^2), но преобразование может быть вычислено в O(n*log(n)): en.wikipedia.org/wiki/Fast_Walsh%E2%80%93Hadamard_transform
Вы определенно захотите открыть код большего размера, например, size = 16 (четыре вектора по 4 double в каждом). Вы уверены, что вам нужен double, а не float? В SIMD половина размера элемента означает вдвое больше работы на вектор. С другой стороны, горизонтальный материал может потребовать большего количества перемещений элементов внутри одного вектора.
x265 (видеокодек) имеет макросы NASM x86 для преобразований Адамара в 8- и 16-битные целочисленные данные. См. github.com/videolan/x265/blob/… для макроса HADAMARD, используемого HADAMARD_2D. Это определено в терминах макроса разное. (x264 и x265 используют макросы больше, чем большинство рукописных проектов asm, за этим может быть довольно сложно следить. Возможно, вам поможет разборка скомпилированной версии.) Возможно, вы сможете адаптировать некоторые идеи в горизонтальном направлении. или вертикальная часть этого кода для вашего одномерного варианта использования.
Ответы 1
Вот экспериментальная реализация, основанная на Быстрое преобразование Уолша – Адамара, с оптимизированным первым проходом. Прекрасно компилируется с clang-3.3 или новее и дает хорошие результаты с clang-4.0 или новее (требуется -O2 для правильного встраивания соответствующих функций). Если у вас нет FMA, вам необходимо соединить два нижних элемента hada2_ с -0.0 в hadamard4 (и использовать обычный _mm256_add_pd).
Все проверенные мной версии gcc потребовали бы замены memcpy с помощью встроенных функций ручной загрузки / сохранения, чтобы получить аналогичные результаты.
Кроме того, в качестве упражнения я оставил обработку кейсов len<16. И, возможно, стоит реализовать hadamard32 и, возможно, hadamard64, аналогичный hadamard16, чтобы лучше использовать доступные регистры и уменьшить доступ к памяти. В C++ это можно сделать с помощью рекурсивной реализации шаблона.
#include <immintrin.h> // avx+fma
#include <assert.h> // assert
#include <string.h> // memcpy
inline __m256d hadamard4(__m256d x0123)
{
__m256d x1032 = _mm256_permute_pd(x0123, 5); // [x1, x0, x3, x2]
__m256d hada2 = _mm256_addsub_pd(x1032,x0123); // [x0+x1, x0-x1, x2+x3, x2-x3]
__m256d hada2_= _mm256_permute2f128_pd(hada2, hada2, 1); // [x2+x3, x2-x3, x0+x1, x0-x1]
// if no FMA is available, this can be done with xoring and adding:
__m256d res = _mm256_fmadd_pd(hada2_, _mm256_set_pd(1.0, 1.0, -1.0, -1.0), hada2);
return res;
}
inline void hadamard8(__m256d data[2])
{
__m256d a = hadamard4(data[0]);
__m256d b = hadamard4(data[1]);
data[0] = _mm256_add_pd(a,b);
data[1] = _mm256_sub_pd(a,b);
}
inline void hadamard16(__m256d data[4])
{
hadamard8(data+0);
hadamard8(data+2);
for(int i=0; i<2; ++i)
{
__m256d tmp = data[i];
data[i] = _mm256_add_pd(tmp, data[i+2]);
data[i+2]= _mm256_sub_pd(tmp, data[i+2]);
}
}
void hadamard(double* p, size_t len)
{
assert((len&(len-1))==0); // len must be power of 2
assert(len>=16); // TODO implement fallback for smaller sizes ...
// first pass: hadamard of 16 values each
for(size_t i=0; i<len; i+=16)
{
__m256d data[4];
memcpy(data, p+i, sizeof(data)); // should get optimized to 4x vmovupd
hadamard16(data);
memcpy(p+i, data, sizeof(data)); // should get optimized to 4x vmovupd
}
for(size_t h=32; h<len; h*=2)
{
for(size_t i=0; i<len; i+=2*h)
{
for(size_t j=i; j<i+h; j+=4)
{
__m256d x = _mm256_loadu_pd(p+j);
__m256d y = _mm256_loadu_pd(p+j+h);
_mm256_storeu_pd(p+j, _mm256_add_pd(x,y));
_mm256_storeu_pd(p+j+h, _mm256_sub_pd(x,y));
}
}
}
}
Бенчмаркинг также оставлен в качестве упражнения ;-)
Отказ от ответственности: я не тестировал это. Глядя на это, возможно, я перепутал hada2_ и hada2 в инструкции _mm256_fmadd_pd ...
Большое спасибо @chtz! Это очень помогло. Мне пришлось изменить три вещи: 1. Я изменил порядок аргументов _mm256_set_pd () на -1.0, -1.0, 1, 1, так как крайнее левое значение заканчивается справа, когда я снова конвертирую результат в двойное. . 2. Аргументы для fmadd_pd необходимо поменять местами (как вы правильно указали). 3. Мне потребовался еще один hada2 = _mm256_permute_pd(hada2,5); после аддуба. Ускорение составляет примерно 2-3 раза. Очень здорово! Однако это снова вызывает у меня несколько вопросов, поскольку я хотел бы узнать из этого примера:
1. Я также придумал версию AVX2, которая тоже работает, но не дает никаких ускорений. Для этого я попытался заменить цикл после двух рекурсивных вызовов в приведенных выше функциях на его аналог AVX, то есть выполняя простые _mm256_add и _mm256_sub на соответствующих частях длины 4 переменной p. 2. Как вы узнали, что хороший подход - начать с преобразования 4, а затем подняться вверх? Почему бы не начать с 8 или 16? 3. Следуя идее @Peter Cordes об использовании чисел с плавающей запятой: можем ли мы ожидать здесь еще лучшей производительности для длин 512?
Мне любопытно узнать больше о рабочем процессе, который вы используете, где вычисление матрицы Адамара является узким местом. Не могли бы вы подробнее рассказать о том, что делаете? Можно было бы, скажем, жестко закодировать некоторые из этих матриц в программе или спрятать их в файлах данных, а затем при необходимости использовать mmap. Или, возможно, вы пересчитываете их слишком много раз и можете просто кэшировать то, что создаете.