Знайомтеся, loop fracking

image
Метою даної роботи є позначення ще однією техніки оптимізації циклів.
При цьому немає завдання орієнтуватися на якусь існуючу архітектуру, а, навпаки,
будемо намагатися діяти по можливості абстрактно, спираючись переважно на здоровий глузд.

Автор назвав цю техніку «loops fracking» за аналогією з, наприклад, "loops unrolling
або «loops nesting». Тим більше, що термін відображає сенс і не зайнятий.


Цикли є основним об'єктом для оптимізації, саме в циклах більшість програм проводять основний час.
Існує достатня кількість технік оптимізації, ознайомитися з ними можна тут.

Основні ресурси для оптимізації це:
  1. Економія на логіці закінчення циклу. Перевірка критерію закінчення циклу призводить до розгалуження, розгалуження «ламає конвеєр», давайте перевіряти рідше. В результаті з'являються прекрасні зразки коду, наприклад, Duf's device.
    void send(int *to, int *from, int count)
    {
    int n = (count + 7) / 8;
    switch (count % 8) {
    case 0: do {
    *to = *from++;
    case 7: *to = *from++;
    case 6: *to = *from++;
    case 5: *to = *from++;
    case 4: *to = *from++;
    case 3: *to = *from++;
    case 2: *to = *from++;
    case 1: *to = *from++;
    } while (--n > 0);
    }
    }
    
    На даний момент провісники переходів (там, де вони є) в процесорах зробили подібну оптимізацію малоефективною.
  2. Винос інваріанти циклу за дужки (hoisting).
  3. Оптимізація роботи з пам'яттю для більш ефективної роботи кеша пам'яті. Якщо в циклі є звернення до кількості пам'яті, яка явно перевищує обсяг кеш-пам'яті, стає важливим, в якому порядку йдуть ці звернення. Компілятор складно крім очевидних випадків розбиратися з цим, іноді для досягнення ефекту фактично пишеться інший алгоритм. Тому дана оптимізація лягає на плечі прикладних програмістів. А компілятор / профілювальник можуть забезпечити статистикою, дати підказки… зворотний зв'язок.
  4. Використання (явного чи неявного) паралелізму процесора. Сучасні процесори вміють виконувати код паралельно.

    У випадку явно паралельної архітектури (EPIC, VLIW) одна інструкція може містити декілька команд (які будуть виконуватися паралельно), які зачіпають різні функціональні блоки.

    Суперскалярні процесори самостійно розбирають потік інструкцій, вишукують паралелізм і по можливості його використовують.
    image
    Принципова схема суперскалярного виконання команд

    Ще один варіант — векторні операції, SIMD.
Зараз ми займемося пошуком способу максимального використання паралелізму процесора(ів).

Що ми маємо.
Для початку розберемо кілька простих прикладів, використовуючи для дослідів MSVS-2013(x64) процесор Intel Core-i7 2600. До речі, GCC вміє все те ж і не гірше, у всякому разі на таких простих прикладах.

Найпростіший цикл — обчислення суми цілочисельного масиву.
int64_t data[100000];
... 
int64_t sum = 0;
for (int64_t val : data) {
sum += val;
}

Ось що робить з цього компілятор:
lea rsi,[data] 
mov ebp,186A0h ;100 000
mov r14d,ebp 
...
xor edi,edi 
mov edx,edi 
nop dword ptr [rax+rax] ; вирівнювання

loop_start:
add rdx,qword ptr [rsi] 
lea rsi,[rsi+8] 
dec rbp 
jne loop_start 

Те ж, але з double'ами (AVX, /fp:precise & /fp:strict — ANSI сумісність )
vxorps xmm1,xmm1,xmm1 
lea rax,[data] 
mov ecx,186A0h 
nop dword ptr [rax+rax] 

loop_start:
vaddsd xmm1,xmm1,mmword ptr [rax] 
lea rax,[rax+8] 
dec rcx 
jne loop_start
Мільйон разів цей код виконується за 85 сек.

