Застосування FPGA для розрахунку деполімеризації мікротрубочки методом броунівської динаміки

Все готово, щоб розповісти Хабр аудиторії про застосування FPGA у сфері наукових високопродуктивних обчислень. І про те, як на цій задачі треба вдалося значно обскакати GPU (Nvidia K40) не тільки в метриці продуктивність на ват, але і просто з точки зору швидкості обчислення. Як FPGA платформи використовувався кристал Xilinx Virtex-7 2000t, підключений за PCIe до хост комп'ютера. Для створення апаратного обчислювального ядра використовувалася мова C++ (Vivado HLS).
Під катом текст нашої оригінальної статті. Там, як зазвичай буває, спочатку йде довгий опис навіщо це все треба і моделі, якщо немає бажання це читати, то можна переходити відразу до реалізації, а модель подивитися потім при необхідності. З іншого боку без хоча б побіжного ознайомлення з моделлю читач не зможе отримати враження про те, які складні обчислення можна реалізувати на FPGA )
Анотація
У даній роботі розглянуто апаратна реалізація розрахунку деполімеризації білкової мікротрубочки методом броунівської динаміки на кристалі програмованої логічної інтегральної схеми (FPGA) Xilinx Virtex-7 з використанням високорівневого транслятора з мови Сі Vivado HLS. Реалізація на FPGA порівнюється з паралельними реалізаціями цього ж алгоритму на багатоядерний процесор Intel Xeon і графічному процесорі Nvidia K40 за критеріями продуктивності та енергоефективності. Алгоритм працює на броунівських часи і тому вимагає великої кількості нормально розподілених випадкових чисел. Оригінальний послідовний код був оптимізований під многоядердную архітектуру з допомогою OpenMP, для графічного процесора — з допомогою OpenCL, а реалізація на FPGA була отримана за допомогою високорівневого транслятора Vivado HLS. В роботі показано, що реалізація на FPGA швидше CPU в 17 разів і швидше GPU в 11 разів. Що стосується енергоефективності (продуктивності на ват), FPGA була краще CPU в 160 разів і краще GPU в 80 разів. Прискорене на FPGA додаток було розроблено з допомогою SDK, що включає готовий проект FPGA, має PCI Express інтерфейс для зв'язку з хост-комп'ютером, і програмні бібліотеки для спілкування хост-додатки з FPGA прискорювачем. Від кінцевого розробника було необхідно лише розробити обчислювальне ядро алгоритму на мові Сі в середовищі Vivado HLS, і не потрібно спеціальних навичок FPGA розробки.
Введення
Високопродуктивні обчислення проводять на процесори (CPU), об'єднаних у кластери та/або мають апаратні прискорювачі – графічні процесори на відеокартах (GPU) або програмовані логічні інтегральні схеми (FPGA) [1]. Сучасний процесор сам по собі є відмінною платформою для високопродуктивних обчислень. До достоїнств CPU можна віднести многоядерную архітектуру з загальної когерентної кеш-пам'яттю, підтримку векторних інструкцій, високу частоту, а також величезний набір рограммных коштів, компіляторів та бібліотек, що забезпечує високу гнучкість програмування. Висока продуктивність платформи GPU ґрунтується на можливості запустити тисячі паралельних обчислювальних потоків на незалежних апаратних ядрах. Для GPU доступні добре зарекомендували себе засоби розробки (CUDA, OpenCL), що знижують поріг використання GPU платформи для прикладних обчислювальних завдань. Незважаючи на це, в останні роки FPGA все частіше стали використовуватися в якості платформи для прискорення завдань, в тому числі використовують речові обчислення [2]. FPGA володіють унікальною властивістю, різко відрізняє їх від CPU і GPU, а саме можливістю побудувати конвеєрну апаратну схему під конкретний обчислювальний алгоритм. Тому, незважаючи на значно меншу тактову частоту, на якій працюють FPGA (порівняно з CPU і GPU), на деяких алгоритмах на FPGA вдається домогтися більшої продуктивності [3]–[5]. З іншого боку, менша частота роботи означає менше енергоспоживання, і FPGA практично завжди більш ефективні, ніж CPU і GPU, якщо використовувати метрику «продуктивність на ват» [5].
Одним із класичних додатків, що вимагають високопродуктивних обчислень є метод молекулярної динаміки, що використовується для розрахунку руху систем атомів або молекул. У рамках цього методу взаємодії між атомами і молекулами описуються в рамках законів Ньютонівської механіки за допомогою потенціалів взаємодії. Розрахунок сил взаємодії проводиться ітеративне і являє істотну обчислювальну складність, враховуючи велику кількість атомів/молекул в системі і велику кількість розрахункових ітерацій. Прискоренню розрахунків молекулярної динаміки було приділено багато уваги в літературі в різних системах: суперкомп'ютерах [6], кластери [7], спеціалізованих під молекулярно-динамічні розрахунки машинах [8]–[10], машинах c прискорювачами на основі GPU [11] і FPGA [12]–[17]. Було продемонстровано, що FPGA може бути конкурентною альтернативою в якості апаратного прискорювача для молекулярно-динамічних обчислень у багатьох випадках, однак на сьогоднішній день не існує консенсусу про те, для яких саме завдань і алгоритмів вигідніше застосовувати платформу FPGA. В даній роботі ми розглядаємо важливий окремий випадок молекулярної динаміки – броуновскую динаміку. Основна особливість методу броунівської динаміки порівняно з молекулярною динамікою полягає в тому, що, молекулярна система моделюється більш грубо, тобто в якості елементарних об'єктів моделювання виступають не окремі атоми, а більш великі частинки, такі як окремі домени макромолекул або цілі макромолекули. Молекули розчинника й інші малі молекули в явному вигляді не моделюються, а їх ефекти враховуються у вигляді випадкової сили. Таким чином вдається значно знизити розмірність системи, що дозволяє збільшити інтервал часу, покриваються модельними розрахунками на порядки. Нам невідомі описані в літературі спроби дослідити ефективність FPGA порівняно з альтернативними платформами для прискорення завдань броунівської динаміки. Тому ми зробили дослідження даного питання на прикладі задачі моделювання деполімеризації мікротрубочки методом броунівської динаміки. Мікротрубочки – це трубки діаметром близько 25 нм і довжиною від декількох десятків нанометрів до десятків мікрон, що складаються білка тубуліну і входять до складу внутрішнього скелета живих клітин. Ключовою особливістю мікротрубочок є їх динамічна нестабільність, тобто можливість спонтанно перемикатися між фазами полімеризації і деполімеризації [18]. Це поведінка важливо насамперед для захоплення і переміщення хромосом мікротрубочками під час клітинного поділу. Крім того, мікротрубочки відіграють важливу роль у внутрішньоклітинному транспорті, рух війок і джгутиків і підтримці форми клітини [19].
Механізми, що лежать в основі роботи мікротрубочок, досліджуються вже декілька десятків років, але лише недавно розвиток обчислювальних технологій дало змогу описувати поведінку мікротрубочок на молекулярному рівні. Найбільш детальна молекулярна модель динаміки мікротрубочок, створена нещодавно нашою групою на основі методу броунівської динаміки, була реалізована базі CPU і дозволяла розраховувати часи полімеризації/деполімеризації мікротрубочок порядку декількох секунд [20]. Це пролило світло на ряд важливих аспектів динаміки мікротрубочок, однак, тим не менш, багато ключові експериментально спостережувані явища залишилися за рамками теоретичного опису, оскільки вони відбуваються в микротрубочках на часи десятків і навіть сотень секунд [21]. Таким чином, для прямого порівняння теорії та експерименту критично важливо досягти прискорення розрахунків динаміки мікротрубочок хоча б на порядок.
У даній роботі ми досліджуємо можливість прискорення розрахунків броунівської динаміки мікротрубочки на FPGA і порівнюємо результати, отримані при реалізації одного і того ж алгоритму динаміки мікротрубочок на трьох різних платформах, за критеріями продуктивності та енергоефективності.
Математична модель
Загальні відомості про структуру мікротрубочки
Структурно микротрубочка являє собою циліндр, що складається з 13 ланцюжків – протофиламентов.

