Применение технологий параллельного программирования для систем с общей памятью при решении гиперболических систем уравнений Текст научной статьи по специальности «Математика»

УДК 519.633

А. М. Иванов, Н. И. Хохлов

Московский физико-технический институт (государственный университет)

Применение технологий параллельного программирования для систем с общей памятью при решении гиперболических систем уравнений

Рассматривается применение технологий параллельного программирования OpenMP и POSIX Threads для решения гиперболических систем уравнений. Данные технологии предназначены для систем с общей памятью. Кроме этого, рассматривается увеличение производительности при использовании векторных инструкциий процессора. Решается задача распространения динамических волновых возмущений в геологической среде в упругой постановке в двумерном случае. Для численного решения используется сеточно-характеристический метод. Исследуется влияние привязки потоков к определенным ядрам процессора в NUMA-системах.

Ключевые слова: математическое моделирование, сеточно-характеристический метод, параллельное программирование, общая память, гиперболические уравнения, векторизация.

A.M. Ivanov1, N.I. Khokhlov1

1 Moscow Institute of Physics and Technology (State University)

Applying parallel programming technologies for shared memory systems to solving hyperbolic PDE

We consider the application of OpenMP and POSIX Threads parallel programming technology to the solution of hyperbolic equations. These technologies are designed for shared memory systems. Also, an increase in productivity due to using vector instructions for processor is considered. The problem is a two-dimensional elastic wave propagation. For the numerical solution the grid characteristic method is used. The effect of the binding threads on specific processor cores in NUMA systems is considered.

Key words: mathematical modelling, grid characteristic method, parallel programming, shared memeory, hyperbolic PDE, vectorization.

1. Введение

За последние годы все чаще используются технологии параллельного программирования для решения разного рода задач. Один из классов ресурсоемких задач - это задачи вычислительной геофизики, сейсмики. Задачи такого рода требуют для своего решения большого количества вычислительных ресурсов, поэтому немаловажен вопрос о написании высокопроизводительных алгоритмов для их решения.

С ростом возможностей современных вычислительных систем возникает вопрос об их эффективном использовании. Так как большое количество алгоритмов было написано довольно давно, не все из них могут использовать параллельные системы с достаточной эффективностью и требуют доработки. Так, в работе [1] последовательный код для моделирования эффектов, происходящих в грунте при землетрясении, преобразуется в параллельный. Для этого используются технологии ОрепМР и ИРР. Авторы сравнивают ускорение и эффективность при выполнении на нескольких процессорах.

В работе [2] решается задача сейсморазведки. Программная реализация задействует инструкции БТМБ. При этом используется несколько ядер для ускорения вычислений с

помощью технологии OpenMP. В некоторых работах [3] задачи сейсмики решаются с помощью технологий параллельного программирования на системах с распределенной памятью. Авторы [4] рассматривают распараллеливание метода конечных разностей применительно к задачам сейсмики. В работе [5] сравниваются различные технологии написания параллельных программ. Рассматривается решение гиперболических уравнений для моделирования физических процессов.

Иногда при решении такого рода задач возникает дисбаланс между потоками вычислений. Это происходит из-за граничных условий, вычисление которых может требовать большего (или меньшего) количества вычислений, чем вся остальная сетка. Для решения этой задачи может использоваться динамическая балансировка нагрузки. В работе [6] рассматривается такой подход при решении задач распространения сейсмических волн.

Помимо задач сейсмики, техники написания параллельных программ активно применяются для решения задач, где требуется решение других уравнений в частных производных. Например, еще одна работа [7] посвящена совместному использованию OpenMP и векторных инструкций AVX процессора. За счет них авторам удалось добиться ускорения в 57 раз на 10 ядрах по сравнению с одним ядром при решении уравнений мелкой воды и уравнений Эйлера. Авторами другой работы [8] был разработан гибридный метод решения задач гидродинамики с использованием различных технологий распараллеливания как для систем с общей памятью, так и для систем с распределенной памятью.

В NUMA-системах обращение к различным областям памяти может занимать разное время. Работа [9] посвящена этому вопросу. Авторы используют технику, которая заключается в размещении в локальной памяти процесса данных, к которым процесс выполняет наибольшее количество обращений. В результате было получено уменьшения времени выполнения, так как время обращения к локальной памяти быстрее.