Не видно тут ніякої роботи компілятора по виявленню паралельності, хоча,
здавалося б, у задачі вона в наявності. Компілятор виявив залежність даних і не зміг її обійти.

Те ж, але (AVX, /fp:fast — без ANSI — сумісності)
vxorps ymm2,ymm0,ymm0 
lea rax,[data] 
mov ecx,30D4h ; 12500, 1/8
vmovupd ymm3,ymm2 

loop_start:
vaddpd ymm1,ymm3,ymmword ptr [rax+20h] ; SIMD
vaddpd ymm2,ymm2,ymmword ptr [rax] 
lea rax,[rax+40h] 
vmovupd ymm3,ymm1 
dec rcx 
jne loop_start

vaddpd ymm0,ymm1,ymm2 
vhaddpd ymm2,ymm0,ymm0 
vmovupd ymm0,ymm2 
vextractf128 xmm4,ymm2,1 
vaddpd xmm0,xmm4,xmm0 
vzeroupper 
Виконується 26 сек, використовуються векторні операції.

Той же цикл, але в традиційному стилі
for (i = 0; i < 100000; i ++) {
sum += data[i];
}

Отримуємо несподівано для (/fp:precise)
vxorps xmm4,xmm4,xmm4 
lea rax,[data+8h] 
lea rcx,[piecewise_construct+2h] 
vmovups xmm0,xmm4 
nop word ptr [rax+rax] 

loop_start:
vaddsd xmm0,xmm0,mmword ptr [rax-8] 
add rax,50h 
vaddsd xmm1,xmm0,mmword ptr [rax-50h] 
vaddsd xmm2,xmm1,mmword ptr [rax-48h] 
vaddsd xmm3,xmm2,mmword ptr [rax-40h] 
vaddsd xmm0,xmm3,mmword ptr [rax-38h] 
vaddsd xmm1,xmm0,mmword ptr [rax-30h] 
vaddsd xmm2,xmm1,mmword ptr [rax-28h] 
vaddsd xmm3,xmm2,mmword ptr [rax-20h] 
vaddsd xmm0,xmm3,mmword ptr [rax-18h] 
vaddsd xmm0,xmm0,mmword ptr [rax-10h] 
cmp rax,rcx 
jl loop_start
Тут немає паралельності, всього лише спроба заощадити на обслуговуванні циклу.
Виконується цей код 87 с.
Для /fp:fast код не змінився.

Підкажемо компілятору, використовуючи loop nesting.
double data[100000];
... 
double sum = 0, sum1 = 0, sum2 = 0;
for (int ix = 0; i < 100000; i+=2)
{
sum1 += data[i];
sum2 += data[i+1];
}
sum = sum1 + sum2;
Одержуємо те, що просили, причому код однаковий для опцій /fp:fast & /fp:precise.
Операції vaddsd на деяких процесорах (Бульдозер від АМД) можуть виконуватися паралельно.
vxorps xmm0,xmm0,xmm0 
vmovups xmm1,xmm0 
lea rax,[data+8h] 
lea rcx,[piecewise_construct+2h] 
nop dword ptr [rax] 
nop word ptr [rax+rax] 

loop_start:
vaddsd xmm0,xmm0,mmword ptr [rax-8] 
vaddsd xmm1,xmm1,mmword ptr [rax] 
add rax,10h 
cmp rax,rcx 
jl loop_start
Мільйон разів цей код виконується за 43 сек, вдвічі швидше, ніж «наївний і точний» підхід.

При кроці в 4 елемента код виглядає так (знову ж таки, однаково для опцій компілятора /fp:fast & /fp:precise)
vxorps xmm0,xmm0,xmm0 
vmovups xmm1,xmm0 
vmovups xmm2,xmm0 
vmovups xmm3,xmm0 
lea rax,[data+8h] 
lea rcx,[piecewise_construct+2h] 
nop dword ptr [rax] 

loop_start:
vaddsd xmm0,xmm0,mmword ptr [rax-8] 
vaddsd xmm1,xmm1,mmword ptr [rax] 
vaddsd xmm2,xmm2,mmword ptr [rax+8] 
vaddsd xmm3,xmm3,mmword ptr [rax+10h] 
add rax,20h 
cmp rax,rcx 
jl loop_start 

