Чисельні методи та прогноз продуктивності

Паралельне розв’язання нелінійних рівнянь

Корені рівняння $f(x)=0$ на відрізку $[a;b]$ шукають у два етапи:

1. відокремлення коренів: відрізок ділять на $M$ малих відрізків і знаходять ті, на кінцях яких $f$ має різні знаки; на кожному такому відрізку є корінь;
2. уточнення кожного кореня бісекцією, методом хорд або Ньютона.

Обидва етапи ідеально паралельні: значення $f$ у $M+1$ вузлах незалежні (паралелізм даних), а корені уточнюються незалежно (паралелізм задач). Крок сітки має бути меншим за найменшу відстань між коренями: якщо на відрізку опинилися два корені, знаки на кінцях однакові, і обидва корені буде втрачено. У прикладі «Корені рівняння» для многочлена Лежандра $P_{1000}$ сітка з 400 тисяч відрізків знаходить усі 1000 коренів, а зі 100 тисяч – лише 996: біля кінців відрізка $[-1;1]$ корені розташовані ближче за крок $2\cdot {10}^{-5}$.

Бісекція (bisection) на кожній ітерації ділить відрізок навпіл і зберігає половину зі зміною знака; після $k$ ітерацій довжина зменшується в $2^k$ разів (для точності ${10}^{-15}$ на відрізку $5\cdot {10}^{-6}$ – 32 ітерації). Одну бісекцію теж можна розпаралелити: $p$-секційний пошук обчислює $f$ одночасно в $p$ внутрішніх точках і зменшує відрізок у $p+1$ разів за ітерацію. Це вигідно лише для дуже дорогої функції: для дешевої накладні витрати на задачі більші за виграш.

Метод Ньютона $x_{k+1}=x_k-f(x_k)/f^{\prime }(x_k)$ збігається квадратично (кількість правильних цифр подвоюється на кожній ітерації), але лише з достатньо близького наближення. Багатостартовий метод (multistart) запускає Ньютона паралельно з багатьох початкових точок і збирає різні знайдені корені. Для комплексних многочленів множина початкових точок, з яких метод збігається до певного кореня, утворює басейн притягання (basin of attraction) з фрактальними межами – класичний приклад паралелізму даних, де кожна точка площини незалежна. Для систем нелінійних рівнянь $F(x)=0$ метод Ньютона на кожній ітерації обчислює матрицю Якобі (її $n^2$ частинних похідних незалежні – паралельно) і розв’язує лінійну систему (метод Гаусса, тема 7).

Метод спряжених градієнтів

Для великих розріджених систем $Ax=b$ із симетричною додатно визначеною матрицею (дискретизація рівнянь Лапласа й Пуассона, задачі теплопровідності, пружності) найчастіше використовують метод спряжених градієнтів (conjugate gradient, CG). На відміну від методу Якобі (тема 7), він рухається не вздовж залишку, а вздовж напрямків, спряжених відносно $A$ ($d_i^{T}A{d}_j=0$), тому в точній арифметиці знаходить розв’язок не більше ніж за $n$ ітерацій, а на практиці досягає потрібної точності за кількість ітерацій порядку $\sqrt{\kappa }$, де $\kappa$ – число обумовленості матриці. Алгоритм:

1. $x_0=0$, $r_0=b$, $d_0=r_0$;
2. для $k=0,1,\dots$: $q=A{d}_k$; $\alpha =(r_k,r_k)/(d_k,q)$; $x_{k+1}=x_k+\alpha d_k$; $r_{k+1}=r_k-\alpha q$;
3. якщо $||r_{k+1}||<ϵ||b||$ – кінець; інакше $\beta =(r_{k+1},r_{k+1})/(r_k,r_k)$, $d_{k+1}=r_{k+1}+\beta d_k$.

Одна ітерація містить одне множення матриці на вектор, два скалярні добутки і три операції axpy ($y\leftarrow \alpha x+y$). Множення розрідженої матриці на вектор паралельне за рядками (горизонтальні смуги), axpy – повністю паралельні, а скалярні добутки – редукції. Саме редукції обмежують масштабованість: кожна потребує глобальної синхронізації (на кластері – MPI_Allreduce, тема 12), і ітерація не може продовжитися, доки не відоме число $\alpha$ чи $\beta$. Тому варіанти CG для кластерів переставляють операції, щоб об’єднати два скалярні добутки в одну редукцію.