В данной работе рассматриваются технологии OpenMP [10] и POSIX Threads. Они предназначены для написания параллельных программ в системах с общей памятью. Это значит, что обмен данными между процессами осуществляется за счет чтения и записи в память, общую для всех процессов. Сравнение этих технологий производилось на задаче решения уравнения линейной динамической теории упругости, которая описывается гиперболической системой уравнений. Эффективность распараллеливания измерялась в отношении времени выполнения последовательной программы к времени выполнения распараллеленной. Для повышения производительности были задействованы векторные инструкции процессора. Они позволяют выполнять несколько арифметических операций над данными за время, требующееся для выполнения только одной операции. Также потокам было запрещено перемещаться между ядрами процессора, чтобы повысить производительность в NUMA-системах.

2. Математическая модель

Для описания поведения среды использовалась модель идеального изотропного линейно-упругого материала. Приведённая ниже система дифференциальных уравнений в частных производных описывает состояние элементарного объёма упругого материала в приближении малых деформаций для двумерного случая:

дьх _ дахх + д(гху дуу _ доху + диуу ^ дЬ дх ду ^ дЬ дх ду ' д(хх + 2 + ^дЦу д(уу _ уд^ + 2 дщ, даху _ /ду^ + ОуЛ дЬ ^ дх ду ' дЬ дх ^ ду ' дЬ ^ \дх ду ) '

где р — плотность среды; — параметры Ламе; ух и уу — горизонтальная и вертикальная составляющие скорости частиц среды; 0ХХ'0УУ,оху — компоненты тензора напряжения. Данную систему можно переписать в матричной форме:

ди„ . диа „ диа , .

иг + А» тих + (1)

где u — вектор из 5 независимых переменных u = (axx,ayy,axy, vx, vy)T. Явный вид матриц Ард, Bpq представлен в [11]. Здесь и далее подразумевается суммирование по повторяющимся индексам. Собственные значения матриц Ард и Bpq таковы: Si = — ср, s2 = — cs, sз = 0, s4 = cs, s5 = ср, где ср и cs — скорости распространения продольных и поперечных волн в среде.

3. Численный метод

Применяя покоординатное расщепление, можно свести задачу построения разностной схемы для системы уравнений (1) к задаче построения разностной схемы для системы вида

дир . диа , .

ИТ + А»а? = °- (2)

Для гиперболической системы уравнений (2) матрицу A можно представить в виде A = RAR-1, где Л — диагональная матрица, элементы которой — собственные значения матрицы A, а R — матрица, состоящая из правых собственных векторов матрицы A. Введём новые переменные: w = R-1u (так называемые инварианты Римана). Тогда система уравнений (2) сведётся к системе из 5 независимых скалярных уравнений переноса.

Приведем схему третьего порядка точности для численного решения одномерного линейного уравнения переноса щ + аих = 0, а > 0,а = ат/h, т — шаг по времени, h — шаг по координате:

=ипт + а(До + Д2)/2 + а2(До — Д2)/2 + ( „ ) (Ai — 2Дс + Д2), (3)

,п+1 — „,п _L _L / О _L ГТ2(_ / О _L _И ,

6

До = ипт-1 —ипт, (4)

Д1 = ит-2 —ит-1, (5)

Д2 =иПт — и™+1- (6)

Схема (3) устойчива для чисел Куранта, не превышающих единицу. Используется сеточно-характеристический критерий монотонности, опирающийся на характеристическое свойство точного решения:

В местах выполнения данного критерия порядок схемы падает до второго.

После того как значения инвариантов Римана на следующем шаге по времени найдены, востанавливается решение: ига+1 = Rw"+1.

4. Постановка задачи

Тестовая модель изображена на рис. 1. Размеры области даны в километрах. На нижней и боковых границах устанавливалось неотражающее граничное условие, на верхней -свободная граница. В качестве источника возмущения использовалась вертикальная сила, приложенная к площадке с 925.7 м по 974.1 м на дневной поверхности, амплитуда которой описывалась импульсом Рикера частоты 40 Гц. Результаты расчета представлены на рис. 2.

5. Измерения

Все тесты производились на одной и той же тестовой задаче. При этом использовалась двумерная сетка размерами 4000 х 4000 узлов. Выполнялось 100 шагов по времени. Отдельно были измерены результаты для вычислений с одинарной и с двойной точностью. Каждый узел содержал по 5 переменных, поэтому в случае вычислений с одинарной точностью вся сетка занимала в памяти 305 Мбайт, а в случае двойной точности - 610 Мбайт.

О 0.5 1 1.5 2

Рис. 1. Геологическая модель антиклиналь- Рис. 2. Результат расчета - волновая карти-

ной ловушки [12] на в момент времени £ = 0.38 с

Т а б л и ц а 1

Процессоры, на которых производились измерения

Название (обозначение на графиках) Частота, ГГц Кеш, Кбайт Процессоры и ядра Компиляторы Производительность ядра/общая, Гфлопс

AMD Opteron 6272 (a64) 2.1 2048 8 процессоров по 8 ядер gcc-4.4.6 16.8/1075.2

Intel Xeon E5-2697 (i24) 2.7 30720 2 процессора по 12 ядер gcc-4.4.7, icc-15.0.0 21.6/518.4

AMD Opteron 8431 (a48) 2.4 512 8 процессоров по 6 ядер gcc-4.4.7 19.2/921.6

Использовались следующие флаги компиляции: -fno-tree-vectorize, чтобы запретить векторизацию там, где мы ее не хотим использовать, -fopenmp и -pthread для задействования той или иной техологии, -O2 для использования стандартных оптимизаций, -msse, -mavx, чтобы явно указать, какой набор инструкций мы хотим использовать при автовекторизации.

Измерения были произведены на процессорах, указанных в табл. 1. В дальнейшем обозначения используемых компиляторов на графиках - gcc и icc. Обозначения для используемого набора векторных инструкций процессора - sse или avx.

6. Процесс оптимизации

В первую очередь была написана последовательная версия алгоритма. Пересчет узлов производился в два шага: по оси X (обход сетки в начале по столбцам, потом по строкам) и по оси Y (обход сетки в начале по строкам, потом по столбцам). Было измерено количество арифметических операций в коде программы и из этого было получено теоретическое значение флопс - 262, требуемое для пересчета одного узла. На графиках эта версия обозначается simple.

Далее была предпринята попытка оптимизировать шаг по Y. На этом шаге был изменен порядок обращения к узлам так, что он совпал с обращением в шаге по X. Таким образом обращение к памяти стало последовательным, что привело к более эффективному использованию кеш-памяти. Кроме того, путем замены арифметических операций на эквивалентные было уменьшено требуемое количество флопс для пересчета узла до 190. Обозначение этой версии на графиках - cache-friendly.

400

200

□ 0 simple О □ cache-friendly

□ Ш

i24, icc i24, gcc a48, gcc a64, gcc а) одинарная точность

400

200

О □ simple D 0 cache-friendly

i24, icc i24, gcc a48, gcc a64, gcc б) двойная точность