vaddsd xmm0,xmm1,xmm0 
vaddsd xmm1,xmm0,xmm2 
vaddsd xmm1,xmm1,xmm3 
Цей код мільйон разів виконується за 34 сек.

Для того ж, щоб гарантувати векторні обчислення, доведеться вдаватися до різним трюкам, як то
  1. писати підказки компілятору у вигляді прагм: #pragma ivdep #pragma loop(ivdep) #pragma GCC ivdep), #pragma vector always, #pragma omp simd… і сподіватися що компілятор не проігнорує їх
  2. intrinsic'і — вказівки компілятору, які інструкції використовувати, наприклад, додавання двох масивів може виглядати так:
Якось все це не дуже в'яжеться зі світлим образом «мови високого рівня».
З одного боку, при необхідності отримати результат, зазначені оптимізації зовсім не в тягар.
З іншого, виникає проблема переносимості. Припустимо, програма писалася і відлагоджувати для процесора з чотирма сумматорами.
Тоді при спробі виконати її на версії процесора з 6-ю сумматорами, ми не отримаємо очікуваного виграшу.
А на версії з трьома — отримаємо уповільнення в 2 рази, а не на чверть, як хотілося б.

Наостанок обчислимо суму квадратів (/fp:precise)
vxorps xmm2,xmm2,xmm2 
lea rax,[data+8h] ; pdata = &data[1]
mov ecx,2710h ; 10 000
nop dword ptr [rax+rax] 

loop_start:
vmovsd xmm0,qword ptr [rax-8] ; xmm0 = pdata[-1]
vmulsd xmm1,xmm0,xmm0 ; xmm1 = pdata[-1] ** 2
vaddsd xmm3,xmm2,xmm1 ; xmm3 = 0 + pdata[-1] ** 2 ; sum 
vmovsd xmm2,qword ptr [rax] ; xmm2 = pdata[0]
vmulsd xmm0,xmm2,xmm2 ; xmm0 = pdata[0] ** 2
vaddsd xmm4,xmm3,xmm0 ; xmm4 = sum + pdata[0] ** 2 ; sum
vmovsd xmm1,qword ptr [rax+8] ; xmm1 = pdata[1]
vmulsd xmm2,xmm1,xmm1 ; xmm2 = pdata[1] ** 2
vaddsd xmm3,xmm4,xmm2 ; xmm3 = sum + pdata[1] ** 2 ; sum
vmovsd xmm0,qword ptr [rax+10h] ; ...
vmulsd xmm1,xmm0,xmm0 
vaddsd xmm4,xmm3,xmm1 
vmovsd xmm2,qword ptr [rax+18h] 
vmulsd xmm0,xmm2,xmm2 
vaddsd xmm3,xmm4,xmm0 
vmovsd xmm1,qword ptr [rax+20h] 
vmulsd xmm2,xmm1,xmm1 
vaddsd xmm4,xmm3,xmm2 
vmovsd xmm0,qword ptr [rax+28h] 
vmulsd xmm1,xmm0,xmm0 
vaddsd xmm3,xmm4,xmm1 
vmovsd xmm2,qword ptr [rax+30h] 
vmulsd xmm0,xmm2,xmm2 
vaddsd xmm4,xmm3,xmm0 
vmovsd xmm1,qword ptr [rax+38h] 
vmulsd xmm2,xmm1,xmm1 
vaddsd xmm3,xmm4,xmm2 
vmovsd xmm0,qword ptr [rax+40h] 
vmulsd xmm1,xmm0,xmm0 
vaddsd xmm2,xmm3,xmm1 ; xmm2 = sum;
lea rax,[rax+50h] 
dec rcx 
jne loop_start
Компілятор знову пиляє цикл на шматки по 10 елементів з метою економії на логіці циклу,
при цьому обходиться 5 регістрами — один під суму і по парі на кожну з двох паралельних гілок множень.

