Векторизація коду перетворення координат у просторі на Intel® Xeon Phi™ з допомогою низькорівневих інструкцій

Введення
При вирішенні задач моделювання руху об'єктів у тривимірному просторі практично завжди потрібне використання операцій просторових перетворень, пов'язаних з множенням матриць перетворень і векторів. Завдання N тіл ця операція використовується багаторазово завдання для повороту і зсуву тіла відносно початку координат. Матриця просторового перетворення має розмірність 4х4, а розмірність вектора, до якого застосовується перетворення, відповідно 4x1. Розглянемо оптимізацію виконання такої операції з великим числом матриць і векторів під архітектуру Intel® Xeon Phi™.

У сопроцессоре Intel® Xeon Phi™ є блок з 32 регістрів довжиною 512 біт (
zmm00-zmm31
). Відповідно, компілятор С++ від Інтел [1] підтримує набір векторних інструкцій, які працюють з цими регістрами, що потрапляє під специфікацію з кодовою назвою KNC [2] (не плутати з AVX-512). Цей набір інструкцій обмежений, порівняно з інструкцій AVX-512, але тим не менш включає операції, необхідні для вирішення розглянутої задачі. Високорівневі обгортки над машинними інструкціями, які підтримуються компілятором, називаються функціями-интринсиками і описані якісної документації від Інтел (Intel® Intrinsic Guide) [3]. Векторні регістри дозволяють при правильній упаковці даних виконувати деяку операцію відразу над набором даних (SIMD — single instruction multiple data). Далі будемо використовувати назви для регістрових змінних в коді, відповідні імена регістрів (
zmm00-zmm31
), і припускати, що вони відображаються на ті ж самі регістри при компіляції.
У просторових трансформаціях поширені дві операції: множення матриці на матрицю перетворення і множення матриці на вектор. Операцію множення матриць перетворення в більшості випадків можна звести до послідовності оптимізованих операцій множення матриці на вектор-стовпці іншої матриці, тому операція множення на вектор є найбільш цікавою, і множення матриць окремо розглядати не будемо. Щоб зрозуміти, як працює векторизація, розглянемо найпростіший приклад множення матриці 4x4 на вектор 4x1.
В цьому випадку відбувається наступне:
A\times \vec{a}=
\begin{pmatrix}
A_{11}\\
A_{21}\\
A_{31}\\
A_{41}\\
\end{pmatrix}\circ
\begin{pmatrix}
a_{1}\\
a_{1}\\
a_{1}\\
a_{1}\\
\end{pmatrix}+
\begin{pmatrix}
A_{12}\\
A_{22}\\
A_{32}\\
A_{42}\\
\end{pmatrix}\circ
\begin{pmatrix}
a_{2}\\
a_{2}\\
a_{2}\\
a_{2}\\
\end{pmatrix}+
\begin{pmatrix}
A_{13}\\
A_{23}\\
A_{33}\\
A_{43}\\
\end{pmatrix}\circ
\begin{pmatrix}
a_{3}\\
a_{3}\\
a_{3}\\
a_{3}\\
\end{pmatrix}+
\begin{pmatrix}
A_{14}\\
A_{24}\\
A_{34}\\
A_{44}\\
\end{pmatrix}\circ
\begin{pmatrix}
a_{4}\\
a_{4}\\
a_{4}\\
a_{4}\\
\end{pmatrix}
,
де для наочності порожнім гуртком показано поелементне множення векторів.
Якщо помістити в регістри стовпці матриці і розмножити по регістрах елементи вектора, операцію можна описати за допомогою функцій-интринсиков для множення і FMA (fused multiply-add – змішаного множення з додаванням). Тобто для отримання у векторному регістрі
zmm00
вектора-результату треба завантажити стовпці матриці в регістри
zmm00-zmm03
, елементи вектора розмножити в регістри
zmm04-zmm07
, так, щоб у кожному лежав багаторазово один елемент вектора, а потім виконати 3 операції:
zmm00=_mm512_mul_ps(zmm00, zmm04);
zmm00=_mm512_fmadd_ps(zmm01, zmm05, zmm00);
zmm00=_mm512_fmadd_ps(zmm02, zmm06, zmm00);
zmm00=_mm512_fmadd_ps(zmm03, zmm07, zmm00);