Рис. 3. Время работы в секундах с оптимизацией кеш-памяти и без нее

Время работы первоначальной и оптимизированной версии алгоритма на различных процессорах представлено на рис. 3.

0

0

I I I I I I I I I I I | I I I I I | I I I I I | I I I I I | I I I I I | I I I I I | I I I I I | I I I I I |

icc, cache friendly, avx icc, cache friendly, sse icc, cache friendly gcc, cache friendly icc, simple gcc, simple

0 20 40 60 80 100 120 140 160 180 Время работы, с

Рис. 4. Время работы при использовании векторизации на Intel Xeon E5-2697, одинарная точность

icc, cache friendly, avx icc, cache friendly, sse icc, cache friendly gcc, cache friendly icc, simple gcc, simple

0 20 40 60 80 100 120 140 160 Время работы, с

Рис. 5. Время работы при использовании векторизации на Intel Xeon E5-2697, двойная точность

После этого были векторизованы циклы, в которых обрабатывались узлы сетки. Это было сделано с помощью директив (#pragma omp simd) технологии OpenMP, появившихся в стандарте 4.0. Явно было указано, какой набор инструкций должен использоваться: SSE или AVX. Измерения производились только на процессоре Intel Xeon E5-2697. Результаты представлены на рис. 4 и 5.

Исходя из частоты ядра каждого конкретного процессора, была определена его пиковая производительность. Таким образом, был вычислен процент от пиковой производительности для каждого теста (рис. 6, 7).

a48, gcc, cache friendly a48, gcc, simple a64, gcc, cache friendly a64, gcc, simple i24, gcc, cache friendly i24, icc, cache friendly i24, gcc, simple i24, icc, simple i24, icc, cache friendly, sse i24, icc, cache friendly, avx

5 10 15 20 Рис. 6. Процент от пиковой производительности последовательной версии, одинарная точность

a48, gcc, cache friendly a48, gcc, simple a64, gcc, cache friendly a64, gcc, simple i24, gcc, cache friendly i24, icc, cache friendly i24, gcc, simple i24, icc, simple i24, icc, cache friendly, sse i24, icc, cache friendly, avx

4 6 8 10 12 14 16 18 Рис. 7. Процент от пиковой производительности последовательной версии, двойная точность

■

7. Параллельная версия

Алгоритм распараллелен с использованием технологий OpenMP и POSIX Threads. При этом отдельно распараллелена каждая из версий, чтобы была возможность пронаблюдать последствия оптимизации последовательного кода на эффективность его распараллеливания.

—*-simple,gcc,aff simple,gcc cf,gcc,aff -^^cf,gcc

—♦— spl-gr,gcc --*--spl-gr,gcc,aff-a--pthread,gcc

Рис. 8. Ускорение и эффективность распараллеливания на AMD Opteron 8431 при вычислениях с одинарной точностью

Щ 20

к я

CD

а §

£ 10

—гаг / Ls / >f LIU T-г7*!-\— ^хЦтТПт V\V\ 1 1 Ml \ \ \ ж 1T hi J-tif J4-FH ^жит

1

л

н о

о

я и к

S 0.5

id

■В

О)

