Математическое моделирование прочности строительных конструкций на гибридных вычислительных системах
Розглядаються алгоритми розв’язування на комп’ютерах гібридної архітектури (на багатоядерних процесорах і графічних прискорювачах) систем лінійних алгебраїчних рівнянь зі стрічковими симетричними додатно-визначеними матрицями, які виникають при чисельному моделюванні з використанням методу скінченни...
Gespeichert in:
| Datum: | 2017 |
|---|---|
| Hauptverfasser: | , , , |
| Format: | Artikel |
| Sprache: | Russian |
| Veröffentlicht: |
Інститут кібернетики ім. В.М. Глушкова НАН України
2017
|
| Schriftenreihe: | Проблемы управления и информатики |
| Schlagworte: | |
| Online Zugang: | https://nasplib.isofts.kiev.ua/handle/123456789/208571 |
| Tags: |
Tag hinzufügen
Keine Tags, Fügen Sie den ersten Tag hinzu!
|
| Назва журналу: | Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| Zitieren: | Математическое моделирование прочности строительных конструкций на гибридных вычислительных / А.Ю. Баранов, А.В. Попов, Я.Е. Слободян, А.Н. Химич // Проблемы управления и информатики. — 2017. — № 4. — С. 68-81. — Бібліогр.: 15 назв. — рос. |
Institution
Digital Library of Periodicals of National Academy of Sciences of Ukraine| id |
irk-123456789-208571 |
|---|---|
| record_format |
dspace |
| spelling |
irk-123456789-2085712025-11-02T01:15:59Z Математическое моделирование прочности строительных конструкций на гибридных вычислительных системах Математичне моделювання міцності будівельних конструкцій на гібридних обчислювальних системах Mathematical modeling of strength of building structures for hybrid computing systems Баранов, А.Ю. Попов, А.В. Слободян, Я.Е. Химич, А.Н. Конфликтно-управляемые процессы и методы принятия решений Розглядаються алгоритми розв’язування на комп’ютерах гібридної архітектури (на багатоядерних процесорах і графічних прискорювачах) систем лінійних алгебраїчних рівнянь зі стрічковими симетричними додатно-визначеними матрицями, які виникають при чисельному моделюванні з використанням методу скінченних елементів, просторових будівельних конструкцій. Наведено результати апробації запропонованих гібридних алгоритмів на розв’язуванні задач статичного аналізу міцності декількох будівельних об’єктів і окремих конструкцій. Algorithms for solving on computers of hybrid architecture (with multi-core processors and graphics accelerators) the systems of linear algebraic equations with band symmetrical positive definite matrices, arising in the numerical modeling of 3D building constructions using the finite element method is considered. Some results of the proposed hybrid algorithms testing for solving the problems of static strength analysis of construction sites and individual constructions are presented. 2017 Article Математическое моделирование прочности строительных конструкций на гибридных вычислительных / А.Ю. Баранов, А.В. Попов, Я.Е. Слободян, А.Н. Химич // Проблемы управления и информатики. — 2017. — № 4. — С. 68-81. — Бібліогр.: 15 назв. — рос. 0572-2691 https://nasplib.isofts.kiev.ua/handle/123456789/208571 519.6:539.3 10.1615/JAutomatInfScien.v49.i7.20 ru Проблемы управления и информатики application/pdf Інститут кібернетики ім. В.М. Глушкова НАН України |
| institution |
Digital Library of Periodicals of National Academy of Sciences of Ukraine |
| collection |
DSpace DC |
| language |
Russian |
| topic |
Конфликтно-управляемые процессы и методы принятия решений Конфликтно-управляемые процессы и методы принятия решений |
| spellingShingle |
Конфликтно-управляемые процессы и методы принятия решений Конфликтно-управляемые процессы и методы принятия решений Баранов, А.Ю. Попов, А.В. Слободян, Я.Е. Химич, А.Н. Математическое моделирование прочности строительных конструкций на гибридных вычислительных системах Проблемы управления и информатики |
| description |
Розглядаються алгоритми розв’язування на комп’ютерах гібридної архітектури (на багатоядерних процесорах і графічних прискорювачах) систем лінійних алгебраїчних рівнянь зі стрічковими симетричними додатно-визначеними матрицями, які виникають при чисельному моделюванні з використанням методу скінченних елементів, просторових будівельних конструкцій. Наведено результати апробації запропонованих гібридних алгоритмів на розв’язуванні задач статичного аналізу міцності декількох будівельних об’єктів і окремих конструкцій. |
| format |
Article |
| author |
Баранов, А.Ю. Попов, А.В. Слободян, Я.Е. Химич, А.Н. |
| author_facet |
Баранов, А.Ю. Попов, А.В. Слободян, Я.Е. Химич, А.Н. |
| author_sort |
Баранов, А.Ю. |
| title |
Математическое моделирование прочности строительных конструкций на гибридных вычислительных системах |
| title_short |
Математическое моделирование прочности строительных конструкций на гибридных вычислительных системах |
| title_full |
Математическое моделирование прочности строительных конструкций на гибридных вычислительных системах |
| title_fullStr |
Математическое моделирование прочности строительных конструкций на гибридных вычислительных системах |
| title_full_unstemmed |
Математическое моделирование прочности строительных конструкций на гибридных вычислительных системах |
| title_sort |
математическое моделирование прочности строительных конструкций на гибридных вычислительных системах |
| publisher |
Інститут кібернетики ім. В.М. Глушкова НАН України |
| publishDate |
2017 |
| topic_facet |
Конфликтно-управляемые процессы и методы принятия решений |
| url |
https://nasplib.isofts.kiev.ua/handle/123456789/208571 |
| citation_txt |
Математическое моделирование прочности строительных конструкций на гибридных вычислительных / А.Ю. Баранов, А.В. Попов, Я.Е. Слободян, А.Н. Химич // Проблемы управления и информатики. — 2017. — № 4. — С. 68-81. — Бібліогр.: 15 назв. — рос. |
| series |
Проблемы управления и информатики |
| work_keys_str_mv |
AT baranovaû matematičeskoemodelirovaniepročnostistroitelʹnyhkonstrukcijnagibridnyhvyčislitelʹnyhsistemah AT popovav matematičeskoemodelirovaniepročnostistroitelʹnyhkonstrukcijnagibridnyhvyčislitelʹnyhsistemah AT slobodânâe matematičeskoemodelirovaniepročnostistroitelʹnyhkonstrukcijnagibridnyhvyčislitelʹnyhsistemah AT himičan matematičeskoemodelirovaniepročnostistroitelʹnyhkonstrukcijnagibridnyhvyčislitelʹnyhsistemah AT baranovaû matematičnemodelûvannâmícnostíbudívelʹnihkonstrukcíjnagíbridnihobčislûvalʹnihsistemah AT popovav matematičnemodelûvannâmícnostíbudívelʹnihkonstrukcíjnagíbridnihobčislûvalʹnihsistemah AT slobodânâe matematičnemodelûvannâmícnostíbudívelʹnihkonstrukcíjnagíbridnihobčislûvalʹnihsistemah AT himičan matematičnemodelûvannâmícnostíbudívelʹnihkonstrukcíjnagíbridnihobčislûvalʹnihsistemah AT baranovaû mathematicalmodelingofstrengthofbuildingstructuresforhybridcomputingsystems AT popovav mathematicalmodelingofstrengthofbuildingstructuresforhybridcomputingsystems AT slobodânâe mathematicalmodelingofstrengthofbuildingstructuresforhybridcomputingsystems AT himičan mathematicalmodelingofstrengthofbuildingstructuresforhybridcomputingsystems |
| first_indexed |
2025-11-02T02:07:21Z |
| last_indexed |
2025-11-03T02:05:19Z |
| _version_ |
1847733014298099712 |
| fulltext |
© А.Ю. БАРАНОВ. А.В. ПОПОВ, Я.Е. СЛОБОДЯН, А.Н. ХИМИЧ, 2017
68 ISSN 0572-2691
УДК 519.6:539.3
А.Ю. Баранов, А.В. Попов, Я.Е. Слободян, А.Н. Химич
МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ
ПРОЧНОСТИ СТРОИТЕЛЬНЫХ КОНСТРУКЦИЙ
НА ГИБРИДНЫХ ВЫЧИСЛИТЕЛЬНЫХ СИСТЕМАХ
Введение
Средства численного моделирования различных физических процессов, ис-
пользующие метод конечных элементов (МКЭ), сегодня применяются при проек-
тировании и разработке в чрезвычайно широком диапазоне — начиная от проек-
тирования корпуса мобильного телефона и заканчивая разработкой сейсмо-
стойких объектов или ядерных реакторов. При этом большую часть времени
моделирования занимает нахождение решения систем линейных алгебраических
уравнений (СЛАУ), возникающих при дискретизации исходных систем диффе-
ренциальных уравнений. Следует заметить, что обычно моделирование предпола-
гает многократное решение СЛАУ, поскольку интерес представляет динамика ис-
следуемого процесса во времени, причем количество итераций может составлять
сотни и тысячи. В случае нелинейных задач — это также многократное построе-
ние и решение СЛАУ в пределах одного шага по времени.
Таким образом, следует отметить чрезвычайное влияние скорости работы про-
граммы решения СЛАУ на общую скорость работы конечно-элементного решателя.
Еще одним важным фактором, определяющим применимость решателя СЛАУ
для использования в современных задачах моделирования методом конечных
элементов, является используемая схема размещения СЛАУ в оперативной памя-
ти. Задачи проектирования сложных конструкций часто приводят к построению
СЛАУ с миллионами (или десятками миллионов) неизвестных. Например, для
системы порядка n = 1217610, построенной для численного анализа напряженно-
деформируемого состояния многоэтажной конструкции, полуширина ленты мат-
рицы после применения обратного алгоритма перенумерации Катхилла–Макки
составляла m = 3628. Для ее хранения необходим объем оперативной памяти в
33 Гбайт (для 4417489080 элементов, с учетом симметричности матрицы, хранит-
ся только ее половина).
При решении СЛАУ большого порядка (105–108) с полушириной ленты мат-
рицы O(√n) (т.е. 102–104) количество арифметических операций с плавающей за-
пятой (при использовании метода Холецкого) составляет O(109–1016). Общее ко-
личество хранимых элементов матрицы системы — 108–1012. Такие требования к
вычислениям намного опережают возможности традиционных параллельных ком-
пьютеров, даже несмотря на многоядерность процессоров.
Решение проблемы ускорения вычислений на многоядерных компьютерах
может достигаться использованием сопроцессоров-ускорителей для выполнения
больших объемов однородных арифметических операций. Раньше в качестве ускори-
телей использовались преимущественно графические процессоры (Graphical Proces-
sing Unit — GPU). В последнее время более высокая производительность чаще
достигается при использовании процессоров Intel Xeon Phi в качестве таких сопро-
цессоров. Ускорение вычислений в ближайшей перспективе, на наш взгляд, возмож-
но при использовании гибридных систем, которые совмещают многоядерные про-
цессоры (Central Processing Unit — CPU) с сопроцессорами-ускорителями
(см. рейтинг Top500 [1]), в частности с GPU.
Международный научно-технический журнал
«Проблемы управления и информатики», 2017, № 4 69
В статье представлены результаты численного моделирования пространст-
венных строительных конструкций с использованием в качестве процессорной
компоненты программного обеспечения на основе нового гибридного алгоритма
решения СЛАУ с ленточной симметричной положительно-определенной матри-
цей, совмещающего вычисления на многоядерных процессорах и графических ус-
корителях. Приведены многочисленные результаты с использованием для пре-
процессорной и постпроцессорной обработки гибридного варианта программного
комплекса Лира–Сапр [2].
1. Моделирование напряженно-деформированного
состояния строительных конструкций
Математические постановки задач. Математически задачи расчета прочно-
сти конструкций с использованием принципа возможных перемещений могут
быть поставлены в виде следующих вариационных задач [3].
Необходимо найти вектор-функцию ,0Uu∈ удовлетворяющую для любой век-
тор-функции 0Uv∈ (любого возможного перемещения) интегральному тождеству:
• для статической задачи
);,(),( vflvua = (1)
• для динамической задачи
),,(),(),(),( vflvucvubvua =′+′′+ (2)
;)(,)( )1(
0
)0(
0 utuutu =′= (3)
• для задачи о собственных колебаниях
),,(),( vubvua λ= (4)
где U0 — бесконечномерное функциональное пространство возможных переме-
щений; симметричные билинейные формы ),(),,(),,( vucvubvua ′′′ пропорцио-
нальны соответственно потенциальной, кинетической энергиям деформации и
энергии торможения, а линейная форма ),( vfl пропорциональна работе прило-
женных (внешних) усилий при нагружении; u′ — первая производная вектор-функ-
ции u по времени; u ′′ — вторая.
Дискретизация задач методом конечных элементов. Теоретической осно-
вой Лира–Сапр является метод конечных элементов [4], реализованный в форме
перемещений. Выбор именно этой формы объясняется простотой ее алгоритмиза-
ции и физической интерпретации, наличием единых методов построения матриц
жесткости и векторов нагрузок для различных типов конечных элементов, воз-
можностью учета произвольных граничных условий и сложной геометрии рас-
считываемой конструкции.
Для получения с помощью МКЭ дискретной задачи область разбивается
на конечные элементы, назначаются узлы и их степени свободы (перемещения
и углы поворота узлов). Степеням свободы соответствуют базисные (коорди-
натные, аппроксимирующие) вектор-функции ,iϕ отличные от нуля только на
элементах, которые содержат соответствующий данной степени свободы узел.
Кроме того, для степеней свободы и базисных функций справедливы следую-
щие соотношения:
,)( ijijL δ=ϕ (5)
где ijδ — символ Кронекера, а результатом операции )( ijL ϕ является значение
компоненты вектор-функции iϕ для степени свободы .jL
70 ISSN 0572-2691
Приближенные решения соответствующих задач ищутся в конечномерном
подпространстве hU0 пространства .0U Вектор-функции из подпространства hU0
являются кусочно-полиномиальными и могут быть представлены в виде удовле-
творяющей главным (кинетическим) условиям линейной комбинации базисных
вектор-функций
∑ χϕχ
N
=i
iih x=u
1
),()( (6)
где ),,2,1( Nii K=ϕ — упомянутый выше кусочно-полиномиальный базис .0
hU
Тогда разрешающие дискретные задачи имеют вид:
• системы линейных алгебраических уравнений
bAx = (7)
для статической задачи (1);
• задачи с начальными условиями
)1(
0
)0(
0 )(,)(),()()()( xtxxtxtbtAxtxCtxB =′==+′+′′ (8)
для динамической задачи (2), (3);
• алгебраической проблемы собственных значений (АПСЗ)
BxAx hλ= (9)
для задачи о собственных колебаниях (4).
Разрешающие задачи (7)–(9) имеют такие особенности:
1) высокий порядок матриц разрешающих дискретных задач — от 100 000 до
десятков миллионов;
2) элементы матриц и векторов разрешающих задач вычислены с погрешно-
стями, которые вызваны погрешностями исходных данных, дискретизации и вы-
числения значений этих элементов в компьютере;
3) матрицы дискретных задач — жесткостей A, масс B, демпфирования C —
являются
• симметричными;
• положительно-определенными или положительно-полуопределенными;
• разреженной структуры (ленточной, профильной, блочной и т.п.).
Матрица называется разреженной, если число ее элементов, отличных от 0,
составляет kn и ,nk << где n — порядок матрицы.
Важным преимуществом МКЭ является то, что элементы матриц и векторов
правых частей задач (7)–(9) получаются суммированием соответствующих эле-
ментов матриц жесткости (масс, демпфирования) и векторов нагрузок, построен-
ных для отдельных конечных элементов.
Такое свойство матриц и векторов дискретных задач МКЭ позволяет эффек-
тивно распараллелить процесс формирования этих матриц и векторов. Например,
область может быть разбита на приблизительно равные подобласти (по числу ис-
пользуемых процессоров), и для конечных элементов, содержащихся в этих по-
добластях, независимо друг от друга формируются соответствующие матрицы и
векторы с использованием обменов между процессорами, формируются глобаль-
ные матрицы и векторы соответствующей дискретной задачи (7)–(9). При этом
глобальные матрицы и векторы могут быть или сформированы в файле (файлах) на
внешнем носителе, или распределены между процессорами в соответствии с требо-
ваниями метода, который будет использоваться для решения дискретной задачи.
Международный научно-технический журнал
«Проблемы управления и информатики», 2017, № 4 71
В целях минимизации количества арифметических операций при решении
дискретных задач, а также необходимой оперативной и внешней памяти, учиты-
вая разреженную структуру матриц задач (7)–(9), в большинстве случаев прово-
дится переупорядочение неизвестных задач. В зависимости от критерия миними-
зации — ширины ленты матрицы, ее профиля, общего количества арифметиче-
ских операций при треугольном разложении матрицы и т.д. — используются
разные методы переупорядочения, например обратный алгоритм Катхилла–Мак-
ки, фактор деревьев, параллельных сечений, минимальной степени. Конкретные
рекомендации выбора метода переупорядочения дать сложно, потому что эффек-
тивность того или другого алгоритма существенно зависит от исходной структу-
ры конкретной матрицы.
Решение дискретных задач. Для решения дискретных задач (7)–(9) исполь-
зуются параллельные алгоритмы [5] традиционных методов решения соответст-
вующих задач.
Решение систем линейных алгебраических уравнений. Для решения
СЛАУ (7) — исходной или переупорядоченной — целесообразно использовать
прямой метод Холецкого. При этом полученное треугольное разложение матрицы
системы может быть сохранено в файле и использовано в последующих расчетах
прочности конструкции для других нагрузок (т.е. правых частей), что существен-
но уменьшит общее время расчета прочности конструкции.
Решение задачи с начальными условиями. В полудискретном МКЭ при-
ближенное решение ищется в виде (6), где коэффициенты ix являются функция-
ми времени t. В результате получаем систему обыкновенных дифференциальных
уравнений второго порядка с начальными условиями (8), где )1()0( ,),( xxtx — век-
торы с элементами ),(txi ),( )0()0( uLx ii = ).( )1()1( uLx ii =
Для решения задачи с начальными условиями (8) используется метод разло-
жения по формам собственных колебаний, т.е. по собственным векторам соответст-
вующей задачи (9). Приближенное решение представляется в виде линейной ком-
бинации нескольких собственных векторов, соответствующих наименьшим собст-
венным значениям этой задачи. Если kk z,λ ( ,T
jkkj Bzz δ= nkj ,,2,1, K= ) —
решение алгебраической проблемы собственных значений (9), то, считая в (8)
,)()(
1
∑
n
=k
kk zty=tx получим (при некоторых предположениях относительно мат-
рицы демпфирования С) систему, которая распадается на независимые относи-
тельно )(tyi уравнения:
,y=tyy=ty
,t>t,tP=ty+ty+ty
kkkk
kk
2
kkkkk
(1)
0
(0)
0
0
)(,)(
)()()(2)(
′
ω′ωξ′′
(10)
где ,10 <ξ< k ,2
kk λ=ω ),()( T tbztP kk = ,)( T)0()0(
kk Bzxy = ,)( T)1()1(
kk Bzxy =
.,,2,1 nk K= Анализ решений системы (10) свидетельствует о том, что сущест-
венный вклад в решение )(tx задачи (8) дают лишь 10–20 слагаемых, соответст-
вующих минимальным собственным значениям.
Решение алгебраической проблемы собственных значений. Для решения
частичной обобщенной АПСЗ (9) используется метод итераций на подпростран-
стве [6], позволяющий вычислить заданное количество минимальных собствен-
ных значений задачи и соответствующих им собственных векторов. Этот метод
72 ISSN 0572-2691
заключается в построении последовательности подпространств ),,2,1( K=tEt
сходящейся к подпространству ,∞E которое содержит искомые собственные век-
торы. Предполагается, что матрицы задачи (9) симметричные и одна из них (обычно
матрица A) положительно-определенная.
На t-й итерации необходимо выполнить следующие операции: (i) вычисление
решения СЛАУ ;1−= tt YAX (ii) вычисление проекции 1
T
−= ttt YXA матрицы А на
подпространство ;tE (iii) вычисление прямоугольной матрицы ;tt BXW = (iv) вы-
числение проекции ttt WXB T= матрицы B на подпространство ;tE (v) решение
полной АПСЗ для проекций ;ttttt ZBZA Λ= (vi) вычисление очередного прибли-
жения ;ttt ZWY = (vii) проверка условий окончания итерационного процесса. Если
эти условия выполняются после t итераций, то приближенным решением задачи (9)
считаем )1(* +λ=λ t
ii ,,,2,1( ri K= в предположении, что собственные значения
упорядочены по возрастанию), .11
*
++= tt ZXX Матричная система 1−= tt YAX
решается методом Холецкого. Причем разложение матрицы выполняется один раз
до начала итерационного процесса, а на каждой итерации решаются две системы с
нижней и верхней треугольными матрицами.
Для сходимости итерационного процесса важно, чтобы матрицы проекций
tA и tB были положительно определены. Однако при математическом модели-
ровании прочности конструкций матрица массы может быть положительно полу-
определена. Поэтому необходимо так выбирать начальное подпространство ,0E
чтобы матрица 1B была положительно определена.
Решение нелинейных задач выполняется итерационно — на каждой итерации
формируется новая СЛАУ с разреженной симметричной матрицей жесткости. Эта
система решается методом Холецкого. После этого проверяются условия оконча-
ния итерационного процесса, и, если заданная точность не достигнута, на основе
полученного приближения к решению модифицируются данные для формирова-
ния новой системы.
Таким образом, решение нелинейных задач, а также задач (8) и (9) базируется
на решении СЛАУ вида (7) с разреженной симметричной матрицей, основой ре-
шения которой служит разложение этой матрицы в произведение нижней и верх-
ней треугольных матриц. Следовательно, основные параметры (схемы распреде-
ления между процессорными устройствами, хранения и обработки данных) и эф-
фективность параллельных алгоритмов решения задач (7)–(9) в первую очередь
определяются качеством параллельных алгоритмов разложения методом Холец-
кого разреженных симметричных матриц (жесткости) различных структур.
В зависимости от структуры и заполненности матрицы, а также архитектуры
компьютера используются
— блочные параллельные алгоритмы для исследования и решения СЛАУ с
узкими ленточными матрицами и блочно-диагональными матрицами с окаймле-
нием [5];
— одномерные блочные циклические параллельные алгоритмы для иссле-
дования и решения СЛАУ с ленточными, профильными или «небоскребными»
(у которых не принимают участие в вычислениях внутренние элементы про-
филя матрицы, остающиеся нулевыми в процессе разложения матрицы) мат-
рицами [5, 7–9].
Международный научно-технический журнал
«Проблемы управления и информатики», 2017, № 4 73
Далее подробно рассмотрим блочные циклические алгоритмы решения на ком-
пьютерах гибридной архитектуры СЛАУ (7) с ленточной симметричной матрицей.
2. Гибридный алгоритм решения СЛАУ
с симметричной положительно-определенной ленточной матрицей
Рассмотрим решение СЛАУ (в общем случае матричной) вида (7)
BAX = (11)
порядка n с ленточной (ширина ленты )12 +m симметричной положительно-
определенной матрицей и q правыми частями.
Как известно, наиболее эффективным прямым методом решения задачи (11)
является метод Холецкого. В этом случае задача (11) разбивается на три под-
задачи:
• факторизация (LLT- или LDLT -разложение) исходной матрицы системы
TLLA = или ;TLDLA = (12)
• решение СЛАУ с нижней треугольной матрицей
;BLY = (13)
• решение СЛАУ с верхней треугольной матрицей
YXL =T или .T YXDL = (14)
Схемы распределения и хранения данных. Проблему эффективного распа-
раллеливания алгоритмов, в основе которых лежит факторизация исходной мат-
рицы, удалось решить с использованием циклических схем распределения и хра-
нения элементов матриц и векторов. Эти схемы были предложены независимо
друг от друга и практически одновременно несколькими группами исследовате-
лей [10, 11]. При использовании строчно-циклической схемы элементы матрицы
распределяются циклически между p процессорами (потоками, процессами), уча-
ствующими в вычислениях: если в (i–1)-м процессоре располагаются элементы
r-й строки матрицы, то в следующем (i–м) — элементы (r+1)-й строки.
Однако более эффективными являются блочно-циклические схемы распреде-
ления и хранения элементов матриц [5, 7, 8]. Например, при использовании одно-
мерной блочно-циклической схемы каждому CPU распределяются не отдельные
строки, а полосы из s подряд идущих строк, а операции с отдельными элементами
матрицы заменяются операциями с квадратными блоками порядка s или с прямо-
угольными блоками из s строк. С одной стороны, такой подход позволяет умень-
шить коммуникационные потери, а с другой, — обеспечивает значительно боль-
ший объем однородных вычислений, что приобретает особое значение при ис-
пользовании GPU.
Пусть гибридный алгоритм решения задачи (11) реализуется на p графиче-
ских процессорах и таком же количестве CPU, виртуальные связи между которы-
ми обеспечивают возможность передачи данных из одного вычислительного уст-
ройства в любое другое. Вычислительные устройства имеют логические номера
.1,,1,0 −pK При этом GPU с логическим номером i является сопроцессором-
ускорителем CPU с таким же логическим номером.
Для простоты изложения предположим, что порядок ленточной симметрич-
ной матрицы системы и полуширина ее ленты кратны порядку матричных бло-
ков s, т.е. Nsn = и .Msm = Параметр s должен выбираться таким образом, чтобы
блок полностью помещался в кеш-памяти как CPU, так и GPU. Таким образом,
74 ISSN 0572-2691
ленточная симметричная матрица разбивается на 2N квадратных блоков размера
,ss× причем каждая строка содержит не более 12 +M ненулевых блоков. По-
скольку матрица системы симметричная, для решения достаточно хранить только
ее диагональные элементы и элементы ее нижнего или верхнего треугольника
(в зависимости от варианта алгоритма), находящиеся в ненулевых блоках. Анало-
гичное представление используется и для нижней треугольной ленточной матри-
цы разложения L (рис. 1) или верхней треугольной матрицы LT.
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎟
⎠
⎞
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎜
⎝
⎛
−
−−
NN,NN,
NN
+M
+M+M
LL
L
L
LL
LL
LL
L
1
11,
2,2
1,21,1
3,23,1
2,22,1
1,1
00
0
0000
0
000
0000
L
LMM
O
MMMM
O
L
L
Рис. 1
Элементы исходной матрицы распределяются между CPU согласно следую-
щей схеме: блок JIA , хранится в CPU с логическим номером pI mod)1( − . Таким
же образом распределяются блоки матрицы L. Прямоугольные матрицы B и X де-
лятся на прямоугольные блоки размера .qs× Эти блоки распределяются между
CPU согласно той же блочно-циклической схеме, что и блоки матриц A и L.
Гибридный алгоритм LLT-разложения ленточной симметричной матрицы.
Блочный алгоритм факторизации состоит из N шагов. Будем называть CPU и GPU
ведущими на данном шаге, если блок IIL , или IIA , хранится в их памяти. Для
обозначения блоков матриц A и L, находящихся в памяти GPU, используется
верхний индекс d.
На I-м шаге вычисляются ненулевые блоки I-го столбца матрицы L и моди-
фицируется квадратная подматрица порядка Ms или меньше, в верхнем левом
углу которой расположен блок .1,1 ++ IIA Таким образом, на I-м шаге алгоритма
факторизации обрабатываются только те блоки, которые находятся в I-м столбце
блоков и в M столбцах блоков, следующих за ним [6]. Из разбиения матрицы
(рис. 1) видно, что в каждом таком столбце имеется не больше 1+M ненулевых
блоков. Поэтому с учетом симметрии на каждом шаге алгоритма достаточно хра-
нить в памяти GPU не более 2/)2)(1( ++ MM блоков. Предполагается, что блоки
подматрицы матрицы A порядка ,)1( sM + расположенной в верхнем левом углу
матрицы, скопированы в память соответствующих GPU согласно правилу: блок
jIA , располагается в памяти GPU с номером .mod)12( pJI −−
Таким образом, алгоритм LLT-разложения ленточной симметричной матрицы
состоит из следующих действий (решений подзадач) для каждого .,,2,1 NI K=
(i) LLT-разложение блока IIA , на ведущем CPU с использованием последо-
вательного алгоритма
.T
,,, IIIIII LLA = (15)
Международный научно-технический журнал
«Проблемы управления и информатики», 2017, № 4 75
(ii) Рассылка блока IIL , из ведущего CPU всем GPU, участвующим в вы-
числениях.
(iii) Вычисление соответствующим GPU ненулевых блоков I-го столбца
,,
d
IJL где },,{min,,1 NMIIJ ++= K путем решения соответствующих матрич-
ных СЛАУ
.T
,
T
,, IJIJII ALL = (16)
(iv) Асинхронная 2s-ранговая модификация на GPU при условии
0mod)( =− pJK
,T1)()(
IK,IJ,
I
KJ,
I
KJ, LLA=A −− (17)
где },,min{,,1 NMIIJ ++= K }.,min{,, NMIJK += K
(v) Рассылка каждым CPU вычисленных при выполнении операции (iii) бло-
ков d
IJL , всем остальным вычислительным устройствам.
(vi) 2s-ранговая модификация на соответствующем GPU (при условии
)0mod)( >− pJK по формуле (17) блоков ,d
KJ,A где },,{min,,1 NMIIJ ++= K
}.,{min,, NMIJK += K
(vii) Асинхронное копирование всех блоков матрицы, которые находятся в
столбце I +M +2, в память соответствующих GPU, если I +M +2 < p.
(viii) Синхронизация вычисления блока d
+I+IA 11, на CPU с логическим номе-
ром .mod pI
(ix) Копирование блока d
+I+IA 11, в память соответствующего CPU.
Заметим, что условие 0mod)( =− pJK при выполнении операции (iv) алго-
ритма означает, что все блоки, необходимые для вычисления по формуле (16), на-
ходятся в одном и том же GPU. А при выполнении операции (vi) алгоритма мо-
дифицируются только те блоки, у которых разность индексов не делится на p.
Это означает, что необходимые блоки были получены при рассылке (v) от дру-
гих устройств.
Следует также отметить, что аналогичный алгоритм может быть исполь-
зован для LDLT-разложения ленточной симметричной матрицы. В этом случае
вычисления проводятся по формулам ,T
,,, IIIIII LGA = T
,
T
,, IJIJII AGL = и =A I
KJ,
)(
,T)1(
IK,IJ,
I
KJ, LGA −= − где ,,, IIJIJ DLG = позволяющим не увеличивать количество
арифметических операций, но требующим дополнительный объем памяти для
временного хранения блоков IJG , .
Решение СЛАУ с нижней треугольной ленточной матрицей. Вычисли-
тельная схема алгоритма решения задачи (13) BLY = во многом аналогична вы-
числительной схеме алгоритма факторизации. Этот алгоритм также состоит из
N шагов. На каждом шаге NI ,,2,1 K= выполняются следующие действия (ре-
шаются подзадачи), в результате которых вычисляется один блок решения.
(i) Решение СЛАУ )1( −I
III,I B=YL с нижней треугольной матрицей IIL , на
ведущем GPU.
(ii) Рассылка из ведущего GPU вычисленного блока решения IY всем осталь-
ным GPU.
76 ISSN 0572-2691
(iii) Модификация (s-ранговая) блоков правой части JB K,1( += IJ
}),{min, NMI +K по формуле IIJ,
I
J
I
J YLB=B −− )1()( на GPU в соответствии с их
распределением.
Решение СЛАУ с верхней треугольной ленточной матрицей. Вычисли-
тельная схема алгоритма решения задачи (14) YXL =T существенно отличается
от предыдущей, поскольку распределенными между процессорными устройства-
ми являются столбцы блоков матрицы системы. Этот алгоритм также состоит из
N шагов. На I-шаге ( NI ,,2,1 K= ) выполняются следующие действия (решаются
подзадачи), в результате которых вычисляется один блок решения системы (здесь
).1+−= INK
(i) Вычисление при условиях 1>I и 0mod)( =− pJK каждым GPU матрич-
ных произведений JKJ, XLT ( },{min,,1 NMKKJ ++= K ) в соответствии с рас-
пределением блоков .JX
(ii) Вычисление (при pI > ) каждым GPU сумм вычисленных ранее матрич-
ных произведений.
(iii) Мультисборка ведущим CPU (при )pI > сумм матричных произведений
и вычисление ∑
=
++
− −=
KM
J
JKKJKK
I
K XLYY
1
T
,
)1( , .},{min KMMM K −=
(iv) Решение СЛАУ 1)(T −I
KKKK, Y=XL с верхней треугольной матрицей T
KK,L
на ведущем GPU ( ).( )0(
II YY ≡
Эффективность алгоритма факторизации. Количество арифметических
операций, выполняемых при факторизации (12) матрицы системы (11), в qm раз
больше количества арифметических операций, выполняемых при решении тре-
угольной системы (13) или (14), поэтому эффективность рассматриваемого гиб-
ридного алгоритма решения СЛАУ с ленточной симметричной матрицей опреде-
ляется эффективностью алгоритма разложения матрицы системы.
Обозначим время выполнения одной арифметической операции на CPU как
,Ct а на GPU — .Gt Количество арифметических операций, которое GPU может
выполнять одновременно, обозначим on . Поскольку GPU реализует архитектуру
SIMD, среднее время, затраченное одним GPU на выполнение Q однородных
арифметических операций, можно оценить величиной .oG ntQ
Количество арифметических операций, выполняемых при решении k-й под-
задачи алгоритма, обозначим kO . Тогда справедливы следующие оценки:
,33
1 sO ≈ ,2
3 msO ≈ ,)(4 pmssp+mO ≈ .)1(2
7 ppsmO −≈ Остальные подзада-
чи связаны с обменами данными между процессорными устройствами. Время пе-
ресылки одного блока между двумя GPU обозначим .Bt Если для рассылки дан-
ных от одного GPU всем остальным используется алгоритм «дерева» [12], то об-
щее время мультирассылки можно оценить величиной .log2 ptB Обозначим
суммарное время копирования (ii) и (ix) блоков между CPU и GPU как .It Время
на выполнение копирования (vii) можно не учитывать, так как это копирование
проводится на фоне вычислений.
Таким образом время выполнения алгоритма на 1 CPU и 1 GPU оценивается
величиной
),2( 11 oGIC nts)+ms(m+t+tON=T (18)
Международный научно-технический журнал
«Проблемы управления и информатики», 2017, № 4 77
а время выполнения алгоритма на p CPU и p GPU величиной
+−+ oGICp ntpmsmms+t+tOp
p
N=T )()(( 1
}).logmax{ 42 oGB ntOp;pMt+ (19)
В большинстве практических задач временем обмена It и факторизации диа-
гонального блока CtO1 можно пренебречь, так как эти времена, как правило, на-
много меньше (при естественном условии )pM > остальных слагаемых в фор-
мулах (18) и (19). При таком предположении оценки ускорения и эффективности
рассматриваемого гибридного алгоритма имеют следующий вид :)log( 2 pMpCB =
.
;,
;,
)2(
log
2
/11
4
4
1
3
2
1
p
S
E
tCntOp
tCntO
t
tn
Ms
pp
M
pMp
T
T
S
p
p
BBoG
BBoG
G
Bo
p
p
=
⎪
⎪
⎩
⎪⎪
⎨
⎧
≥
<⎟
⎟
⎠
⎞
⎜
⎜
⎝
⎛
+
+
+
+
−
≅=
−
(20)
Оценки (20) свидетельствуют, что эффективность алгоритма не зависит от
порядка матрицы, а при определенных условиях определяется количеством нену-
левых блоков в строке блоков и их размером. На эффективность алгоритма также
влияет отношение времени пересылки одного блока между двумя GPU Bt и зна-
чения .oG nt Также следует отметить, что хотя время факторизации диагональ-
ных блоков существенно не влияет на эффективность распараллеливания, его
сокращение позволяет уменьшить общее время разложения матрицы системы
(учитывая, что эти факторизации выполняются в последовательном режиме). На-
пример, можно перенести эти вычисления на GPU.
Экспериментальное исследование гибридных алгоритмов. Описанный
гибридный алгоритм решения СЛАУ с ленточной симметричной матрицей
реализован в программной среде MPI на компьютерах с многоядерными и
графическими процессорами — кластерах семейства СКИТ Института кибер-
нетики им. В.М. Глушкова НАН Украины [13] и семейства Инпарком (совместная
разработка Института кибернетики и ГНПП «Электронмаш» [14]). Для организа-
ции вычислений на GPU использована технология CUDA. При реализации гиб-
ридного алгоритма для матрично-векторных и матрично-матричных операций ис-
пользованы функции библиотек Intel MKL (для вычислений на CPU) и CUBLAS
(для GPU).
Для лучшего использования вычислительных ресурсов гибридного компью-
тера необходимо асинхронное выполнение некоторых операций, при этом для вы-
числений на GPU применяются дополнительные потоки CUDA. Это позволяет
выполнять одновременно операцию синхронизации (viii), разложение на CPU
диагонального блока (i) и 2s-ранговую модификацию по формуле (17) на GPU (iv).
Исследования производительности предложенного гибридного алгоритма
LLT-разложения ленточной симметричной матрицы в зависимости от величины
выбранного блока и количества GPU проводились для двух матриц СЛАУ со сле-
дующими параметрами: первая — порядок n = 1151112, полуширина ленты
m = 5357, вторая — n = 1438890, m = 6693. Результаты (производительность в
Gflops — 109 арифметических операций с плавающей запятой в секунду) пред-
78 ISSN 0572-2691
ставлены в виде диаграмм (А — для первой матрицы, Б — второй) на рис. 2. Эти
исследования позволяют экспериментально определить оптимальный размер бло-
ка. Как видно из приведенных диаграмм, для матриц с указанными выше пара-
метрами оптимальным является размер блока s = 192, так как увеличение размера
дает или незначительный рост, или уменьшение производительности, в то же
время увеличение размера блока с 96 до 192 позволило увеличить производитель-
ность в 1,5–2 раза.
100
0
200
300
400
500
600
96 128 192 256 512 640
2 GPU
3GPU
Размер блока
G
flo
ps
А
100
0
200
300
400
500
600
96 128 192 256 512 640
2 GPU
3GPU
Размер блока
Б
700
G
flo
ps
Рис. 2
Также проведено сравнение предложенного гибридного алгоритма и парал-
лельных алгоритмов для многоядерных компьютеров, в частности алгоритма,
программно реализованного в библиотеке Intel MKL. На рис. 3 представлены ре-
зультаты этого сравнения — ускорение гибридного алгоритма по отношению к
используемому библиотекой Intel MKL в зависимости от количества используе-
мых ядер CPU. Результаты получены при решении СЛАУ порядка 1151112 (по-
луширина ленты матрицы — 5357) и свидетельствуют о высокой эффективности
предложенного алгоритма.
10
0
20
30
40
50
60
2 4 6 8 16
2 GPU
3GPU
Количество
У
ск
ор
ен
ие
Рис. 3
В заключение этого раздела следует отметить, что общая схема предложен-
ных алгоритмов (разбиение на блоки, распределение между процессорными уст-
ройствами, обработка только блоков, соответствующих ненулевым блокам мат-
рицы разложения) может быть успешно применена и для решения СЛАУ с дру-
гими структурами разреженных матриц — профильной, блочно-диагональной с
окаймлением и т.п.
3. Численное моделирование
Для исследования возможностей предложенного параллельного алгоритма
решения СЛАУ решен ряд сложных практических задач проектирования — чис-
ленные исследования ряда объектов, в том числе 27-этажного здания (рис. 4).
В таблице приведены результаты исследований ускорений, получаемых при
решении на компьютере гибридной архитектуры СЛАУ и возникающих при ста-
Международный научно-технический журнал
«Проблемы управления и информатики», 2017, № 4 79
тическом анализе прочности различных строительных объектов и отдельных кон-
струкций. При дискретизации элементов конструкции использованы трехмерные
конечные элементы (ЗD КЭ).
Таблица
Время решения СЛАУ
Объект
Порядок
матрицы
жесткости
Полуширина
ленты
матрицы
Лира–Сапр*
(мин.)
Инпарком_G**
(мин.)
Здание, 12 эт. 189956 22 585 1 0,274
Фундамент 283031 19 530 2 0,384
Балка, 3D КЭ 300000 335 1 0,048
Здание, 27 эт. 661590 34 242 2 0,948
Плита, 3D КЭ 666000 1 004 1 0,138
Плита, 3D КЭ 1000332 1 004 2 0,207
Балка, 3D КЭ 1200000 1 265 16 0,371
Примечание:
*Лира–Сапр [2] — многопоточная версия 2015 г., GPU не использовался.
**Инпарком_G [14] — использован 1 GPU, гибридный алгоритм.
На рис. 4 представлена 3D-модель МКЭ 27-этажного здания. Рассмотрим
общую характеристику и входные данные численной моде-
ли пространственных несущих конструкций этого здания:
• сваи основания моделируются одноузловыми ко-
нечными элементами (КЭ): Rz = 50990 т/м, Rx = 50990 т/м,
Ry = 50990 т/м;
• фундаментная плита на отметке h = – 4 м моделирует-
ся КЭ оболочки; модуль деформации бетона Е = 3,059×
×106 т/м2; коэффициент Пуассона ν = 0,25; толщина плиты
H = 150 см;
• колонны моделируются стержневыми КЭ; Е = 3,059×
×106 т/м2; ν = 0,25; ширина и высота сечения колонн ( b×h )
25×90 см, 25×120 см.
• перекрытия моделируются КЭ оболочки; Е = 3,059×
×106 т/м2; ν = 0,25; толщина плит перекрытия Н = 18 см.
Рассмотрено шесть вариантов нагружений:
(1) собственный вес 27-этажного здания;
(2) полы, перегородки, кровля;
(3) полезные, снег;
(4) трапециевидная нагрузка от грунта на фундаментно-подвальную часть здания;
(5) ветер по направлению X;
(6) ветер по направлению Y.
Таким образом, разрешающая СЛАУ (11) имеет шесть правых частей, т.е.
q = 6.
Некоторые результаты численного анализа прочности высотного здания
представлены детальными картинами распределения перемещений по Y, Z (рис. 5)
и усилий N (рис. 6) для рассматриваемых вариантов нагружений.
Полученные результаты численных исследований поведения строительных
конструкций при статических нагружениях показали многократное сокращение
времени решения СЛАУ с ленточными симметричными положительно-опреде-
ленными матрицами на многопроцессорных (многоядерных) компьютерах с гра-
фическими ускорителями при использовании предложенных гибридных алгоритмов.
Рис. 4
80 ISSN 0572-2691
X
Z Y
– 20,3
– 17,8
– 15,2
– 12,7
– 10,2
– 7,62
– 5,08
– 0,203
0
– 2,54
X
Z Y
– 0,328
– 0,287
– 0,246
– 0,207
–0,164
– 0,123
– 0,0819
0,0810
0,845
0,041
– 0,041
– 0,00844
0,00844
Z
Y
0
0,131
1,63
3,27
4,9
6,53
8,16
9,8
13,1
11,4
Рис. 5
X
Z Y
– 275
– 241
– 206
– 172
– 137
– 103
– 68,7
– 34,4
– 0,134
0,134
134
X
Z Y
– 116
– 87,1
– 72,6
– 58,1
– 43,6
– 29
– 68,7
– 14,5
– 0,0618
0,0618
6,19
X
Z Y
– 32,4
– 28,3
–24,2
– 20,12
– 16,2
– 12,1
– 4,04
– 8,08
– 0,287
0,287
4,04
8,08
12,01
16,2
20,2
24,2
28,3
28,7
Рис. 6
Заключение
Общая схема параллельных алгоритмов (15)–(17), а также ее реализация (разбие-
ние на блоки, распределение между процессорными устройствами, обработка блоков,
соответствующих ненулевым блокам матрицы разложения) могут быть успешно
применены и для решения СЛАУ с другими структурами разреженных матриц —
профильной, блочно-диагональной с окаймлением и т.п. (см., например, [15]).
Оптимизация вычислений для компьютеров гибридной архитектуры лежит,
очевидно, в такой модификации алгоритмов из [8], которая позволяет более полно
учитывать разреженную структуру матрицы жесткости и ее LLT-разложения.
Проведенная апробация предложенных гибридных алгоритмов на решении
практических задач свидетельствует о перспективности этого направления в раз-
витии алгоритмического и программного обеспечения для математического моде-
лирования строительных конструкций.
Особенно существенный эффект ускорения можно ожидать при решении прак-
тических задач проектирования ответственных особо сложных пространственных
конструкций при статических и динамических нагружениях, в том числе при иссле-
довании жизненного цикла высотных сооружений с учетом нелинейных свойств ма-
териалов конструкций и грунтов оснований с использованием вышеизложенных пря-
мых методов решения СЛАУ (с количеством неизвестных 106–109).
Международный научно-технический журнал
«Проблемы управления и информатики», 2017, № 4 81
Использование приведенных высокопроизводительных параллельных алго-
ритмов позволяет обеспечить высокое качество проектирования и безопасность
возведения, эксплуатации и реконструкции особо сложных строительных объектов.
А.Ю. Баранов, О.В. Попов, Я.О. Слободян, О.М. Хіміч
МАТЕМАТИЧНЕ МОДЕЛЮВАННЯ МІЦНОСТІ
БУДІВЕЛЬНИХ КОНСТРУКЦІЙ НА ГІБРИДНИХ
ОБЧИСЛЮВАЛЬНИХ СИСТЕМАХ
Розглядаються алгоритми розв’язування на комп’ютерах гібридної архітектури
(на багатоядерних процесорах і графічних прискорювачах) систем лінійних ал-
гебраїчних рівнянь зі стрічковими симетричними додатно-визначеними матри-
цями, які виникають при чисельному моделюванні з використанням методу
скінченних елементів, просторових будівельних конструкцій. Наведено результа-
ти апробації запропонованих гібридних алгоритмів на розв’язуванні задач статич-
ного аналізу міцності декількох будівельних об’єктів і окремих конструкцій.
A.Yu. Baranov, A.V. Popov, Ya.E. Slobodyan, A.N. Khimich
MATHEMATICAL MODELING OF STRENGTH
OF BUILDING STRUCTURES FOR HYBRID
COMPUTING SYSTEMS
Algorithms for solving on computers of hybrid architecture (with multi-core proces-
sors and graphics accelerators) the systems of linear algebraic equations with band
symmetrical positive definite matrices, arising in the numerical modeling of 3D building
constructions using the finite element method is considered. Some results of the pro-
posed hybrid algorithms testing for solving the problems of static strength analysis of
construction sites and individual constructions are presented.
1. https://www.top500.org/lists/2016/11/
2. http://www.liraland.ua/company/person.php
3. Тимошенко С.П., Гудьер Дж. Теория упругости. — М. : Наука, 1975. — 575 с.
4. Стренг Г., Фикс Дж. Теория метода конечных элементов. — М. : Мир, 1973. — 349 с.
5. Параллельные алгоритмы решения задач вычислительной математики / А.Н. Химич,
И.Н. Молчанов, А.В. Попов, Т.В. Чистякова, М.Ф. Яковлев. — Киев : Наук. думка, 2008. — 248 с.
6. Парлетт Б. Симметричная проблема собственных значений. — М. : Мир, 1983. — 318 с.
7. Попов А.В., Химич А.Н. Параллельный алгоритм решения систем линейных алгебраических
уравнений с ленточной симметричной матрицей // Компьютерная математика. — 2005. —
№ 2. — С. 52–59.
8. Химич А.Н., Попов А.В., Полянко В.В. Алгоритмы параллельных вычислений для задач ли-
нейной алгебры с матрицами нерегулярной структуры // Кибернетика и системный анализ.
— 2011. — № 6. — С. 159 –174.
9. Хіміч О.М., Баранов А.Ю. Гібридний алгоритм розв’язування лінійних систем із стрічкови-
ми матрицями прямими методами // Компьютерная математика. — 2013. — № 2. — С. 80–87.
10. Численные методы для многопроцессорного вычислительного комплекса ЕС / Михале-
вич В.С., Молчанов И.Н., Сергиенко И.В. и др. / Под ред. И.Н. Молчанова. — М. : ВВИА
им. проф. Н.Е. Жуковского, 1986. — 401 с.
11. Ортега Дж. Введение в параллельные и векторные методы решения линейных систем. —
М. : Мир, 1991. — 368 с.
12. Воеводин В.В., Воеводин Вл.В. Параллельные вычисления. — СПб. : БХВ-Петербург, 2002.
— 608 с.
13. http://icybcluster.org.ua
14. http://www.inparcom.com
15. Хіміч О.М., Сидорук В.А. Плитковий гібридний алгоритм факторизації розріджених блочно-
діагональних матриць з обрамленням // Компьютерная математика. — 2016. — № 1. — С. 72−79.
Получено 29.12.2016
|