Щоб розташувати стовпців матриць в регістрах можна використовувати транспонування, якщо матриця зберігається по рядках, але це трудомістка операція, тому краще всього зберігати матриці по стовпцях. І комбінувати перетворення так, щоб виходили матриці по стовпцях.
У цьому прикладі з речовими 32-бітними числами (суфікс ps) регістри використовуються тільки на чверть, так як в 512-бітний регістр поміщається 16 елементів, а в нашому випадку тільки 4, тобто має місце програш в 4 рази. Якщо використовувати подвійну точність (еквівалентний код з суфіксом pd), то регістри будуть заповнені наполовину. Навіть такі варіанти часткового заповнення регістрів можуть дати прискорення якщо використовувати з 4 елементів векторного регістра з результатом для збереження, так як виконуються відразу 4 операції, замість однієї.
Особливості завантаження і вивантаження даних
Для підвищення ефективності завантаження даних з пам'яті в регістри і вивантаження з регістрів у пам'ять виділена в програмі пам'ять повинна бути вирівняна по межі в 64 байти. В цьому випадку можна використовувати операції
_mm512_load
та
_mm512_store
з відповідними суфіксами для типів даних.
Для динамічного виділення вирівняною пам'яті реалізовані в компіляторі спеціальна функція
_mm_malloc
, яка використовується в зв'язці з
_mm_free
. Від звичайної функції
malloc
відмінність в тому, що доданий другий аргумент – кратність кордону вирівнювання пам'яті, для 512-бітних регістрів повинен повертатися вказівник на пам'ять з адресою, кратним 64 байтам. Для типу Type синтаксис виділення вирівняною пам'яті наступний:
Type* pointer = (Type*)_mm_malloc(sizeof(Type)*N, 64);

Для статичних масивів даних також передбачені директиви вирівнювання пам'яті, наприклад, для визначення вирівняного масиву в компіляторі Інтел передбачений наступний синтаксис:
__declspec(align(64)) Type mem[N];

також Можна використовувати підтримуваний *nix – варіант:
Type mem[N] __attribute__((aligned(64)));

У разі неможливості вирівняти пам'ять необхідно використовувати интринсики для роботи з невыровненной пам'яттю
_mm512_loadu
та
_mm512_storeu
. З незрозумілих причин в компіляторі не реалізовані функції завантаження в регістри з невыровненной пам'яті і вивантаження з регістрів у не вирівняну пам'ять для KNC, для AVX-512 такі функції є.
Для KNC ці функції можна реалізувати самостійно на основі интринсиков для часткової завантаження до кордону вирівнювання:
inline __m512d _mm512_loadu_pd(const double* a)
{
__m512d reg = _mm512_setzero_pd();
reg =_mm512_loadunpacklo_pd(reg, a);
reg =_mm512_loadunpackhi_pd(reg, a+8);

return reg;
}

Аналогічно виглядає функція вивантаження:
inline void _mm512_storeu_pd(double* a, __m512d reg)
{
reg =_mm512_packstorelo_pd(a, reg);
reg =_mm512_packstorehi_pd(a+8,reg);
return reg;
}