На малюнку ліворуч – схема моделі мікротрубочки. Сірим показано субодиниці тубуліну, чорними крапками – центри взаємодії між ними. Праворуч — вид енергетичних потенціалів взаємодії між тубулинами.
Кожен протофиламент побудований з димерів білка тубуліну. Сусідні протофиламенты пов'язані один з одним бічними зв'язками і зрушені відносно один одного на відстань 3/13 довжини одного мономеру, так що микротрубочка має спиральность. При полімеризації димери тубуліну приєднуються до кінців протофиламентов, причому протофиламенты мікротрубочки прагнуть брати пряму конформацію. При деполімеризації бічні зв'язки між протофиламентами на кінці мікротрубочки розриваються, і протофиламенты закручуються назовні. При цьому від них випадковим чином відриваються олігомери тубуліну.
Моделювання деполімеризації мікротрубочки методом броунівської динаміки
Використовувана тут молекулярна модель мікротрубочки була вперше представлена в статті [20]. Оскільки завданням цього дослідження було порівняння продуктивності різних обчислювальних платформ, ми обмежилися моделюванням тільки деполімеризації мікротрубочки.
Коротко, микротрубочка моделювалася як набір сферичних частинок, що представляють собою мономери тубуліну. Мономери могли рухатися лише у відповідній їм радіальної площини, тобто площини, що проходить через вісь мікротрубочки і відповідний протофиламент. Таким чином, положення та орієнтація кожного мономеру повністю визначалися трьома координатами: двома декартовими координатами центру мономеру і кутом орієнтації. Кожен мономер мав чотири центру взаємодії на своїй поверхні: два центру бічного взаємодії і два центру поздовжнього взаємодії. Енергія тубулін-тубулинового взаємодії залежала від відстані r між сайтами взаємодії на поверхні сусідніх субодиниць і від кута нахилу між сусідніми мономерами тубуліну в протофиламенте. Бічні і поздовжні взаємодії між диммерами тубуліну визначалися потенціалом, має наступний вигляд:

де A і b визначали глибину потенційної ями і висоту енергетичного бар'єру, r0 і d – параметри, які визначають ширину потенційної ями і форму потенціалу в цілому. Параметр A приймав різні значення для бічних і поздовжніх зв'язків, так що бічні взаємодії були слабші поздовжніх, всі інші параметри збігалися для обох типів зв'язків (повний список параметрів і їх значень представлений в Таблиці 1 [20]). Поздовжні взаємодії всередині димеру моделювалися як нерозривні пружини з квадратичним енергетичним потенціалом u_r:

де k — жорсткість зв'язку тубулін-тубулинового взаємодії. Енергія вигину g(χ) пов'язана з поворотом мономерів один щодо одного і також описувалася квадратичних нерозривному функцією:

де χ — кут між сусідніми мономерами тубуліну в протофиламенте, χ 0 — рівноважний кут між двома мономерами, B — изгибная жорсткість. Повна енергія мікротрубочки записувалася наступним чином:

де n — номер протофиламента, i — номер мономеру в n-му протофиламенте, Kn — кількість субодиниць тубуліну в n-му протофиламенте, бічного взаємодії v_k_lateral — енергія бічного взаємодії між мономерами, v_k_longitudinal — енергія поздовжнього взаємодії між диммерами.
Еволюція системи розраховувалося за допомогою методу Броунівської динаміки [22]. Початковою конфігурацією мікротрубочки була коротка «затравка», містить 12 мономерів тубуліну в кожному протофиламенте. Ми розглядали тільки деполімеризацію МТ і моделювали всі тубулины з рівноважним кутом χ 0 = 0.2 радий. Координати всіх мономерів системи на i-ої ітерації виражалися наступним чином:

де dt — крок по часу, U total виражається через (4), k_B — постійна Больцмана, T — температура, N(0,1) – випадкове число з нормального розподілу, створений за допомогою алгоритму вихор Мерсенна [23]. γ q γ τ — в'язкісні коефіцієнти опору для зсуву і повороту відповідно, розраховані для сфер радіуса r = 2 нм.
Похідна повної енергії за незалежними координатами q{k,n} виражалася через бічну, поздовжню складові енергії взаємодії між сусідніми диммерами і всередині димера, а також енергію вигину:

Для прискорення розрахунків були використані аналітичні вирази для
всіх градієнтів енергії:

Слід зазначити, що розмір цієї задачі порівняно малий. Ми розглядали тільки 12 шарів мономерів, що дає повне число частинок дорівнює 156. Однак, це аніскільки не зменшує значущість обчислень, оскільки в реальних розрахунках достатньо обчислювати положення крайніх декількох (близько 10) слів мономерів, т. к. при зростанні мікротрубочки далекі від кінця мікротрубочки молекули тубуліну утворюють стійку циліндричну конфігурацію, і брати їх в розрахунок немає сенсу.
Псевдокод алгоритму розрахунку
Алгоритм є ітеративним по часу з кроком 0.2 нс. Існують масив тривимірних координат молекул, а також масиви сил поперечного (латерального) і поздовжнього (лонгитудального) взаємодій. На кожній ітерації по часу послідовно виконуються два вкладених циклу з молекул, у першому виробляється обчислення сил взаємодії по відомим координатам, у другому – оновлюються самі координати. У циклі обчислення сил взаємодії необхідно прочитати координати трьох молекул, одна центральна і дві сусідні («ліва» і «верхня», див. Рис.2), а результатом обчислення буде сила поперечного взаємодії між центральною і лівої молекулами і сила поздовжнього взаємодії між центральною і верхній.

Рис. 2. Схема розташування взаємодіючих субодиниць в моделі мікротрубочки
У підсумку після цього циклу виявляються обчисленими всі сили взаємодії між усіма молекулами. У циклі оновлення координат за відомим силам обчислюються зміни координат, а також беруться до уваги випадкові добавки для обліку Броунівського руху. Таким чином, псевдокод алгоритму можна записати наступним чином.
Вхід: масив координат молекул M = {x, y, teta}. Граничні умови на сили
взаємодії.
Вихід: масив координат M після K кроків по часу

for t in {0.. K-1} do
for i in {0.. 13} // кількість протофиламентов
for j in {0.. 12} do // кількість шарів молекул
Mc <- M[i,j]
Ml <- M[i+1,j]
Mu <- M[i,j+1]
// за формулами (7, 8, 9, 10)
F_lat[i,j] <- calc_calteral(Mc, Ml)
F_long[i,j] <- calc_long(Mc, Mu)
end for
end for
for i in {0.. 13}
for j in {0.. 12} do
// за формулами (5)
M[i,j] <- update_coords(F_lat[i,j], F_long[i,j])
end for
end for
end for