Або ось так для /fp:fast
vxorps ymm4,ymm0,ymm0 
lea rax,[data] 
mov ecx,30D4h ;12500 1/8

loop_start:
vmovupd ymm0,ymmword ptr [rax] 
lea rax,[rax+40h] 
vmulpd ymm2,ymm0,ymm0 ; SIMD
vmovupd ymm0,ymmword ptr [rax-20h] 
vaddpd ymm4,ymm2,ymm4 
vmulpd ymm2,ymm0,ymm0 
vaddpd ymm3,ymm2,ymm5 
vmovupd ymm5,ymm3 
dec rcx 
jne loop_start

vaddpd ymm0,ymm3,ymm4 
vhaddpd ymm2,ymm0,ymm0 
vmovupd ymm0,ymm2 
vextractf128 xmm4,ymm2,1 
vaddpd xmm0,xmm4,xmm0 
vzeroupper 


Зведена таблиця:
MSVC,/fp:strict, /fp:precise, с MSVC,/fp:fast, с
foreach 85 26
C style цикл 87 26
C style nesting X2 43 43
C style nesting X4 34 34
Як пояснити ці числа?
Варто відзначити, що справжню картину того, що відбувається, знають тільки розробники процесора, ми ж можемо тільки робити припущення.

Перша думка, що прискорення відбувається за рахунок декількох незалежних суматорів, по всій видимості помилкова.
Процесор i7-2600 має один векторний суматор, який не вміє виконувати незалежні скалярні операції.