Правильне використання вирівняною пам'яті може призводити до прискорення до 3 разів, тому далі будемо розглядати тільки роботу з вирівняною пам'яттю.
Оптимізація зі спеціальним зберіганням матриць в пам'яті
У найбільш загальному випадку для чисел одинарної точності можна помітити, що матриці завантажуються повністю в один регістр. Якщо елементи матриць зберігати послідовно, то при завантаженні отримуємо регістр з 16 елементами однієї матриці. Це добре для операцій векторного множення, але погано для додавання, тому що складати елементи необхідно всередині регістра. Операції додавання елементів всередині регістра пов'язані з повільними інструкціями перемішування, зсуву та перестановки елементів, або редукції, так як операцій горизонтального додавання в наборі інструкцій KNC немає. Якщо написати складний алгоритм, який робить це на интринсиках, виходить більш повільний код, що згенерував компиляторным оптимізатором.
Для мінімізації числа інструкцій нам необхідно, щоб обчислення робилися наступним чином:
\begin{pmatrix}
A\times\vec{a}\\
B\times\vec{b}\\
C\times\vec{c}\\
D\times\vec{d}\\
\end{pmatrix}=
\begin{pmatrix}
A_{11}\\
A_{21}\\
A_{31}\\
A_{41}\\
B_{11}\\
\cdots\\
D_{41}
\end{pmatrix}\circ
\begin{pmatrix}
a_{1}\\
a_{1}\\
a_{1}\\
a_{1}\\
b_{1}\\
\cdots\\
d_{1}
\end{pmatrix}+
\begin{pmatrix}
A_{12}\\
A_{22}\\
A_{32}\\
A_{42}\\
B_{12}\\
\cdots\\
D_{42}\end{pmatrix}\circ
\begin{pmatrix}
a_{2}\\
a_{2}\\
a_{2}\\
a_{2}\\
b_{2}\\
\cdots\\
d_{2}
\end{pmatrix}+
\cdots+
\begin{pmatrix}
A_{14}\\
A_{24}\\
A_{34}\\
A_{44}\\
B_{14}\\
\cdots\\
D_{44}\end{pmatrix}\circ
\begin{pmatrix}
a_{4}\\
a_{4}\\
a_{4}\\
a_{4}\\
b_{4}\\
\cdots\\
d_{4}
\end{pmatrix}
,
Щоб порахувати відразу 4 вектора результату, має сенс розташувати дані в пам'яті так, щоб рядка матриці лежали по черзі, і за один набір інструкцій виконувалося множення відразу 4 матриць на 4 вектора. Для цього необхідно записувати нові елементи матриць в наступному порядку:
A_{11},A_{21},A_{31},A_{41},B_{11},B_{21},B_{31},B_{41},\ldots,D_{14},D_{24},D_{34},D_{44}
Звідси видно, що необхідно зберігати спочатку перші стовпці, потім другі, потім треті і в кінці четверті одразу для 4 матриць. За одну операцію завантаження в регістр будуть поміщатися 16 елементів:
zmm00=(A_{11},A_{21},A_{31},A_{41},B_{11},B_{21},B_{31},B_{41},\ldots,D_{41})
За аналогією вектора теж слід розташувати за першим, другим, третім і четвертим елементів відразу 4 векторів. Але, як правило, це незручно, так як вектора завжди використовуються як послідовності елементів для складання, визначення положення тіл у просторі, в той час як матриці повороту зазвичай обчислюються з тригонометричних перетворень, а орієнтація тіла задається системою кутів (наприклад, кути Ейлера) або кватернионами. Тобто в обчислювальному коді вектора повинні зберігатися послідовно для швидкого звернення до них. Матриці можуть зберігатися довільно, так як вони обчислюються через параметри, що описують орієнтацію, а єдину необхідну операцію з ними ми як раз векторизуем спеціальним чином з урахуванням цього особливого зберігання.
Виходячи з послідовного зберігання векторів розглянемо, як можна отримати необхідні нам послідовності елементів. Тобто, якщо ми мали в пам'яті вектор
vec=(a_1,a_2,a_3,a_4,b_1,b_2,b_3,b_4,c_1,c_2,c_3,c_4,d_1,d_2,d_3,d_4 ),
то необхідно одержати наступний набір регістрів
\begin{array}l
zmm04=(a_1,a_1,a_1,a_1,b_1,b_1,b_1,b_1,c_1,c_1,c_1,c_1,d_1,d_1,d_1,d_1 ),\\
zmm05=(a_2,a_2,a_2,a_2,b_2,b_2,b_2,b_2,c_1,c_2,c_3,c_4,d_1,d_2,d_3,d_4 ),\\
zmm06=(a_3,a_3,a_3,a_3,b_3,b_3,b_3,b_3,c_3,c_3,c_3,c_3,d_3,d_3,d_3,d_3 ),\\
zmm07=(a_4,a_4,a_4,a_4,b_4,b_4,b_4,b_4,c_4,c_4,c_4,c_4,d_4,d_4,d_4,d_4 ).
\end{array}
Використовуючи інструкцію розмноження елементів, ми можемо завантажити в 4 регістра елементи a1, b1, c1, d1:
zmm08 = _mm512_extload_ps(vec + 0, 
_MM_UPCONV_PS_NONE, _MM_BROADCAST_1X16, _MM_HINT_NONE);
zmm09 = _mm512_extload_ps(vec + 4, 
_MM_UPCONV_PS_NONE, _MM_BROADCAST_1X16, _MM_HINT_NONE);
zmm10 = _mm512_extload_ps(vec + 8, 
_MM_UPCONV_PS_NONE, _MM_BROADCAST_1X16, _MM_HINT_NONE);
zmm11 = _mm512_extload_ps(vec + 12, 
_MM_UPCONV_PS_NONE, _MM_BROADCAST_1X16, _MM_HINT_NONE);