Програмна реалізація на CPU і GPU
Реалізація на CPU
Була зроблена спроба максимально розпаралелити код на CPU Intel Xeon E5-2660 2.20 GHz під управлінням ОС Ubuntu 12.04 з допомогою бібліотеки OpenMP. Паралельна секція починалася до циклу за часом. Цикли розрахунку сил взаємодії та оновлення координат були распараллелены з допомогою директиви omp for schedule(static), між циклами була вставлена бар'єрна синхронізація. Масиви, що містять сили взаємодії та координати молекул, були оголошені як private для кожного потоку.
При реалізації розрахунків на CPU було виявлено, що розмір завдання не дозволяв її ефективно розпаралелювати. Залежність часу виконання однієї ітерації від кількості паралельних потоків була немонотонна. Мінімальний час розрахунку однієї ітерації по часу було отримано при використанні всього 2 потоків (ядер CPU). Пояснюється це тим, що зі збільшенням кількості потоків зростає час на копіювання даних між потоками і їх синхронізацію. При цьому розмір завдання дуже малий, щоб виграш від збільшення кількості ядер перевищив ці накладні витрати. При цьому експерименти показали, що завдання слабо масштабується при збільшенні розміру (weak scaling), тобто при дновременном збільшення розміру задачі і числа паралельних потоків час обчислення залишалося приблизно однаковим. У підсумку кращим результатом на даному CPU було 22 мкс на одну ітерацію з часу при використанні двох ядер CPU. Код не був векторизован з-за складності обчислень сил взаємодії.
Реалізація на GPU
Ми запускали OpenCL реалізацію на граф процесорі Nvidia Tesla K40. Цикли, обчислюють сили взаємодії та оновлення координат були распараллелены, головний цикл за часом був ітеративним. Були реалізовані два варіанти – з однією і кількома робочими групами (work groups). В першому випадку було виділено по одному робочого потоку (work item) на кожну молекулу. В кожному потоці був цикл за часом, у якому обчислювалися сили і координати молекули потоку. При цьому застосовувалася бар'єрна синхронізація після обчислення сил і після оновлення координат. У цьому випадку участь хоста не вимагалося для обчислень, він тільки займався управлінням та запуском ядер.
У другому випадку були два типи потоків, в одному просто обчислювалися сили для однієї молекули, у другому – оновлювалися координати. Головний цикл за часом був на хості, який керував запуском і синхронізацією ядер на кожній ітерації циклу.
Найбільша продуктивність була отримана в розрахунках з однією групою потоків і бар'єрної синхронізації між ними. Без використання генераторів псевдовипадкових чисел одна ітерація обчислювалася протягом 5 мкс, якщо використовувати один генератор чисел на всі потоки, то час роботи зростало до 9 мкс, а при максимальному заповненні загальної пам'яті (shared memory) вдавалося включити 7 незалежних генераторів, при цьому час обчислення однієї ітерації по часу склало 14 мкс, що було в 1.57 рази швидше реалізації на CPU.
Завантаженість ядер GPU склала 7% від одного мультипроцесора (SM), при цьому загальна пам'ять, де розміщувалися масиви сил, координат і буфери даних генераторів псевдовипадкових чисел, була заповнена на 100%. Тобто з одного боку розмір задачі був явно замалий для повного завантаження GPU, з іншого боку, при збільшенні розміру задачі довелося б використовувати глобальну DDR пам'ять, що могло б призвести до обмеження зростання продуктивності.
Реалізація на FPGA
Опис платформи
Обчислення на FPGA проводилися на платформі RB-8V7 виробництва фірми НВО «Зростання». Вона являє собою 1U блок для установки у стійку. Блок складається з 8 кристалів Xilinx FPGA Virtex-7 2000T. Кожна FPGA має 1 GB зовнішньої пам'яті DDR3 і PCI Express 2.0 x4 інтерфейс до внутрішнього PCIe комутатора. Блок має два інтерфейсу PCIe 3.0 x4 до хост-комп'ютеру через оптичні кабелі, які повинні бути з'єднані зі спеціальним адаптером, встановленим в хост-комп'ютер.
як хост-комп'ютера був використаний сервер з CPU Intel Xeon E5-2660 2.20 GHz, що працює під управлінням ОС Ubuntu 12.04 LTS – такий же як і для обчислень просто на CPU з допомогою OpenMP. Програмне забезпечення, що працює на CPU хост-комп'ютера «бачить» блок RB-8V7 як 8 незалежних FPGA пристроїв, підключених по шині PCI Express. Далі буде описуватися взаємодія CPU тільки з одного FPGA XC7V72000T, при цьому система дозволяє використовувати FPGA незалежно і паралельно.
Прискорене за допомогою FPGA додаток було розроблено з допомогою SDK з наступною моделлю. На CPU хост-комп'ютера (далі просто CPU) працює основна програма, яка використовує прискорювач FPGA для найбільш обчислювально ємних процедур. CPU передає дані в прискорювач і назад через зовнішню пам'ять DDR, підключену до FPGA, а також керує роботою обчислювального ядра в FPGA. Обчислювальне ядро складається заздалегідь на мові C/C++, перевіряється і транслюється в RTL код за допомогою засобу Vivado HLS. RTL код обчислювального ядра вставляється в основний FPGA проект, в якому вже реалізована необхідна логіка управління і передачі даних, що включає PCI Express ядро, DDR контролер і шину на кристалі (Рис. 3). Основний FPGA проект іноді називають Board Support Package (BSP), він розробляється виробником обладнання, і від користувача не потрібно його модифікації. Обчислювальне ядро HLS після запуску сама звертається до DDR пам'ять, зчитує звідти вхідний буфер даних для обробки і записує туди ж результат обчислень. На рівні мови C++ звернення в пам'ять відбувається через аргумент функції верхнього рівня обчислювального ядра типу покажчик.