1.

К

10 20 30 40 Число потоков

10 20 30 40 Число потоков

simple,gcc,aff simple,gcc cf,gcc,aff -

spl-gr,gcc --*--spl-gr,gcc,aff-a--pthread,gcc

- cf,gcc

Рис. 9. Ускорение и эффективность распараллеливания на AMD Opteron 8431 при вычислениях с двойной точностью

40

§ 30

я

id

сР g 20

скУ

10 0

1

0.5

0 10 20 30 40 50 Число потоков

60

Число потоков

simple,gcc,aff simple,gcc -в— cf,gcc,aff -spl-gr,gcc --*--spl-gr,gcc,aff-o--pthread,gcc

cf,gcc

Рис. 10. Ускорение и эффективность распараллеливания на AMD Opteron 6272 при вычислениях с одинарной точностью

40

■d 30

и

я

id

а §

скУ

10

20

0

1

л

т с о

я

в и

кт 0.5

ID

■В

О)

0 10 20 30 40 50 Число потоков

60

0

%

*!

0 10 20 30 4

0 50 60

Число потоков

simple,gcc,aff simple,gcc cf,gcc,aff ^^cf,gcc

spl-gr,gcc --*--spl-gr,gcc,aff-i>-pthread,gcc

0

0

0

0

0

Рис. 11. Ускорение и эффективность распараллеливания на AMD Opteron 6272 при вычислениях с двойной точностью

15

10

1

♦

i r * 4

k' ■

ï \

rV \ S

w П| \

f Л m - i >

M tS & s

V*'

S *

s

✓

J г

✓

я

✓

и"

10 15 20 Число потоков

1

Л

т с о

я

в и

S 0.5

<d

■В

О)

25

I

Я \ \ *

V \ \ \ ■

ч \ \

> Л a- *

1

10 15 Число потоков

20

25

5

0

0

0

5

0

5

—•- simple,gcc -■- cf,gcc —®—simple,icc—*— cf,icc —♦— spl-gr,icc --•-- cf,icc,sse -t>-spl-gr,gcc-«-- cf,icc,avx -^-spl-gr,icc,avx--»--spl-gr,icc,sse

