Пытался тут в очередной раз ускорить расчет векторного произведения для звуковых данных. С подачи Ивана в комментариях к моему посту о том, что нужно выполнять сразу четыре произведения, что должно дать существенный эффект.
Однако, для того, что бы оптимизировать, надо найти основу, с чего начинать. В качестве основы я взял, да и развернул циклы по векторному произведению для того простейшего случая, когда длина вектора равна 4. С включенным оптимизатором.
Естественно, выглядит такой код странно. Переписал его же на цикл. Скорость выполнения снизилась почти в два раза. Чему, честно сказать, я удивился. Ни конвейер, ни суперскалярная архитектура и прочие достижения современного процессоростроения ничего не смогли противопоставить старой, как мир, оптимизации в лоб.
После чего я перешел к оптимизации под SSE. Основная проблема тут состоит в том, что изначально данные расположены не так, как было бы удобно их вычислять с использованием SSE. Но так как мне было любопытно, какой выигрыш даст SSE в том случае, если данные уже расположены правильно, то я сделал вид, что они лежат типа правильно.
Результат оказался забавным. Умножение сразу четырех векторов было почти в три раза быстрее, чем скалярный расчет с использованием цикла, и всего лишь в полтора раза быстрее скалярного варианта с развертыванием циклов.
Вот тут уж современные достижения в скалярных вычислениях продемонстрировали себя полностью. Также складывается впечатление, что для скалярных и векторных операции используется одни и те же блоки арифметических операций с плавающей точкой.
Но засада состоит все же в том, что в реальности данные организованы не так, как нужно для корректных вычислений. Поэтому на пути к векторизации вычислений лежит преграда в виде транспонирования векторов. Это потребует дополнительных операций и пока не понятно, даст ли в результате SSE какой-то выигрыш и если даст, то какой?
Набор команд SSE меня слегка разочаровал. Нет, конечно, очень важно сделать векторные команды. Но не менее важно иметь и команды, которые готовят данные к векторным операциям. В данном случае меня интересовали команды сбора данных из разных участков памяти в один вектор, а также команда умножения вектора на скаляр.
Первая из требуемых мне команд появилась лишь в наборе AVX2. Что касается второй команды, то ее нет вовсе. Что, в принципе понятно, так как умножение на скаляр можно заменить скалярным умножением на вектор, заполненный одним и тем же числом. Однако, такая команда, которая умеет весь вектор заполнять одним числом, также появляется лишь в AVX2.
Беда же состоит в том, что мой процессор, хоть и не самый старый, тем не менее AVX2 не поддерживает. Поэтому вкусить немного радости от его использования не удалось. Пришлось ограничится скалярным вариантом. В этом и состоит разочарование и вопрос: неужели в большинстве прикладных задач, требующих векторных вычислений, данные уже лежать так, как надо?
В результате выигрыш от векторизации есть. И для случая векторного произведения он составляет целых 21.5%!
Вот более детальные результаты:
| Тест | Результат, такты |
| Развернутый скалярный цикл | 642 |
| Скалярный цикл | 1149 |
| SSE с формирование вектора из скаляра | 436 |
| SSE с использованием готового вектора, сформ. из скаляра | 425 |
| SSE+Транспонирование с формирование вектора из скаляра | 504 |
| SSE+Транспонирование с использованием готового вектора, сформ. из скаляра | 630 |
Хотя, возможно, я просто не умею их готовить (вектора в смысле)?
На всякий случай выложу фрагмент кода, с помощью которого получил указанные выше результаты. Код писан для случая wavelet-преобразования, поэтому на самом деле рассчитывается V1*V2 и V1*V3, то есть реально рассчитывается 8 векторных произведений, каждый вектор длиной 4.
Попробую на днях и я посчитать :) Только вроде DP4x4InLineSimple и DP4x4InLineLoop работают по разному: первая для каждого нового векторного произведения смещает индекс у V1, а во второй i не участвует во внутреннем цикле... Или я туплю?
ОтветитьУдалитьНет, ты все правильно мыслишь. Это я накасячил. ) Но, на производительность исправление никакого эффекта не оказало.
Удалить> Ни конвейер, ни суперскалярная архитектура и прочие достижения современного процессоростроения ничего не смогли противопоставить старой, как мир, оптимизации в лоб.
ОтветитьУдалитьПомоему это компилятор должен оптимизировать, а процессор должен быть более менее тупым как пробка :) и решать не разрешимые на программном уровне задачи. Современные компиляторы с этим не плохо справляются, просто Delphi похоже совсем никому не нужна стала :( Переходите на поддерживаемые большими компаниями языки :)
Попробовал переписать DP4x4InLineLoop на C и скормить его clang. Правда я заменил v1[j] на v1[i * 2 + j] во вложенном цикле. Теперь он работает также как InLineSimple. Это правильно?
Правильно ли я понимаю что массив v1 - огромный, а вот v2 и v3 - единичные векторы?
В лоб:
static inline void inline_loop(float *v1, float *v2, float *v3, int len) {
for (int i = 0; i < len; i++) {
float l = 0;
float h = 0;
for (int j = 0; j < 4; j++) {
l += v1[i * 2 + j] * v2[j];
h += v1[i * 2 + j] * v3[j];
}
v1[i * 2] = l;
v1[i * 2 + 1] = h;
}
}
Компилятор развернул внутренний цикл, раскидал все по куче регистров, в перемежку грузил данные, умножал и складывал регистры xmmX, но в них только 1 float лежит.
Дальше я решил подсказать компилятору что одновременное вычисление 8 скалярных произведений должно быть эффективнее. Никакого ассемблера, только немного странный код:
static inline void inline_loop_2(float *v1, float *v2, float *v3, int len) {
for (int i = 0; i < len; i+=4) {
float l[4] = {0};
float h[4] = {0};
for (int j = 0; j < 4; j++) {
for (int k = 0; k < 4; k++) {
l[k] += v1[(i + k) * 2 + j] * v2[j];
h[k] += v1[(i + k) * 2 + j] * v3[j];
}
}
for (int j = 0; j < 4; j++) {
v1[(i + j) * 2] = l[j];
v1[(i + j) * 2 + 1] = h[j];
}
}
}
И этой подсказки clang оказалось достаточно чтобы ускорить метод почти в 2 раза! Достаточно хитро он переложил данные в регистры. Попробую описать:
1. За раз считаем 8 произведений чего и хотелось.
2. Самое интересное, до чего я сходу не додумался, а компилятор сделал: раскладывать по регистрам нужно так чтобы итоговые 8 произведений положить назад в V1 двумя movups без перекладывания по кусочкам.
3. Трансформируем V2 и V3 в 4 xmm регистра (делаем всего один раз на весь V1): xmm{i|0..3} = {V2[i], V3[i], V2[i], V3[i]}: для каждого регистра надо 2*movss + 3*unpcklps
4. Трасформируем V1 (снова 2*movss + 3*unpcklps на регистр) и сразу в процессе mulps и addps:
- {V1[0], V1[0], V1[2], V1[2]}, {V1[1], V1[1], V1[3], V1[3]} ... {V1[7], V1[7], V1[9], V1[9]}
Надеюсь я правильно расшифровал ассемблер :)
rdtscp показала следующее:
1. 4 пары произведений: 32 такта для обеих реализаций, но я не думаю что этому можно верить сильно маленькое число и наверное большая погрешность.
2. 100 пар произведений: 884 такта простое решение, 468 после SSE магии.
Железо и софт: i7-4870HQ CPU @ 2.50GHz, OS X, Apple LLVM version 7.3.0 (clang-703.0.29)
Интересно почему результаты по тактам у Вас и у меня на столько разные?
Упс, не нужно было вчера сидеть до 2-х ночи :)
ОтветитьУдалить1. Отступы в коде потерялись :( Но вроде более-менее читаемо...
2. V1 трансформация требует не 3, а только 2 unpcklps на 4 исходных float.
> Помоему это компилятор должен оптимизировать
ОтветитьУдалитьНу так-то да. Но, как выясняется, не всегда.
> Правильно ли я понимаю что массив v1 - огромный, а вот v2 и v3 - единичные векторы?
Да, все так и есть. V2 и V3 - короткие векторы, они всегда постоянны. Но они не единичны.
Замерял я, однако, не расчет для большого массива, а для маленького. Как раз на 4 вектора, правда, требуется вычислить 8 произведений, ну это из-за особенностей вейвлет-преобразования. Сразу не сообразил, можно было и для 4 произведений сделать, для наглядности.
Ага, компилятор реально прикольную оптимизацию провел. Если не трудно, скинь ассемблерный фрагмент, если он читабельный, конечно.
А вот почему настолько разные такты получаются - загадка. Может быть, единицы измерения у интеловского и амдешного процессора разные? 100 пар дает примерно такой же результат, как и одиночный замер. Вообще, 32 такта - очень мало. При расчете 8 векторных произведений требуется произвести 32 умножения и 24 сложения, плюс пересылки. Тут либо интеловский процессор более чем в 10 раз быстрее амдешного. Либо ошибки в измерении у одного из нас. )
> Замерял я, однако, не расчет для большого массива, а для маленького. Как раз на 4 вектора
УдалитьВ таком случае львиная доля уходит на трансформацию V2 & V3, которые хоть для 4-х векторов, хоть для 10К - один раз делать. SSE все таки на поток рассчитаны...
> Если не трудно, скинь ассемблерный фрагмент, если он читабельный, конечно.
Вроде более-менее читабельно: http://pastebin.com/VxrFGeut там считаю 10К пар произведений.
> А вот почему настолько разные такты получаются - загадка.
Это точно. Может у меня и не такты вовсе... Давайте реальным временем померяемся :) У меня на расчет 10К пар произведений (len = 10000 в моих методах) и расчет их в цикле 10000 раз (дабы увеличить время) заняло 0.23 сек и 0.12 сек соответственно.
Как-то коряво я время описал. Я обернул мои методы еще в один цикл, чтобы считать 10К пар произведений 10К раз. Первая реализация расчитываеет за 0.23 сек, а вторая за 0.12 сек.
УдалитьЯ подумал про такты еще маленько :), вроде у меня как раз они получались:
ОтветитьУдалить- 100 пар произведений: 468 тактов
- 10К * 10К пар произведений: 468 * 10^6 тактов
- процессор у меня 2.5GHz, но Turbo Boost частота 3.7GHz. Метод выполняется на одном ядре, не долго, проц. перегреться не должен был успеть, так что я думаю именно с этой частотой все выполнялось
- и того 468 * 10^6 / (3700 * 10^6) = 0.126 сек
- как в аптеке :) по прямому измерению времени у меня выходило именно столько же 0.12 сек.
> При расчете 8 векторных произведений требуется произвести 32 умножения и 24 сложения, плюс пересылки.
У нас ведь SSE :) и того только 2 * (4 + 3) инструкции на вычисление + трансформация V1 и запись результата, плюс как я понимаю за один так может выполнится 2 SSE инструкции или даже больше, см: http://stackoverflow.com/questions/15655835/flops-per-cycle-for-sandy-bridge-and-haswell-sse2-avx-avx2:
Intel Core 2 and Nehalem:
8 SP FLOPs/cycle: 4-wide SSE addition + 4-wide SSE multiplication
На Intel Haswell/Broadwell/Skylake еще круче, но FMA инструкции - это уже другая история :)
32 SP FLOPs/cycle: two 8-wide FMA (fused multiply-add) instructions
Удивил гугл, ага... Неделю назад писал комментарий, а коварный гугл его потерял. Придется восстанавливать по памяти...
УдалитьХотел уточнить. 100 пар произведений = 100*10К? Не очень понял, как во второй строчке появился коэффициент 10^6.
> У нас ведь SSE :) и того только 2 * (4 + 3) инструкции на вычисление + трансформация V1 и запись результата, плюс
> как я понимаю за один так может выполнится 2 SSE инструкции или даже больше, см: http://stackoverflow.com/questions
> /15655835/flops-per-cycle-for-sandy-bridge-and-haswell-sse2-avx-avx2:
Так-то все правильно. Но! Сколько инструкций выполнится за такт, зависит от процессора. Если процессор реально выполняет 2 инструкции за такт, то значит у него должно быть 8 АЛУ для плавающей точки. Я думаю, далеко не все процессоры могут этим похвастаться.
Так же не стоит забывать и про латентность.
Замутил у себя тест на время тоже. Но только пока не стал делать размер вектора 10К. Дело в том, что смысл Wavelet преобразования не в ВП длинных векторов. А в коротких ВП, но их много.
Что было какое-то основание для сравнения, заменил умножение 10К векторов на 2500 векторов длиной 4 элемента. У меня получились вот такие цифры: 0.56 с для FPU (развернутый цикл) и 0.39 c для моего варианта SSE.
Не совсем корректно, конечно, так сравнивать, поэтому еще неделю назад хотел переделать для сравнения на 10К элементные векторы, взяв за основу твой ассемблерный код, но так руки и не дошли. Отвлекся пока на андроид немного.
> Неделю назад писал комментарий, а коварный гугл его потерял.
УдалитьБывает, eventual consistency и все такое, через годик-другой сообщение найдется и добавится :)
> Не очень понял, как во второй строчке появился коэффициент 10^6.
Это чтобы Вас запутать :)
Когда я такты считал, то я рассчитывал 1 раз для 100 пар произведений:
- len = 100 в моих методах что я писал выше
- V1: 2*100 + 2 = 202 элемента, если я правильно понял Ваш код, то для расчета каждого нового вектора мы смещаемся на 2, а не на 4 элемента, поэтому у меня такая математика.
- V2, V3 всегда по 4 элемента
А когда считал время то я зачем-то, наверное чтобы запутать, стал считать 10К пар произведений (размерность V1: 2 * 10K + 2) и еще обернул это в цикл по 10К элементам.
И того: 100 пар * 10^6 = 10K пар * 10K раз :)
> Так-то все правильно. Но! Сколько инструкций выполнится за такт, зависит от процессора. Если процессор реально выполняет 2 инструкции за такт, то значит у него должно быть 8 АЛУ для плавающей точки.
Я вот явно про 8 АЛУ не нашел, но нашел что у моего Haswell полно разных 256-bit unit-ов (как по русски-то, блоков?) которые за один такт переваривают 256-битные регистры (скажем с 8 float-ми), а 128-битные и подавно. Вот дока с картинками и описанием интересная: http://www.realworldtech.com/haswell-cpu/4
> Так же не стоит забывать и про латентность.
Это да, и судя по месиву инструкций clang активно борется с этим, но на сколько успешно не знаю.
> Дело в том, что смысл Wavelet преобразования не в ВП длинных векторов. А в коротких ВП, но их много.
Напишите реалистичный пример когда считается чего-нибудь много и я попробую скормить его clang :)
> заменил умножение 10К векторов на 2500 векторов длиной 4 элемента
Что-то я не понял :) Можете код выложить.
> Я вот явно про 8 АЛУ не нашел
УдалитьНу, они это на приведенном тобою рисунке называют 256-bit unit. Один такой модуль как раз и составляет 8 ALU, правда, весьма специализированных.
На приведенном рисунке таких приведено 6. Некоторые, правда, не показаны для экономии места. Тем не менее, судя по рисунку, процессор может параллельно выполнять одно векторное умножение, а вот два - уже нет. Кроме того, утверждение, что эти регистры переваривают за один такт 256 бит не совсем точное. Они переваривают за один такт, как я понял, только если конвейер полностью заполнен, а операнды уже находятся во внутренних регистрах (явных или неявных) и независимы друг от друга. В противном случае операции могут выполнятся чуток помедленнее.
> Напишите реалистичный пример когда
Легко. Что-то типа
var
Data : array[0..DataLen-1] of single;
i : integer;
begin
i := 0;
while i < DataLen-6 do
begin
DP4x4InLineSimple(@Data[i], V1, V2);
// процедура, выполняющая векторное произведение Data[i] размером в 4 элемента на V1 и V2, а затем
// Data[i+2], также на V1 и V2
i := i+4;
end;
Код примерный, не отлаживал.
>> заменил умножение 10К векторов на 2500 векторов длиной 4 элемента
> Что-то я не понял :) Можете код выложить.
Ага. Тут я был косоязычен. ) Имелось в виду, что заменил умножение одного вектора в 10К элементов на 2500 умножений вектора в 4 элемента.
> Ну, они это на приведенном тобою рисунке называют 256-bit unit.
УдалитьАга это я понял, но что там внутри я не знаю. Не думаю что там прям 8 независимых сумматора для float и еще 4 для double. Скорее 1 хитрый который может в разных режимах работать.
> процессор может параллельно выполнять одно векторное умножение, а вот два - уже нет.
Вроде может 2 умножения или умножение + сложение: "One interesting point is that while Haswell can execute two 256-bit FMAs or FMULs per cycle, there is still only a single 256-bit FADD."
> Они переваривают за один такт, как я понял, только если конвейер полностью заполнен
Согласен, плюс на реальных задачах скорее всего память будет узким местом.
> Код примерный, не отлаживал.
Вроде я как раз этот "реалистичный" пример и считал выше, Ваш DataLen = 2 * мой len.
> Имелось в виду, что заменил умножение одного вектора в 10К элементов на 2500 умножений вектора в 4 элемента.
Я не смог clang перехитрить, он упорно трансформировал V1 и V2 векторы один раз и месил их в цикле :), скорее всего это именно ключевой момент оптимизации.
Наткнулся на интересную статью http://sci.tuomastonteri.fi/programming/sse там есть несколько примеров как лучше данные разложить, может что-то полезное найдете в ней :)
> Не думаю что там прям 8 независимых сумматора для float и еще 4 для double.
УдалитьИ я не знаю, чего там. Но блоки у них разбиты по функциональному назначению. Например, один блок может выполнять только сдвиги и умножения. А другой - сложения. И т. д.
> Вроде может 2 умножения или умножение + сложение
Я судил по количеству блоков векторного умножения, которых на приведенном рисунке вроде один. Хотя, возможно, я что-то не правильно понял.
> Вроде я как раз этот "реалистичный" пример и считал выше
Да? Ну тогда, я значит, в коде просто не разобрался. )
Спасибо за статью. Обязательно гляну. Но сейчас что-то со временем совсем напряг.
Извиняюсь за очередной коммент от меня, ну уж очень у Вас задачка интересная :)
ОтветитьУдалитьВ Linux есть Perf Tool. Этот инструмент анализирует разные CPU счетчики и даже позволяет профилировать ассемблер.
Позапускал я его для задачки все того же размера: 10K итераций по подсчету 20K произведений (другими словами размерность массива 20K + 2) на Xeon W3565 (не Haswell!) и собрал статистику для non-SSE / SSE версий (“perf stat” команда):
- task-clock (msec): 299 / 129
- cycles: 1,005,987,980 / 433,782,823
- instructions: 2,890,554,213 / 1,167,327,335
- insns per cycle: 2.87 / 2.69
Получается даже 4-х летний Xeon выполняет несколько SSE инструкций за один такт!
Еще немного математики:
* non-SSE:
- на 1 векторное произведение: 2,890,554,213 / (10,000 * 20,000) ~ 14 инструкций.
- в теории нужно: 4 (загрузки в память из V1) + 4 (умножение) + 3 (сложение) + 1 (запись результата) = 12, что очень близко к 14-ти.
* SSE:
- на 4 векторных произведения: 1,167,327,335 / (10,000 * 20,000) * 4 ~ 23 инструкции.
- в теории (по сгенерированному clang коду): (2 загрузки + 2 упаковки) * 4 + 4 (умножение) + 3 (сложение) + 1 (запись результата) = 24, опять-же почти 23.