Детермінованість особливо важлива в ітераційних методах: якщо скалярні добутки залежать від порядку завершення потоків, кількість ітерацій може змінюватися від запуску до запуску. У лабораторній роботі (приклад 2) метод спряжених градієнтів для п’ятидіагональної матриці рівняння Пуассона на сітці $500\times 500$ виконує однакові 986 ітерацій і дає побітово однаковий результат для 1–16 потоків, бо скалярні добутки обчислюються за 64 фіксованими частинами.

Системи звичайних диференціальних рівнянь

Задачу Коші для системи $y^{\prime }=f(t,y)$, $y(t_0)=y_0$, де $y$ – вектор з $m$ компонент, найчастіше розв’язують методом Рунге–Кутти 4-го порядку (RK4) з кроком $h$:

$$k_1=f(t,y),k_2=f(t+\frac{h}{2},y+\frac{h}{2}{k}_1),$$$$k_3=f(t+\frac{h}{2},y+\frac{h}{2}{k}_2),k_4=f(t+h,y+h{k}_3),$$$$y(t+h)\approx y+\frac{h}{6}(k_1+2k_2+2k_3+k_4).$$

Похибка на відрізку – $O(h^4)$, її контролюють правилом Рунге (порівняння кроків $h$ і $h/2$) або вкладеними формулами з автоматичним вибором кроку.

Кроки за часом послідовні: $y(t+h)$ залежить від $y(t)$. Тому паралелізм шукають в інших вимірах:

• паралелізм за траєкторіями (параметричні розрахунки, parameter sweep): систему розв’язують для тисяч наборів параметрів чи початкових умов, і кожна траєкторія незалежна. Це найпростіший і найефективніший випадок, але тривалість траєкторій різна (одна зупиняється за секунду, інша – за хвилину), тому потрібне динамічне балансування. На кластері й у гріді кожна група траєкторій стає окремим завданням (масиви завдань Slurm, тема 13);
• паралелізм за компонентами: для великої системи (метод прямих для рівняння теплопровідності, задача $N$ тіл з мільйонами частинок) обчислення $f(t,y)$ ділять між потоками за компонентами $y$. Між стадіями $k_1,\dots ,k_4$ потрібна синхронізація (чотири бар’єри на крок), тому компонент має бути досить багато, щоб робота стадії переважала бар’єр;
• паралелізм за часом (алгоритм Parareal та подібні): грубий послідовний розв’язок уточнюється паралельно на відрізках часу; використовують на дуже великих кластерах, коли інші виміри вичерпано.

У лабораторній роботі (приклад 3) 3072 траєкторії маятника з опором для різних коефіцієнтів опору й початкових швидкостей тривають від кількох до сотень мілісекунд. Блочний розподіл дав прискорення лише 2,3, циклічний – 11,0, а динамічний зі спільним лічильником – 11,9 на 16 логічних процесорах.

Аналітичне прогнозування продуктивності

Закони Амдала й Густафсона (тема 1) оцінюють прискорення за часткою послідовного коду. Для конкретного алгоритму точніший прогноз дає модель часу виконання, складена з етапів PCAM:

$$T_p=T_{\text{обч}}(n,p)+T_{\text{обм}}(n,p)+T_{\text{синхр}}(p)+T_{\text{прост}}(n,p),$$

де $T_{\text{обч}}$ – обчислення найзавантаженішого потоку, $T_{\text{обм}}$ – обміни (на кластері – модель α–β чи LogP, на спільній пам’яті – переміщення даних з пам’яті й між кешами), $T_{\text{синхр}}$ – бар’єри, запуск циклів і задач, $T_{\text{прост}}$ – простої через дисбаланс. Параметри моделі отримують мікровимірюваннями: час однієї операції на одному ядрі, час порожнього паралельного циклу чи бар’єра, пропускну здатність пам’яті або мережі. Прогнозоване прискорення $S_p=T_1/T_p$ порівнюють з виміряним: збіг означає, що модель враховує головні витрати, а розбіжність вказує, що саме модель пропустила.

Приклад: множення матриці на вектор. Для горизонтальних смуг на спільній пам’яті

$$T_p=max(\frac{T_1}{p},\frac{8n^2}{B})+t_{\text{синхр}},$$