Рис. 12. Ускорение и эффективность распараллеливания на Intel Xeon E5-2697 при вычислениях с двойной точностью без использования привязки потоков к определенным ядрам

15

10

10

15

20

25

Число потоков

0.5

10

15

20

25

Число потоков

1

5

0

0

0

5

0

5

—•- simple,gcc -■- cf,gcc —®—simple,icc—*— cf,icc —♦— spl-gr,icc --•-- cf,icc,sse --n--spl-gr,gcc-^- cf,icc,avx -^-spl-gr,icc,avx--»--spl-gr,icc,sse

Рис. 13. Ускорение и эффективность распараллеливания на Intel Xeon E5-2697 при вычислениях с двойной точностью с использованием привязки потоков к определенным ядрам

Была сделана оптимизация уже параллельного кода. В ней каждый поток выделяет память, необходимую только для пересчета узлов, отданных в распоряжение именно этого потока. Теоретически это должно привести к росту производительности в NUMA-системах, так как память, выделенная конкретным потоком, согласно стандарту OpenMP, выделяется в локальной памяти потока, обращение к которой занимает наименьшее количество времени. Обозначение этой версии на графиках - spl-gr.

Отдельно были произведены измерения, в которых потокам было запрещено перемещаться между ядрами процессора. Это сделано для уменьшения временных затрат. Без этого оптимизация, связанная с разделением сетки на отдельные части для различных процессов, перестает иметь смысл. Обозначение на графиках - aff. Результаты тестов представлены на рис. 8-13.

7.2. POSIX Threads

Далее ОрепМР-версия кода с разделением сетки между потоками была изменена таким образом, чтобы использовать POSIX Threads в качестве технологии распараллеливания (рис. 14).

—•- pthread,gcc -■- pthread,icc,avx pthread,icc —к-- spl-gr,gcc —♦— spl-gr,icc,avx spl-gr,icc

Рис. 14. Сравнение ускорения и эффективности распараллеливания технологий OpenMP и POSIX Threads на Intel Xeon E5-2697

8. Заключение

Наибольшая производительность от пиковой для данного процессора была получена в последовательной версии 22% для вычислений с одинарной точностью и 17% для вычислений с двойной точностью с оптимизацией доступа к кеш-памяти и использованием векторных инструкций avx.

Максимальные полученные ускорения: на AMD Opteron 6272 - в 37 раз на 64 ядрах, на AMD Opteron 8431 - в 25 раз на 48 ядрах и на Intel Xeon E5-2697 - в 17 раз на 24 ядрах.

Таким образом, удалось продемонстрировать, что в системах с общей памятью основным ограничением для ускорения с ростом числа потоков является максимальная скорость взаимодействия с памятью. Разделение памяти между потоками позволяет не обращаться ядрам процессора к нелокальной памяти, а значит, доступ к памяти осуществляется быстрее. Привязка потоков к определенным ядрам уменьшает количество обращений в память, используемых для поддержания когерентности кешей ядер, и также ведет к увеличению производительности. Следует заметить, что реализация параллельной версии с использованием POSIX Threads оказывается сложнее, чем при использовании OpenMP, но при этом не дает значительного прироста производительности или эффективности распараллеливания.

Исследование выполнено при финансовой поддержке РФФИ в рамках научного проекта № 15-37-20673.

Литература

1. Caserta A., Ruggiero V., Lanucara P. Numerical modelling of dynamical interaction between seismic radiation and near-surface geological structures: a parallel approach // Computers & geosciences. 2002. V. 28, N 9. P. 1069-1077.

2. Araya-Polo M. [et al.] 3D seismic imaging through reverse-time migration on homogeneous and heterogeneous multi-core processors // Scientific Programming. 2009. V. 17, N 1-2. P. 185-198.

3. Петров И.Б., Хохлов Н.И. Моделирование задач 3D сейсмики на высокопроизводительных вычислительных системах // Математическое моделирование. 2014. Т. 6, № 1. С. 83-95.

4. Bohlen T. Parallel 3-D viscoelastic finite difference seismic modelling / /Computers & Geosciences. 2002. V. 28, N 8. P. 887-899.