Рис. 3. Блок-схема проекту FPGA. Синім і жовтим кольорами позначені блоки, що входять в BSP. Зеленим позначені обчислювальні HLS ядра. Також позначено розбиття ядер проекту на блоки (pBlocks) для накладання пространсвенных обмежень при трасуванні.
Для створення прискореного програми була розроблена методологія, що складається з декількох кроків. По-перше, оригінальний послідовний код компилировался в середовищі Vivado HLS, і перевірялося, що скомпільований таких код чином не змінює вихідних даних опорного послідовного коду. По-друге, з цього коду виділялася основна обчислювальна та підходяща для прискорення частина; ця частина відокремлювалася від основного коду за допомогою функції-обгортки. Після чого створювалося дві копії такої функції і логіка перевірки на відповідність результатів обох частин. Перша копія була опорною реалізацією алгоритму в Vivado HLS, а друга була оптимізована для трансляції в RTL код. Оптимізації включали у собі переписування коду, такі як використання статичних масивів замість динамічних, використання спеціальних функцій вводу/виводу в HLS ядро, методи економії пам'яті та перевикористання результатів обчислень. Після кожної зміни результат функції порівнювався з результатом опорної реалізації. Іншим методом оптимізації було використання спеціальних директив Vivado HLS, не змінюють логічне поведінку, але впливають на кінцеву продуктивність RTL коду. На цій стадії слід залишатися до тих пір, поки не будуть отримані задовільні попередні результати трансляції C RTL, такі як продуктивність схеми та займані ресурси.
Наступна стадія – це імплементація розробленого обчислювального ядра в системі Vivado поза контекстом основного проекту. Тут завдання домогтися відсутності тимчасових помилок вже розведеного дизайну всередині розробленого обчислювального ядра. Якщо на цьому етапі спостерігаються тимчасові помилки, то можна застосовувати інші параметри імплементації, або повертатися на попередню стадію і намагатися змінити С++ код або використовувати інші директиви.
На наступній стадії необхідно імплементувати обчислювальне ядро вже разом з основним проектом і його часовими та просторовими обмеженнями. На цій стадії також необхідно домогтися відсутності тимчасових помилок. Якщо вони спостерігаються, то можна змінити частоту роботи обчислювальної схеми, накласти інші просторові обмеження на розміщення схеми на кристалі, чи знову зайнятися
зміною C++ коду і/або використовувати інші директиви.
Остання стадія розробки – це перевірка на відповідність результатів, отриманих на реальному запуску в залозі і з допомогою опорної моделі на CPU. Проходить вона на невеликому проміжку часу, при цьому вважається, що на більш тривалих запусках (коли порівняти з CPU вже проблематично) FPGA рішення видає правильні результати.
Робота в середовищі Vivado HLS
У роботі використано два Vivado HLS ядра (Рис. 3): основне ядро, що реалізує алгоритм молекулярної динаміки мікротрубочок (MT ядро), і ядро для генерації псевдовипадкових чисел (RAND ядро). Нам довелося розділити алгоритм на два обчислювальних ядра з наступної причини. Кристал FPGA Virtex-7 2000T – це найбільший кристал FPGA сімейства Virtex-7 на ринку. Він насправді складається з чотирьох кристалів кремнію, сполучених на підкладці безліччю з'єднань і об'єднаних в один корпус мікросхеми. По термінології Xilinx кожен такий кристал називається SLR (Super Logic Region). При використанні таких великих FPGA завжди виникають проблеми з ланцюгами, що перетинають кордони SLR. Xilinx рекомендує вставляти регістри на такі ланцюги з обох сторін кордону SLR.
Повне HLS ядро, що включає і MT і RAND ядра, вимагали апаратних ресурсів більше, ніж було доступно в одному SLR, тому були ланцюги, які перетинали кордон незалежних кристалів кремнію. На стадії трансляції з мови C++ в RTL Vivado HLS нічого «не знає» про те, які ланцюга згодом будуть перетинати кордон, і тому не може заздалегідь вставити додаткові регістри синхронізації. Тому ми прийняли рішення розділити ядра на два, просторово обмежити їх в різні SLR і вставити регістри синхронізації на інтерфейсні ланцюги між ядрами на рівні RTL.
Ядро MT
Даний алгоритм дуже добре підходить для реалізації на FPGA, тому для порівняно невеликої кількості даних з пам'яті (координати двох молекул) необхідно обчислити складну функцію сил взаємодії і вдається побудувати довгий обчислювальний конвеєр.