Після цього є можливість змішати їх необхідним чином нам з допомогою бітових масок:
zmm06 = _mm512_mask_blend_ps(mask32_0, zmm08, zmm09);
zmm07 = _mm512_mask_blend_ps(mask32_1, zmm10, zmm11);
zmm04 = _mm512_mask_blend_ps(mask32_2, zmm06, zmm07);

Для наочності на малюнку показано, як працюють завантаження і перемішування елементів в регістрах.


Змішувати можна тільки 2 регістру, тому операція повторюється 3 рази з різним значенням маски. Після цього маємо необхідне нам вміст в регістрі
zmm04
. Аналогічні операції завантаження-перемішування дозволяють отримати необхідні елементи векторів в регістрах
zmm05-zmm07
, після чого можна виконувати операції множення та додавання для отримання результату. Такий спосіб можна застосувати на процесорах з підтримкою інструкцій AVX-512 і KNC.
Можна помітити, що елементи векторів необхідно частково дублювати четвірками у відповідних чвертях регістрах. Для цього в наборі інструкцій KNC для Intel® Xeon Phi™ є спеціальна операція
swizzle
, що дозволяє переставляти і розмножувати елементи в чвертях регістрів. У разі завантаження 4 векторів з послідовним проходженням елементів наступний код перетворює його в необхідні множники:
zmm04 = _mm512_swizzle_ps(zmm08, _MM_SWIZ_REG_AAAA);
zmm05 = _mm512_swizzle_ps(zmm08, _MM_SWIZ_REG_BBBB);
zmm06 = _mm512_swizzle_ps(zmm08, _MM_SWIZ_REG_CCCC);
zmm07 = _mm512_swizzle_ps(zmm08, _MM_SWIZ_REG_DDDD);

Як працює
swizzle
, показано на малюнку.


Якщо ми мали в регістрі
zmm08
вектор послідовних елементів vec, то в результаті роботи цих інструкцій, сформуються необхідні нам
zmm04-zmm07
. Останній варіант містить істотно менше інструкцій, але не застосуємо для набору інструкцій AVX-512, тому слід вибирати максимально відповідні інструкції для конкретної архітектури і використовувати їх.
Для демонстрації переваги за часом виконання, отриманого в результаті оптимізації коду, було проведено тестування множень великого числа матриць перетворення на вектори, заповнені випадковими числами. Для порівняння використовувався тест множення 300000 матриць перетворення на сопроцессоре Intel® Xeon Phi™ 5120D.
Порівняння результатів оптимізації проводилося після попереднього тренувального запуску цього ж коду. Після прогріву вибиралося середній час з 10 наступних запусків. Щоб уникнути оптимізації невикористаного коду, також була створена штучна залежність за даними (результати операції записувались у вектор-множник для наступної операції). Результати для 32 і 64-бітних дійсних чисел показані в таблиці 1, де в стовпці авто знаходиться результат для автовекторизации.
Таблиця 1 – Виміри часу і прискорення оптимізованого коду





тип час, з авто час, з blend час, з swizzle прискорення blend прискорення swizzle 32 1.15 0.34 0.2 3.38 5.75 64 1.41 0.55 0.45 2.56 3.13
Оптимізація зі звичайним зберіганням матриць
У деяких випадках матриці впорядкувати не можна, тому необхідно робити їх завантаження особливим чином. Розглянемо алгоритм перетворення матриць для подвійної точності (для одинарної все аналогічно, але дуже громіздко). У випадку подвійної точності нам необхідно в кожен умножаемый регістр з елементами матриць помістити 2 набори по 4 елементи. У наборі інструкцій KNC є операція склеювання зі зрушенням вправо, яка поєднує 2 регістра в послідовність із 1024 біт, після чого зміщує його на задане число елементів вправо і записує в результуючий регістр молодші 512 біт. Операція представлена тільки для цілочисельного типу даних, хоча в інструкціях AVX-512 вона присутня і для речових, але з-за однакового бітового представлення нам тип даних не важливий. Потрібно мати на увазі тільки те, що число сдвигаемых елементів повинно бути в 2 рази більше для чисел подвійної точності. Шлях є дві матриці, які зберігаються як послідовність стовпчиків:
mtx=(A_1,A_2,A_3,A_4,B_1,B_2,B_3,B_4)
Для початку завантажуємо матриці, які зберігаються по колонках послідовно:
__m512i izmm10 = _mm512_load_epi64(pmtx); //A2-A1
__m512i izmm11 = _mm512_load_epi64(pmtx+8); //A4-A3
__m512i izmm12 = _mm512_load_epi64(pmtx+16); //B2-B1
__m512i izmm13 = _mm512_load_epi64(pmtx+24); //B4-B3