Тактова частота процесора — до 3.8 гГц.
85 сек простого циклу дають нам (мільйон разів по 100 000 складань) ~3 такту на ітерацію. Це добре узгоджується з даними (1, 2 про 3 тактах на виконання векторної інструкції vaddpd (нехай навіть ми складаємо скаляри). Оскільки є залежність за даними, швидше 3 тактів ітерацію виконати неможливо.

У разі нестинга (Х2) відсутня залежність за даними всередині ітерації і є можливість завантажити конвеєр суматора з різницею в такт. Але в наступній ітерації залежності за даними виникають так само з різницею в такт, в результаті маємо прискорення в 2 рази.

У разі нестинга (Х4) конвеєр суматора завантажується також з кроком в такт, але прискорення в 3 рази (з-за довжини конвеєра) не відбувається, втручається додатковий фактор. Наприклад, ітерація циклу перестало вміщатися в рядку кешу L0m і отримує штрафний такт(и).

Отже:
  • з моделлю компіляції /fp:fast найпростіший вихідний текст дає найшвидший варіант коду, що логічно оскільки ми маємо справу з мовою високого рівня.
  • ручні оптимізації в дусі nesting'а дають хороший результат для моделі /fp:precise, але тільки заважають компілятор при використанні /fp:fast
  • ручні оптимізації перенесені в більшій мірі, ніж векторизованный код


Трохи про компілятори.
Реєстрові архітектури дають простий і універсальний спосіб отримати прийнятний код
з стерпного тексту на мовах високого рівня.
Компіляцію умовно можна розбити на кілька кроків:
  1. Синтаксичний аналіз. На цьому етапі здійснюється синтаксично керована трансляція, виробляються статичні перевірки. На виході ми маємо дерево (DAG) синтаксичного розбору.
  2. Генерація проміжного коду. При бажанні, генерація проміжного коду може бути об'єднана з синтаксичним аналізом.
    Тим більше, що при використанні в якості проміжного коду трехадресных інструкцій, даний крок стає тривіальним бо "трехадресный код є линеаризованным поданням синтаксичного дерева або дагу, в якому внутрішніх вузлів графа відповідають явні імена".
    В сутності, трехадресный код призначений для віртуального процесора з нескінченною кількістю регістрів.
  3. Генерація коду. Результатом даного кроку є програма для цільової архітектури. Оскільки реальна кількість регістрів обмежено і невелика, на цьому етапі ми повинні визначити, які тимчасові змінні в кожен момент будуть перебувати в регістрах, і розподілити їх по конкретних регістрів. Навіть у чистому вигляді ця задача є NP-повною, крім того, справа ускладнюється тим, що існують різноманітні обмеження по використанню регістрів. Тим не менш, для розв'язання даної задачі розроблено прийнятні евристики. Крім того, трехадресный (або його еквіваленти) код дає нам формальний апарат для аналізу потоків даних, оптимізації, видалення непотрібного коду,…
Вимальовуються проблеми:
  1. Для розв'язання NP-повної задачі розміщення регістрів використовують евристики і це дає код прийнятної якості. Ці евристики не люблять додаткових обмежень на використання пам'яті або регістрів. Наприклад, розщепленої (interlaced) пам'яті, неявного використання регістрів інструкціями, векторних операцій, регістрових перснів… Не люблять аж до того, що евристика може перестати працювати і побудова близького до оптимального коду перестає бути завданням, розв'язуваної універсальним чином.
    В результаті (векторні?) можливості процесора можуть бути задіяні лише тоді, коли компілятор розпізнав типову ситуацію з набору, яким його навчили.
  2. Проблема масштабування. Розміщення (allocation) регістрів робиться статично, якщо ми спробуємо виконати скомпільований код на процесорі з тією ж системою команд, але великою кількістю регістрів, ніякого виграшу не вийде.

    Це стосується і SPARC з його стеком регістрових вікон, для якого більшу кількість регістрів призводить лише до того, що більша кількість фреймів викликів вміщується в пулі регістрів, що зменшує частоту звернень до пам'яті.

    EPIC — зроблена спроба в напрямку масштабування — “Кожна група з декількох інструкцій називається bundle. Кожен bundle може мати стоповий біт, що позначає, що наступна група залежить від результатів даної роботи. Такий біт дозволяє створювати майбутні покоління архітектури з можливістю паралельного запуску декількох bundles. Інформація про залежності обчислюється компілятором, і тому апаратурі не доведеться проводити додаткову перевірку незалежності операндів." Передбачалося, що незалежні bundles можуть виконуватися паралельно, причому, чим більше виконуючих пристроїв в системі, тим ширше може використовуватися внутрішній паралелізм програми. З іншого боку, далеко не завжди такі особливості дають виграш, у всякому разі, для підсумовування масиву це бачиться автору марним.

    Суперскалярні процесори вирішують проблему, вводячи «регістри для нас» і «регістри в собі». На перших орієнтується компілятор, коли розписує(аллоцирует) регістри. Число друге може бути будь-яким, зазвичай їх в рази більше, ніж перших. У процесі декодування суперскалярный процесор заново розписує регістри виходячи з їх фактичної кількості в межах вікна в тілі програми. Розмір вікна визначається складністю логіки, яку процесор може собі дозволити. Крім числа регістрів, звичайно, масштабування підлягають і функціональні пристрої.
  3. Проблема сумісності. Особливо відзначимо X84-64 і лінійку технологій — SSE — SSE2 — SSE3 — SSSE3 — SSE4 — AVX — AVX2 — AVX512 —…

    Сумісність зверху вниз (тобто код скомпільований під старшу технологію, а хочеться виконувати на молодших процесорах) може бути забезпечена одним способом — генерацією коду під кожну згадану технологію і вибір runtime потрібної гілки для виконання. Звучить не дуже привабливо.

    Сумісність знизу вгору забезпечує процесор. Ця сумісність гарантує исполнимость коду, але не обіцяє його ефективного виконання. Наприклад, якщо код скомпільований під технологію з двома незалежними сумматорами, а виконується на процесорі з 4-ма, фактично будуть використані лише два з них. Генерація декількох гілок коду під різні технології не вирішує проблеми майбутніх технологій, планованих чи ні.


Погляд на цикли
Розглянемо ту ж задачу, підсумовування масиву. Уявімо це підсумовування обчислення одного виразу. Оскільки ми використовуємо бінарне додавання, вираз може бути представлено у вигляді бінарного дерева, причому за асоціативності підсумовування, таких дерев багато.

Обчислення проводиться при обході дерева вглиб зліва направо.
Звичайне підсумовування виглядає як ліво-зростаюче вироджене до списку дерево:
image
double data[N];
...
double sum = 0;
for (int i = 0; i < N; i++)
{
sum += data[i];
}

Максимальна глибина стека (обхід в глибину передбачає постфиксное підсумовування, тобто стек), яка тут може знадобитися — 2 елемента. Ніякої паралельності не передбачається, кожне підсумовування (крім першого) має дочекатися результатів попереднього підсумовування, тобто очевидна залежність за даними.
Зате нам достатньо 3 регістрів (сума і два для емуляції верхівки стека) щоб підсумувати масив будь-якого розміру.

Нестинг циклу в два потічка виглядає так:
image
double data[N];
...
double sum = 0;
double sum1 = 0, sum2 = 0;
for (int i = 0; i < N/2; i+=2)
{
sum1 += data[i];
sum2 += data[i + 1];
}
sum = sum1 + sum2;
Для обчислень потрібно вдвічі більше ресурсів, 5 регістрів на все,
але тепер частина підсумовувань можна робити паралельно.

Найжахливіший з обчислювальної точки зору варіант — право-зростаюче, вироджене до списку дерево, для його обчислення потрібен стек розміром з масив при відсутності паралелізму.

Який варіант дерева дасть максимальний паралелізм? Очевидно, збалансований (в межах можливого) дерево, де звернення до вихідними даними є тільки в листових підсумовуючих вузлах.
image
// Алгоритм обчислення суми такий:
double data[N];
for (unsigned ix = 0; ix < N; ix++) {
unsigned lix = ix;
push(data[ix]);
while (lix & 1) {
popadd();
lix >>= 1;
}
}
for (unsigned ix = 1; ix < bit_count(N+1); ix++) {
popadd();
}

У цьому псевдо-коді використані наступні функції:
  1. push(val) — поміщає значення на вершину стека, збільшує стек. Передбачається, що стек організований на пулі регістрів.
  2. popadd() — підсумовує два елемента на вершині стека, поміщає результат в другій зверху, видаляє верхній елемент.
  3. bit_count(val) — підраховує кількість біт в цілочисельному значенні
Після роботи цього псевдокода в стеку залишається єдиний елемент, рівний шуканої суми.

Як це працює? Звернемо увагу, що номер елемента в двійковому поданні кодує шлях до нього від вершини дерева виразу (від старших до молодших бітів). При цьому 0 позначає рух ліворуч, 1 — праворуч (аналогічно коду Хаффмана).

Зазначимо, що кількість безперервно йдуть взведенных молодших біт дорівнює кількості підсумовувань, які треба зробити для обробки поточного елемента. А загальна кількість взведенных біт в номері означає кількість елементів знаходяться в стеку до початку роботи з цим елементом.

На що слід звернути увагу:
  • Розмір стека тобто потрібну кількість регістрів для нього є log2 від розміру даних. З одного боку, це не дуже зручно т. к. розмір даних може бути величиною обчислюваної, а розмір стека хотілося б визначити при компіляції. З іншого боку, ніхто не заважає компілятору пиляти дані на тайли розмір яких визначається виходячи з числа доступних регістрів.
  • У такій інтерпретації завдання досить паралелізму, щоб автоматично задіяти будь-яку кількість доступних незалежних суматорів. Завантаження елементів з масиву може робитися незалежно і паралельно. Підсумовування одного рівня також виконуються паралельно.
  • з паралелізмом є проблеми. Нехай на момент компіляції у нас було N суматорів. Для того, щоб ефективно працювати з відмінним від N їх числом (заради чого все і затівалося), доведеться задіяти апаратну підтримку
    • для архітектур з явним паралелізмом все непросто. Допоміг би пул суматорів і дозвіл виконувати паралельно кілька незалежних «широких» інструкцій. Для підсумовування береться не конкретний суматор, а перший в черзі. Якщо широка інструкція намагається виконати три підсумовування, з пулу будуть витребувані три суматора. Якщо такого до-ва вільних суматорів немає, інструкція буде блокована аж до їхнього звільнення.
    • для суперскалярних архітектур потрібно відстежувати стан стека. Є унарні (ex: зміна знака) і бінарні операції (Ex: popadd) зі стеком. А також листові операції без залежностей (Ex:push). З останніми найпростіше, їх можна ставити на виконання в будь-який момент. Але якщо у операції є аргумент(и), вона перед виконанням повинна дочекатися готовності своїх аргументів.
      Так, обидва значення для операції popadd повинні знаходитися на вершині стека. Але до того моменту як підійде чергу їх підсумувати, вони вже можуть бути зовсім не на вершині стека. Може статися, що вони фізично розташовані не підряд.
      Виходом буде розміщення (allocation) відповідного регістра з пулу стека в момент постановки інструкції на виконання а не в момент, коли результат готовий і його треба розмістити.
      Наприклад push(data[i]) означає, що з пам'яті треба вираховувати значення і розмістити його в регістр. Номер цього регістра береться з черги вільних регістрів до початку читання. Регістр втечуть з черги.
      Для операції popadd на момент постановки на виконання номери регістрів з аргументами вже відомі, хоча, можливо і не обчислені. При цьому береться новий регістр під результат, а при готовності аргументів popadd виконується, регістри з під аргументів віддаються в чергу вільних.
      Постановка інструкцій на виконання відбувається аж до вичерпання пулу регістрів і поступово просувається вперед.
      Таким чином досягається масштабованість і за розміром пулу регістрів і за кількістю функціональних пристроїв. Платою є підтримання стану стека/черги регістрів і залежностей між інструкціями. Що не виглядає непереборною перешкодою.

  • В режимі ANSI сумісності з-за проблем з точністю для виразів з double'ами формально неможливі оптимізаційні алгебраїчні перетворення. Тому існує режим /fp:fast, який розв'язує компілятору руки. Фактично підсумовування пірамідкою навіть більш дбайливо ставиться до точності в силу того, що відбувається підсумовування (більш ймовірно) порівнянних чисел, а не додавання до наростаючою сумою.
  • Можна реалізувати описану схему на існуючих архітектур, без апаратної підтримки стека? Так, для фіксованого розміру масиву. Компілятор пиляє цикл на тайли розміром, скажімо, 64 елемента і явно розписує стекову машину по регістрах. Заодно буде потрібно і код для решти ступенів 2 менше 64 — для виходу з циклу. Це буде досить громіздко, але спрацює.
    
    lea rax,[data] 
    vxorps xmm6,xmm6,xmm6 ; тут будемо тримати 0
    
    vaddsd xmm0,xmm6,mmword ptr [rax] ; 0
    vaddsd xmm1,xmm0,mmword ptr [rax+8] ; 1
    vaddsd xmm0,xmm6,mmword ptr [rax+10h] ; 1 0
    vaddsd xmm2,xmm0,mmword ptr [rax+18h] ; 1 2
    vaddsd xmm0,xmm1,xmm2 ; 0
    vaddsd xmm1,xmm6,mmword ptr [rax+20h] ; 0 1
    vaddsd xmm2,xmm1,mmword ptr [rax+28h] ; 0 2
    vaddsd xmm1,xmm6,mmword ptr [rax+30h] ; 0 2 1
    vaddsd xmm3,xmm1,mmword ptr [rax+38h] ; 0 2 3
    vaddsd xmm1,xmm2,xmm3 ; 0 1
    vaddsd xmm0,xmm0,xmm1 ; 0 
    
    vaddsd xmm1,xmm6,mmword ptr [rax+40h] ; 0 1
    vaddsd xmm2,xmm1,mmword ptr [rax+48] ; 0 2
    vaddsd xmm1,xmm6,mmword ptr [rax+50h] ; 0 2 1
    vaddsd xmm3,xmm1,mmword ptr [rax+58h] ; 0 2 3 
    vaddsd xmm1,xmm2,xmm3 ; 0 1
    vaddsd xmm2,xmm6,mmword ptr [rax+60h] ; 0 1 2
    vaddsd xmm3,xmm2,mmword ptr [rax+68h] ; 0 1 3
    vaddsd xmm2,xmm6,mmword ptr [rax+70h] ; 0 1 3 2
    vaddsd xmm4,xmm2,mmword ptr [rax+78h] ; 0 1 3 4
    vaddsd xmm2,xmm4,xmm3 ; 0 1 2
    vaddsd xmm3,xmm1,xmm2 ; 0 3
    vaddsd xmm1,xmm0,xmm3 ; 1 
    
    Так може виглядати код для 16 -елементної піраміди.


Що далі?
Ми розібрали знаходження суми масиву — дуже просту задачу. Візьмемо що-небудь складніше.

  • Сума квадратів. Знаходиться принципово також, відмінність в тому, що в листових вузлах дерева вираження значення потрібно звести в квадрат.
  • Скалярний добуток векторів. Аналогічно, в листковому вузлі спочатку знаходиться добуток значень з двох масивів.
  • Додавання векторів. Тут немає залежності за даними і природний паралелізм може бути використаний і без особливих хитрувань.
  • Твір матриць. Твір матриць AB складається з усіх можливих комбінацій скалярних творів вектор-рядків матриці A і вектор-стовпців матриці B.
    • Наївний алгоритм оперує скалярними добутками, які ми вже рахувати вміємо. Кубичен.
    • Алгоритм Штрассена рекурсивно розбиває матриці 2Х2, потім перемножує їх, використовуючи 7 множень замість 8. Розбиття продовжується до тих пір, поки підматриці не стануть малого розміру (< 32 ...128), далі вони перемножуються наївним чином. Має складність imageі страждає нестійкістю.
    • Найкращий досягнутий результатimageі він теж рекурсивен. Втім, це суто теоретичний результат без перспектив реального застосування.
  • Наївна свертка.
    for (int i = 0; i < N; i++) {
    for (int j = 0; j < M; j++) {
    result[i + j] += x[i] * h[M - 1 - j];
    }
    }
    З внутрішнім циклом ми вміємо працювати, загальну оптимізацію краще робити через ДПФ.
  • ДПФ. У наївному вигляді не використовується. Алгоритм БПФ знову ж заснований на рекурсії. Наприклад, data-flow діаграма для популярного алгоритму Кулі-Тьюки (Cooley–Tukey) з основою 2 виглядає так (16-точкове):
    image
    Тут достатньо природного паралелізму щоб вдаватися до особливих зусиль.


Останній приклад досить показовий. В цілях оптимізації рекурсія зазвичай переводиться в ітерації, в результаті типовий текст (основний цикл) виглядає так:
nn = N >> 1;
ie = N;
for (n=1; n<=LogN; n++) {
rw = Rcoef[LogN - n];
iw = Icoef[LogN - n];
if(Ft_Flag == FT_INVERSE) iw = -iw;
in = ie >> 1;
ru = 1.0;
iu = 0.0;
for (j=0; j<in; j++) {
for (i=j; i < N; i+=ie) {
io = i + in;
rtp = Rdat[i] + Rdat[io];
itp = Idat[i] + Idat[io];
rtq = Rdat[i] - Rdat[io];
itq = Idat[i] - Idat[io];
Rdat[io] = rtq * ru - itq * iu;
Idat[io] = itq * ru + rtq * iu;
Rdat[i] = rtp;
Idat[i] = itp;
}
sr = ru;
ru = ru * rw - iu * iw;
iu = iu * rw + sr * iw;
}
ie >>= 1;
}

Що можна зробити в даному випадку? В дусі описаної оптимізації циклів, мабуть, що нічого. Тут може бути корисним описаний апаратний стек, питання цікаве. Втім, це вже зовсім інша історія.

PS: окреме спасибі Таситу Мурки (Felid) за консультації з SIMD і не тільки.

PPS: ілюстрація для заголовка взята з відеоряду до King Crimson — Fracture — Live in Boston 1974.

Джерело: Хабрахабр

0 коментарів

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