де $8n^2$ – байти матриці double, $B\approx 44$ ГБ/с – виміряна в темі 7 пропускна здатність пам’яті, $t_{\text{синхр}}$ – час порожнього Parallel.For (кілька мікросекунд). Для вертикальних смуг і шахової схеми додаються редукція часткових векторів і ще один паралельний цикл. Для $n=8000$ матриця читається з пам’яті щонайменше $512\cdot {10}^6/44\cdot {10}^9\approx 11,6$ мс, а $T_1\approx 54$ мс, тому прогноз $S\le 4,6$ незалежно від кількості потоків – і виміряне прискорення 4,3–4,8 з ним збігається. Для $n=1000$ дані в кеші L3, модель прогнозує майже лінійне прискорення, а виміряно не більше 3,7: 0,8 мс роботи мало, і пробудження потоків пулу, яке модель не враховує, стає помітним.

Приклад: множення матриць. У лабораторній роботі прогноз $T_p=T_1/min(p,8)$ враховує лише 8 фізичних ядер (16 логічних процесорів SMT не подвоюють обчислювальних пристроїв), а для алгоритму Кеннона – ще й виміряний час копіювання блоків і бар’єрів. Цей внесок виявився мізерним (прогноз 7,8 замість 8,0), а виміряне прискорення – 4–5,5 (рис. 8.11). Розбіжність показує витрати, яких модель не має: одне ядро працює на вищій частоті, ніж вісім одночасно, а 9 блоків на 8 ядрах дають дисбаланс відображення.

Знімок екрана

Windows Terminal: dotnet run -c Release of the lab example 1 project (MatMul); the table Схема / p / Час, мс / S / Прогноз S / E / max|ΔC| for смуги, шахова, Кеннон

Рис. 8.11. Прогноз і вимірювання прискорення

Ізоефективність

Ефективність можна записати через накладні витрати $T_o=p{T}_p-T_1$ – сумарний час усіх процесорів, не витрачений на корисну роботу:

$$E=\frac{T_1}{p{T}_p}=\frac{1}{1+T_o/T_1}.$$

Зі збільшенням $p$ накладні витрати зростають (більше обмінів і синхронізацій), і ефективність падає. Але $T_o$ зазвичай зростає повільніше за корисну роботу $T_1=W$, якщо збільшувати розмір задачі. Ефективність залишиться сталою, якщо $W$ зростатиме так, щоб

$$W=K\cdot T_o(W,p),K=\frac{E}{1-E}.$$

Залежність $W(p)$, яка це забезпечує, називають функцією ізоефективності (isoefficiency function, Грама, Гупта, Кумар, 1993). Що повільніше вона зростає, то краще масштабується алгоритм (табл. 8.6).

Таблиця 8.6. Функції ізоефективності типових алгоритмів

Алгоритм (розподілена пам’ять)	Накладні витрати $T_o$	Ізоефективність $W$
сума $n$ чисел деревом	$\mathrm{\Theta }(p\log p)$	$\mathrm{\Theta }(p\log p)$
матриця на вектор, смуги	$\mathrm{\Theta }(np)$	$\mathrm{\Theta }(p^2)$
матриця на вектор, шахова	$\mathrm{\Theta }(n\sqrt{p}\log p)$	$\mathrm{\Theta }(p{\log }^{2}p)$
множення матриць Кеннона	$\mathrm{\Theta }(n^2\sqrt{p})$	$\mathrm{\Theta }(p^{3/2})$

Наприклад, для суми $n$ чисел на $p$ процесорах кожен процес додає $n/p$ чисел, а потім ${\log }_{2}p$ кроків редукції; $T_o\approx 2p{\log }_{2}p$. Щоб зберегти ефективність 80 % ($K=4$), потрібно $n\approx 8p{\log }_{2}p$: для 16 процесорів – 512 чисел, для 1024 – 82 тисячі. Для множення матриці на вектор смугами $W=n^2\approx p^2$, тобто подвоєння кількості процесорів вимагає подвоєння $n$, а шахова схема – лише трохи більшого за $\sqrt{2}$ збільшення $n$. Функцію ізоефективності можна знайти й експериментально: для кожного $p$ збільшувати розмір задачі, доки виміряна ефективність не досягне заданого значення (варіант 30 лабораторної роботи). Для аналізу вже виміряних даних зручна метрика Карпа–Флатта (тема 1).