В результаті операції
izmm10-izmm11
– перша матриця, в
izmm12-izmm13
– друга. Нам необхідно помістити перші рядки обох матриць в регістр
zmm00
, другі – в
zmm01
, треті – в
zmm02
, четверті – в
zmm03
. Для цього використовуємо склеювання зі зсувом на 256 біт (8 цілих чисел або 4 дійсних числа подвійної точності):
izmm14 = _mm512_alignr_epi32(izmm10, izmm12, 8); //A1-B2
zmm00 = (__m512d)_mm512_alignr_epi32(izmm12, izmm14, 8); //B1-A1
zmm01 = (__m512d)_mm512_alignr_epi32(izmm14, izmm10, 8); //B2-A2
izmm14 = _mm512_alignr_epi32(izmm11, izmm13, 8); //A3-B4
zmm02 = (__m512d)_mm512_alignr_epi32(izmm13, izmm14, 8); //B3-A3
zmm03 = (__m512d)_mm512_alignr_epi32(izmm14, izmm11, 8); //B4-A4

На малюнку показано, як це працює.

В результаті операції у регістрах
zmm00-zmm03
отримуємо необхідну послідовність стовпців матриць. Після цього використання
swizzle
для завантаження векторів в регістри
zmm04-zmm07
і подальше застосування FMA-інструкцій, як це робилося раніше, дозволяє отримати результат для двох операцій множення матриць перетворень на два вектора в регістрі
zmm00
.
zmm04 = _mm512_swizzle_pd(zmm08, _MM_SWIZ_REG_AAAA);
zmm05 = _mm512_swizzle_pd(zmm08, _MM_SWIZ_REG_BBBB);
zmm06 = _mm512_swizzle_pd(zmm08, _MM_SWIZ_REG_CCCC);
zmm07 = _mm512_swizzle_pd(zmm08, _MM_SWIZ_REG_DDDD);

zmm00=_mm512_mul_pd(zmm00, zmm04);
zmm00=_mm512_fmadd_pd(zmm01, zmm05, zmm00);
zmm00=_mm512_fmadd_pd(zmm02, zmm06, zmm00);
zmm00=_mm512_fmadd_pd(zmm03, zmm07, zmm00);

Аналогічно можна реалізувати код для чотирьох операцій множення матриць на чотири вектора з 32-бітними речовими числами. Якщо порівняти виміри часу, можна оцінити втрати, пов'язані з попереднім перетворенням матриці (таблиця 2), проте навіть цей код дає прискорення більше ніж у 2 рази, порівняно з кодом, згенерованим компілятором.
Таблиця 2 – Виміри часу і прискорення з різних зберіганням матриць





зберігання матриць час, з авто ,
swizzle
прискорення swizzle черезстрочное 1.41 0.45 3.13 звичайне 1.25 0.56 2.23
Висновок
Оптимізація коду для просторових перетворень з матрично-векторних множенням для співпроцесора Intel® Xeon Phi™ дозволяє отримати більш швидкий код, порівняно з компиляторной оптимізацією, якщо грамотно використовувати доступний набір векторних інструкцій і організувати ефективне збереження даних в пам'яті.
Посилання
[1]. Intel® C++ Compiler 16.0 User and Reference Guide//
https://software.intel.com/en-us/intel-cplusplus-compiler-16.0-user-and-reference-guide
[2]. Intel® Xeon Phi™ Coprocessor Instruction Set Architecture Reference//
https://software.intel.com/sites/default/files/forum/278102/327364001en.pdf
[3]. Intel® Intrinsic Guide//
https://software.intel.com/sites/landingpage/IntrinsicsGuide/
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.