Рис. 4. Блок-схема апаратної обчислювальної процедури ядра MT. Зеленим позначені апаратні блоки пам'яті для зберігання координат молекул. Позначені обчислювальні конвеєри сил і оновлення координат, а також блок Save Regs для зберігання проміжних результатів обчислень. Псевдовипадкові числа надходять на конвеєр оновлення координат з іншого HLS ядра.
Кожна молекула, тобто мономер тубуліну, взаємодіє тільки з чотирма своїми сусідами (Рис. 2). На кожній ітерації по часу треба спочатку обчислити сили взаємодії, а потім оновити координати молекул. Функції обчислення сил взаємодії включають в себе безліч арифметичних, експоненціальних і тригонометричних операторів. Нашою першою задачею було синтезувати конвеєр для цих функцій. Робочим типом даних був дійсний тип float. Vivado HLS синтезувала такі функції у вигляді конвеєрів, що працюють на частоті 200 МГц, з латентністю близько 130 тактів. При цьому конвеєри були однотактовые (або, як кажуть, з інтервалом ініціалізації дорівнює 1), що означає, що на вхід вони могли приймати координати нових молекул кожен такт, а потім після початкової затримки (латентність) – видавати оновлені значення сил також кожен такт. Вихідні сили взаємодії використовувалися для оновлення координат, що теж було конвейеризовано. Для оновлення кожної координати кожної молекули були необхідно незалежні псевдовипадкові нормально розподілені числа, одержувані з іншого HLS ядра. Якщо взяти три молекули («поточна», «ліву» та «верхню») вийшло можливим об'єднати конвеєри обчислення сил і оновлення координат в один конвеєр, що реалізує всі обчислення для однієї молекули. Такий конвеєр мав латентність рівну 191 такт (Рис. 4).
Алгоритм проходить по всіх молекул в циклі. На кожній ітерації циклу необхідно мати координати трьох молекул: одна молекула розглядається як «поточна», також є «ліва» і «права» молекули. Відповідно розраховуються сили взаємодії між цими трьома молекулами. Далі при оновленні координат поточної молекули ліва і верхня компоненти сил взаємодії бралися з розрахунку на поточній ітерації, а нижня і права компоненти бралися або з граничних умов, або з попередніх ітерацій з локального регістрового файлу Save Regs (Рис. 4).
Кількість молекул N у системі було невеликим (13 протофиламентов х 12 молекул = 156 молекул). На кожну молекулу потрібно 12 байт. Схема використовувала два масиву координат m1 і m2, загальним обсягом менше 4 КБ, відповідно ці дані легко поміщалися у внутрішню пам'ять FPGA – BRAM, реалізовану всередині HLS ядра. Схема була влаштована таким чином, що на парних ітераціях за часом координати зчитувалися з масиву m1 (і записувалися в m2), а на непарних – навпаки. З точки зору алгоритму можна було читати і писати в один масив координат, але Vivado HLS не могла створити схему, здатну на одному і тому ж такті читати і писати один і той же апаратний масив, що потрібно для роботи однотактового конвеєра. Тому було прийнято рішення подвоїти кількість незалежних блоків пам'яті.

Рис. 5. Схема конвеєрного розрахунку взаємодій тубулінів у микротрубочке.
Виявилося можливим реалізувати три повних паралельних конвеєра, здатних оновлювати координати трьох молекул кожен такт (Рис. 5). Тоді, щоб уникнути простою конвеєрів необхідно було збільшити пропускну здатність до локальної пам'яті і читати координати семи молекул кожен такт. Це проблема легко зважилася, практично не змінюючи вихідний C++ код, а лише за рахунок використання спеціальної директиви,
фізично розбиває вихідний масив даних по чотирьом незалежним апаратним
блоків пам'яті. Т. к. пам'ять BRAM в FPGA є двупортовой, то з чотирьох блоків пам'яті можна прочитати 8 значень за такт. Але, так як три конвеєра вимагають координати 7ми молекул за такт (див. рис 5), це вирішило проблему.
#pragma HLS DATA_PACK variable=m1, m2
#pragma HLS ARRAY_PARTITION variable=m1, m2 cyclic factor=4 dim=2






Період L II BRAM DSP FF LUT Утилізація 5 нс 191 такт 1 такт 52 498 282550 331027 Абсолютна 2 % 23 % 11 % 27 % Відносна
Табл. 1: Продуктивність та утилізація схеми HLS з трьома повними конвеєрами
В табл. 1 наводиться утилізація схеми HLS (тобто кількість споживаних їй апаратних ресурсів FPGA, в абсолютних і відносних одиницях для кристала Virtex-7 2000T) та її продуктивність. L – це затримка або латентність схеми, тобто кількість тактів між подачею в конвеєр перших вхідних даних та одержанням перших вихідних даних, II – це інтервал ініціалізації (або пропускна здатність) конвеєра, що означає через скільки тактів на вхід конвеєра можна подавати такі дані.
Утилізація наводиться як в абсолютних величинах (скільки потрібно тригерів FF або таблиць LUT) для реалізації схеми, так і у відносних до повного кількістю даного ресурсу в кристалі. Як видно з табл. 1 латентність L повного конвеєра була дорівнює 191 такту, при цьому кожен конвеєр повинен був обробити третю частину молекул, що дає теоретичну оцінку часу обчислення однієї ітерації дорівнює T(FPGA) = (L+N/3)*5нс=1.2 мкс
З табл. 1 також видно, що в кристалі залишилося ще багато невикористаної логіки, але далі збільшувати кількість паралельних конвеєрів непрактично. Буде зменшуватися тільки другий доданок, а початкова затримка все одно буде давати значний внесок під час роботи. При цьому збільшення кількість логіки ускладнить розміщення і трасування схеми на наступних стадіях розробки проекту в Vivado.
Ядро RAND
Як було зазначено, алгоритм враховує Броунівський рух молекул, одним з методів розрахунку якого є збільшення нормальної випадкової добавки до зміни координат на кожній ітерації по часу. Необхідно дуже багато нормально розподілених випадкових чисел, на кожну ітерацію – по N3 чисел, що дає потік 420 10 6 чисел/с. Такий потік не може бути завантажений з хоста, тому його необхідно генерувати всередині FPGA «на льоту». Для цього, як і в опорному коді для CPU, був обраний генератор вихор Мерсенна, дає рівномірно розподілені псевдослучаные числа. Далі до них застосовувалося перетворення Боксу-Мюллера і на виході виходили нормально розподілені послідовності. Вихідний відкритий код вихор Мерсенна був модифікований для отримання апаратного конвеєра з інтервалом ініціалізації 1 такт. Алгоритм вимагає 9 нормальних цифр кожен такт, тому ядро RAND включало в себе 10 незалежних генераторів вихор Мерсенна, т. к. перетворення Боксу-Мюллера вимагає два рівномірно розподілених чисел для отримання 2х нормально розподілених. В табл. 2 наводиться утилізація ядра RAND.