5. Rostrup S., De Sterck H. Parallel hyperbolic PDE simulation on clusters: Cell versus GPU // Computer Physics Communications. 2010. V. 181, № 12. P. 2164-2179.

6. Keller Tesser R. [et al.] Improving the performance of seismic wave simulations with dynamic load balancing // Parallel, Distributed and Network-Based Processing (PDP), 2014 22nd Euromicro International Conference on. IEEE, 2014. P. 196-203.

7. Liu J. Y. [et al.] Hybrid OpenMP/AVX acceleration of a Split HLL Finite Volume Method for the Shallow Water and Euler Equations // Computers & Fluids. 2015. V. 110. P. 181188.

8. Schive H. Y., Zhang U. H., Chiueh T. Directionally unsplit hydrodynamic schemes with hybrid MPI/OpenMP/GPU parallelization in AMR // International Journal of High Performance Computing Applications. 2012. V. 26, N 4. P. 367-377.

9. Lof H., Holmgren S. Affinity-on-next-touch: increasing the performance of an industrial PDE solver on a cc-NUMA system // Proceedings of the 19th annual international conference on Supercomputing. ACM, 2005. P. 387-392.

10. Dagum L., Enon R. OpenMP: an industry standard API for shared-memory programming // Computational Science & Engineering, IEEE. 1998. V. 5, N 1. P. 46-55.

11. LeVeque R.J. Finite volume methods for hyperbolic problems. Cambridge university press, 2002. V. 31.

12. Jose M. [et al] Review Article: Seismic modeling // Geophysics. 2002. V. 67, N 4. P. 13041325.

References

1. Caserta A., Ruggiero V., Lanucara P. Numerical modelling of dynamical interaction between seismic radiation and near-surface geological structures: a parallel approach. Computers & geosciences. 2002. V. 28, N 9. P. 1069-1077.

2. Araya-Polo M. [et al.] 3D seismic imaging through reverse-time migration on homogeneous and heterogeneous multi-core processors. Scientific Programming. 2009. V. 17, N 1-2. P. 185-198.

3. Petrov I.B., Khokhlov N.I. Modeling 3D seismic problems using high-performance computing systems. Mathematical Models and Computer Simulations. 2014. V. 6, N 4. P. 342-350. (in Russian).

4. Bohlen T. Parallel 3-D viscoelastic finite difference seismic modelling. Computers & Geosciences. 2002. V. 28, N 8. P. 887-899.

5. Rostrup S., De Sterck H. Parallel hyperbolic PDE simulation on clusters: Cell versus GPU. Computer Physics Communications. 2010. V. 181, N 12. P. 2164-2179.

6. Keller Tesser R. [et al.] Improving the performance of seismic wave simulations with dynamic load balancing. Parallel, Distributed and Network-Based Processing (PDP), 2014 22nd Euromicro International Conference on. IEEE, 2014. P. 196-203.

7. Liu J. Y. [et al.] Hybrid OpenMP/AVX acceleration of a Split HLL Finite Volume Method for the Shallow Water and Euler Equations. Computers & Fluids. 2015. V. 110. P. 181-188.

8. Schive H.Y., Zhang U.H., Chiueh T. Directionally unsplit hydrodynamic schemes with hybrid MPI/OpenMP/GPU parallelization in AMR. International Journal of High Performance Computing Applications. 2012. V. 26, N 4. P. 367-377.

9. Lof H., Holmgren S. Affinity-on-next-touch: increasing the performance of an industrial PDE solver on a cc-NUMA system. Proceedings of the 19th annual international conference on Supercomputing. ACM, 2005. P. 387-392.

10. Dagum L, Enon R. OpenMP: an industry standard API for shared-memory programming. Computational Science & Engineering, IEEE. 1998. V. 5, N 1. P. 46-55.

11. LeVeque R.J. Finite volume methods for hyperbolic problems. Cambridge university press, 2002. V. 31.

12. Jose M. et al. Review Article: Seismic modeling. Geophysics. 2002. V. 67, N 4. P. 1304-1325.

Поступила в редакцию 10.05.2016