BRAM DSP FF LUT Утилізація 30 41 48395 64880 Абсолютна 1.2 % 9 % 0.1 % 5.3 % Відносна
Табл. 2: Утилізація ядра RAND
Видно, що таке ядро вимагає значну частину DSP ресурсів кристала, і це ядро було б складно розмістити в одному SLR з ядром MT, т. до. сума утилизаций двох ядер хоча б по DSP ресурсу 31% більше, ніж може вмістити один SLR (25 %).
Більш детально про те, як ми вибирали генератор псеводслучайных чисел і синтезували його в Vivado HLS можна почитати в статті мого колеги https://habrahabr.ru/post/266897/
Створення битстрима
Після інтеграції обчислювальних ядер в Vivado на проект були накладені просторові обмеження на розміщення IP блоків. Використовувана FPGA Virtex-7 2000T має 4 незалежних кристала кремнію (SLR0, SLR1, SLR2, SLR3). Було показано, що ядро MT не вміщувалося в один SLR, тому було вирішено створити два регіону розміщення (pBlock): pBloch_hls для розміщення тільки MT ядра і pBlock_base для розміщення інших ядер проекту (Рис. 3). Регіон розміщення pBlock_hls включав в себе SLR0 і SLR1, pBock_base – SLR2. Такий підхід дозволив розмістити логіку оптимальним чином, вставити регістри синхронізації на інтерфейси, що перетинають регіони розміщення (а значить і SLR) і добитися позитивних времянных результатів після трасування проекту.
Формат хабра дозволяє підключати багато картинок, тому ось ще один малюнок, як в результаті пройшла імплементація проекту.

Червоним на картинці підсвічені елементи Board Support Package (PCIe core, DDR3 Interface, Internal AXI Bus), блакитним — MT HLS ядро, а пурпуровим — елементи ядра RAND.
Результати
Продуктивність
Результати роботи всіх трьох реалізацій (CPU, GPU і FPGA) були логічно верифіковані відносного оригінального коду і визнані заможними. Порівняння продуктивності вироблялося виміром часу роботи програм на 10 7 ітерацій алгоритму та обчисленням часу, потрібного для розрахунку однієї ітерації. При цьому продуктивність GPU і FPGA платформ брали до уваги час передачі даних між хост-комп'ютером і прискорювачем.
Для оцінки продуктивності зазвичай використовується метрика операцій в секунду. Для даного алгоритму нам виявилося складним обчислити точне значення речових операцій, тому ми просто порівнюємо часи роботи алгоритму для обчислення однієї ітерації, визначаючи продуктивність CPU платформи дорівнює 1. Результати порівняння наведено в табл. 3, у другому стовпці якої наводяться часи обчислення однієї ітерації алгоритму в мікросекундах, а в третьому – відносна продуктивність платформ.






Платформа Час, мкс Продуктивність CPU 22 1 GPU 14 1.6 FPGA 1.3 17
Табл. 3: Порівняння продуктивності трьох платформ
З таблиці видно, що реалізація на GPU швидше CPU всього в 1.6 рази, в той час як FPGA швидше CPU в 17 разів. Це означає, що FPGA швидше GPU в 11 разів. Отримане експериментально час роботи FPGA одно 1.3 мкс на ітерацію більше розрахункового часу в 1.2 мкс через урахування накладних витрат на передачу даних по шині PCI Express.
Енергоефективність
Для вимірювання енергоспоживання ми використовували такі кошти. Для CPU платформи – утиліту Intel Power Gadget. Для GPU платформи — утиліту Nvidia-smi. Для FPGA – спеціальні програмно-апаратні засоби, включені до складу блоку RB-8V7. У всіх випадках замірялася різниця в споживанні всього чіпа до запуску завдання і під час обчислень. Результати наведені в таблиці 4.






Платформа Потужність, Вт Ex Ex_rel CPU 89.6 0.011 1 GPU 67 0.023 2 FPGA 9.6 1.77 160
Табл. 4: Порівняння енергоспоживання обчислювальних платформ
У другому стовпці таблиці наведено потужність, що виділяється при розрахунку на різних платформах. У третьому стовпці наводяться значення абсолютної енергоефективності (продуктивності на Вт) для даної задачі, що визначається за формулою:

У четвертому стовпці наводяться значення відносної енергоефективності різних платформ для даної задачі, що обчислюється за формулою:

Для обох формул, x = {CPU, GPU, FPGA}.
Видно, що у FPGA є велика перевага в енергоефективності перед іншими платформами, що може зіграти роль в середньо та довгостроковій перспективі використання FPGA прискорювачів в датацентрах при оплаті рахунків за електроенергію. Досягається це в першу чергу за рахунок того, що FPGA працюють на порядок меншою частотою.
Обговорення
В раніше опублікованих роботах технологія FPGA неодноразово застосовувалася до вирішення завдань молекулярної динаміки [12]–[17]. Дослідникам з лабораторії CAAD Бостонського Університету вдалося розробити ефективне ядро для розрахунку короткодействующих міжмолекулярних сил, яке було реалізовано на платі ProcStar-III (виробництво фірми Gidel), з встановленим кристалом FPGA Altera Stratix-III SE260. Плата мала PCI Express інтерфейс до хост-комп'ютеру. Було показано, що розроблене прискорене рішення було в 26 разів швидше чистої реалізації на CPU на бенчмарку Apoal. У роботі [24] автори перенесли частину пакета для розрахунку молекулярної динаміки LAMMPS на FPGA. Експрес частина включала в себе обчислення дальнодействующих взаємодій. Розроблене програмне ядро складалося з чотирьох однакових незалежних конвеєрів, що працюють паралельно. Задача була виконана на суперкомп'ютері Maxwell, кожен вузол якого складається з одного процесора Intel Xeon і двох кристалів Xilinx FPGA Virtex-4 [25]. Автори заявили, що розроблене прискорене рішення легко масштабировалось на безліч вузлів суперкомп'ютера Maxwell. З аналізу продуктивності тільки прискорювача випливало, що на двох вузлах комп'ютера можна було отримати прискорення в 13 разів у порівнянні з чисто програмним рішенням. Однак повне час роботи гібридного рішення було гірше чисто програмного із-за того, що час на пересилання даних між CPU і зовнішньою пам'яттю SDRAM, підключеної до FPGA займало 96% часу роботи алгоритму. Але в роботі стверджується, що якщо поліпшити інтерфейс передачі даних, то можна отримати повний виграш у швидкості у 8-9 разів.
У цій роботі ми застосували FPGA до розрахунку руху ансамблю білкових молекул методом броунівської динаміки. Наша програмно-апаратна реалізація алгоритму деполімеризації мікротрубочки показала, що продуктивність FPGA при розрахунку однієї траєкторії мікротрубочки в 17 разів перевищувала продуктивність CPU і в 11 разів продуктивність GPU. Отримане прискорення при розрахунку деполімеризації мікротрубочки методом броунівської динаміки дозволяє здійснити розрахунок на часах порядку декількох десятків і навіть сотень секунд. Це дозволить передбачати поведінку реальних мікротрубочок на експериментально доступних часи і проаналізувати механізми динамічної нестабільності мікротрубочок, що буде предметом майбутньої роботи в даному напрямку.
Отриманий виграш на завданні броунівської динаміки дозволяє говорити про перспективність застосування FPGA для вирішення даного типу завдань. Наскільки нам відомо, це перша спроба порівняти продуктивність і енергоефективність різних типів апаратних прискорювачів на даному алгоритмі.
Довгий час головною проблемою використання FPGA була відсутність високорівневих засобів програмування. Традиційні мови опису апаратури завжди вимагають значного часу для реалізації алгоритму, в той час як перші високорівневі транслятори [26] генерували RTL код низької якості. Однак, кілька років тому компанії Altera і Xilinx стали приділяти значні ресурси цієї проблеми і випустили на ринок свої високорівневі засоби програмування (Altera SDK for OpenCL і Xilinx Vivado HLS). Дані транслятори генерують набагато більш ефективний код і дозволяють прикладному програмісту використовувати мови C/C++ (Xilinx) і OpenCL (Altera) для створення якісних апаратних обчислювальних схем. Останнім часом з'явилося безліч робіт, у яких використовувалися кошти високорівневого синтезу для розробки FPGA прискорювачів [27]–[29]. Наприклад, в роботі [28] за допомогою засобу Vivado HLS реалізований алгоритм оптичного потоку на платформі Xilinx Zynq-7000. Розроблена система мала продуктивність, порівнянну з реалізацією на CPU, при цьому споживання енергії було в 7 разів менше. Автори особливо підкреслювали, що використання коштів HLS порівняно з традиційними RTL мовами значно скоротило термін розробки. Використання засобу Vivado HLS у ході виконання цієї роботи також дозволило значно скоротити час і трудомісткість розробки і залучати до програмування розробників, які не володіють спеціальними навичками роботи з FPGA. Все це дозволяє говорити про FPGA як про що відбулася платформі для високопродуктивних обчислень в області молекулярної і броунівської динаміки.
Фуф, спасибі за увагу! Вдячності і посилання на літературу можна знайти на оригінальному pdf-документі, доступному на сайті Праць ІСП РАН
» Програмний код лежить тут
Джерело: Хабрахабр

0 коментарів

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