Файл: Разработка программного модуля для автоматического выбора решателей систем линейных алгебраических уравнений для прочностного анализа.pdf

Программные продукты и системы / Software & Systems
№ 4 (108), 2014
162
УДК 519.6, 539.3
Дата подачи статьи: 11.04.2014
DOI: 10.15827/0236-235X.108.162-166
РАЗРАБОТКА ПРОГРАММНОГО МОДУЛЯ
ДЛЯ АВТОМАТИЧЕСКОГО ВЫБОРА РЕШАТЕЛЕЙ
СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ
ДЛЯ ПРОЧНОСТНОГО АНАЛИЗА
Н.Е. Стёпин
, аспирант,
N_i_k_i1989@mail.ru
(Московский государственный университет им. М.В. Ломоносова,
Ленинские горы, г. Москва, 119991, Россия)
В работе реализован программный модуль, объединивший в себе различные алгоритмы и методы: прямые и ите- рационные решатели для симметричных и несимметричных матриц системы, различные предобуславливатели в ите- рационных методах, различные способы хранения матрицы в памяти, параллельные вычисления с использованием технологий OpenMP и CUDA. В программном модуле реализован метод решения задач для несжимаемых материа- лов на основе алгоритма Узавы. Для программного модуля разработан и реализован алгоритм оптимального выбора решателя в зависимости от механической постановки задачи, ее размерности и возможностей компьютера. При же- лании пользователь может сам ограничивать некоторые возможности выбора и задавать параметры, влияющие на выбор решателя, или даже указать явно, какой решатель он хочет использовать. По сути программный модуль явля- ется некоторой оболочкой над отдельными решателями, которая принимает матрицу системы, правую часть и неко- торые параметры настройки, а затем в рамках содержащегося в ней алгоритма определяет, какой именно решатель необходимо запускать, настраивает его и приводит матрицу к соответствующему виду (разные решатели могут иметь разные оптимальные форматы хранения для матриц). Проведен ряд численных экспериментов, подтверждаю- щих обоснованность используемых в алгоритме критериев.
Ключевые слова: теория упругости, несжимаемые материалы, метод конечных элементов, алгоритм Узавы,
итерационные методы решения систем линейных алгебраических уравнений.
В процессе решения задач механики деформи- руемого твердого тела при конечных деформациях для сжимаемых и несжимаемых материалов ме- тодом конечных элементов получаются системы нелинейных уравнений. Их решение методом
Ньютона–Рафсона сводится к решению на каждой итерации системы линейных алгебраических урав-
нений (СЛАУ) с разреженной матрицей (размер- ность матрицы зависит от размерности задачи и того, насколько мелкая сетка используется в мето- де конечных элементов) с нерегулярной структу- рой. Решение задачи тем точнее, чем мельче сетка, поэтому актуальна возможность решения системы уравнений максимально возможной размерности при заданных ресурсах ЭВМ [1–3].
В работе реализован программный модуль, объединивший различные алгоритмы и методы: прямые и итерационные решатели для симмет- ричных и несимметричных матриц системы, раз- личные предобуславливатели в итерационных методах и способы хранения матрицы в памяти, параллельные вычисления с использованием тех- нологий OpenMP и CUDA [4–7]. В программном модуле реализован метод решения задач для не- сжимаемых материалов на основе алгоритма Уза- вы [8–12].
Для программного модуля разработан и реали- зован алгоритм оптимального выбора решателя в зависимости от механической постановки задачи, ее размерности и возможностей компьютера [13].
При желании пользователь может сам ограничи- вать некоторые возможности выбора и задавать параметры, влияющие на выбор решателя, или даже указать явно, какой решатель он хочет ис- пользовать. По сути программный комплекс явля- ется некоторой оболочкой над отдельными реша- телями, которая принимает от программы, выпол- няющей механическую задачу (будем называть ее ядром), матрицу системы, правую часть и при же- лании некоторые параметры настройки, а затем в рамках содержащегося в ней алгоритма определя- ет, какой именно решатель необходимо запускать, настраивает его и приводит матрицу к соответст- вующему виду (разные решатели могут иметь разные оптимальные форматы хранения для мат- риц). После окончания работы решателя результат решения возвращается в ядро для дальнейшей обработки. На рисунке 1 схематично показан принцип работы модуля и его взаимодействия с ядром.
Основные учитываемые критерии и их влияние на выбор решателя:
1) если матрица малой размерности (критерий основывается на доступных ресурсах используе- мого компьютера), выбирается решатель, исполь- зующий прямой метод;
2) если задача имеет симметричную матрицу, в качестве итерационного метода целесообразно использовать метод сопряженных градиентов
(CG), в противном случае – бисопряженных гра- диентов (BiCG) или обобщенных минимальных невязок (GMRes);
3) используется наиболее эффективный (бы- стрый, но необязательно самый «сильный») пре- добуславливатель, имеющийся для выбранного решателя;

Программные продукты и системы / Software & Systems
№ 4 (108), 2014
163 4) если компьютер позволяет использовать технологию CUDA, выбирается решатель, исполь- зующий эту технологию;
5) если в задаче используется несжимаемый материал, выбирается метод на основе алгоритма
Узавы.
Проведен ряд численных экспериментов, под- тверждающих обоснованность приведенных выше критериев.
Если по какой-то причине после выбора реша- теля задача не была решена (например, если вы- бранным методом и с выбранным предобуславли- вателем итерационный процесс не сошелся), мо- жет быть запущен другой вариант (например, с более «сильным» предобуславливателем).
Программный модуль уже был использован в программном комплексе прочностного анализа
CAE FIDESYS [14]. Пользовательский интерфейс для задания настроек программного модуля соз- дан так, что некоторые параметры (например, ука- зывающий на решение задачи с несжимаемым ма- териалом или на симметричность задачи) скрыты от пользователя и определяются программно в за- висимости от постановки задачи или параметров, заданных пользователем для материала.
Выбор решателя 1. Всегда, если есть возмож- ность, стоит использовать прямые методы. Их преимущество в том, что они дают точное реше- ние (итерационные методы будут давать прибли- женное), конечно, с учетом машинной точности.
Другим преимуществом прямых методов является то, что они находят решение за один проход, по- этому очень быстры. Но у них есть и серьезный недостаток: требуется большой объем памяти для проведения численной факторизации, поэтому при больших размерах матрицы их использование ста- новится невозможным. Если, помимо прочего, из- вестно, что матрица системы симметрична, то стоит использовать это для оптимизации (напри- мер, можно экономить память, храня только верх- нюю часть матрицы) [4, 5].
Выбор решателя 2. В большинстве реальных
3D-задач матрица имеет слишком большую раз- мерность и оперативной памяти не хватает для того, чтобы можно было пользоваться прямыми методами. В таких случаях используются итера- ционные методы (проекционные методы на со- пряженных подпространствах пространства Кры- лова, основанные на минимизации невязки).
Здесь симметричность может служить уже не просто для оптимизации. Существуют специаль- ные методы только для симметричных СЛАУ, которые являются более стабильными и эффек- тивными, нежели методы для несимметричных систем. В симметричном случае хорошо зареко- мендовал себя метод сопряженных градиентов
(Conjugate Gradient). Для несимметричных задач можно использовать его модификацию (Biconju- gate Gradient Stabilized) или модификацию метода обобщенных минимальных невязок (Flexible Gene- ralized Minimal Residual) [4–7].
Выбор решателя 3. При сложной геометрии задачи и использовании неструктурированных расчетных сеток портрет матрицы может быть не- регулярным (портрет – множество пар индексов, соответствующих ненулевым элементам), а число ее обусловленности – очень большим (от этой величины зависит, насколько хорошо сходится решение задачи итерационными методами; чем больше эта величина, тем решение сходится ху- же).
Для уменьшения числа обусловленности в программном модуле реализованы различные возможности использования предобуславливате- лей [4–7]:
– без предобуславливателя;
– предобуславливатель Якоби (диагональ- ный);
– SSOR (symmetric successive overrelaxation – симметричный метод последовательной верхней релаксации);
– AINV – параметризуемый, поэтому реали- зовано несколько видов под разные наборы пара- метров;
– неполное LU-разложение и его разновидно- сти: ILU(0), ILU(2), ILU
T
;
– AMG (алгебраический многосеточный).
Стоит отметить, что некоторые из них подхо- дят только для симметричных случаев, так как ис- пользование такого предобуславливателя может нарушить симметрию СЛАУ. Некоторые предобу- славливатели реализованы только для CPU, а дру- гие эффективны только для GPU.
Обратим также внимание на то, что использо- вание предобуславливателя вынуждает нас вы- полнять некоторые дополнительные действия, связанные с созданием предобуславливателя и домножением на него. Чтобы убедиться в оправ-
Интерфейс
Решатель #1
Решатель #N
Ядро
…
Рис. 1. Схема передачи СЛАУ решателям
и получения результата
Fig. 1. A
system of linear equations transfer circuit
and results obtaining

Программные продукты и системы / Software & Systems
№ 4 (108), 2014
164 данности этих действий и в том, что уменьшение количества итераций перекрывает все дополни- тельные затраты на использование предобуслав- ливателя, были проведены некоторые численные эксперименты. Как правило, чем сложнее в по- строении предобуславливатель, тем больше вре- мени требуют эти дополнительные действия. Но есть и обратная сторона: чем сложнее предобу- славливатель, тем сильнее он уменьшает число обусловленности матрицы, а значит, дает лучшую сходимость. Необходимо искать компромисc.
Выбор решателя 4. Итерационные методы хо- рошо поддаются распараллеливанию. Поэтому, если позволяет архитектура компьютера, очень важно использовать такую возможность. Много- поточность на центральном процессоре (CPU) есть уже практически на любом современном компьютере. Но если компьютер позволяет осу- ществлять распараллеливание счета с использова- нием видеокарты (GPU), это может давать гораздо более заметный прирост производительности и не стоит упускать такую возможность. Однако это касается не всех видеокарт. Одними из наиболее популярных являются видеокарты NVIDIA, в современных моделях которых имеется поддерж- ка технологии CUDA. Именно эта технология и позволяет перенести решение СЛАУ с CPU на
GPU [13].
На рисунке 2 сравниваются результаты расчета на центральном процессоре (Host) и на графиче- ском устройстве (Device) на двух различных пер- сональных компьютерах (AMD Athlon 64 X2 Dual
Core с графическим устройством NVIDIA GF GTX
260 и Intel Pentium Dual Core E6500 с графическим устройством NVIDIA GF GTX 295). Расчет делал- ся для задачи теории упругости с симметричной матрицей системы линейных алгебраических уравнений, содержащей 300 тыс. строк; расчеты производились методом сопряженных градиентов без предобуславливателя и с предобуславливате- лем Якоби (на рисунке обозначено как Jacobi). Как видно из рисунка, использование графического устройства дает заметно большее преимущество, чем увеличение количества процессоров. Предо- буславливание тоже дает очень заметный прирост.
А если их совмещать, время расчета значительно уменьшается (в 40–75 раз).
Выбор решателя 5. В программном модуле реализован особый метод решения задач для не- сжимаемых материалов на основе алгоритма Уза- вы. Дело в том, что получаемые в этих задачах системы уравнений являются системами с седло- вой точкой (собственные значения матрицы будут разных знаков) и их матрицы имеют особый вид: в правом нижнем углу они содержат блок из нулей.
Описанные выше общепринятые и наиболее по- пулярные методы решения СЛАУ на таких зада- чах чаще всего не сходятся из-за плохой обуслов- ленности этих матриц. Алгоритм Узавы позволяет с помощью релаксационного процесса свести ре- шение такой СЛАУ к решению на каждой итера- ции СЛАУ меньшей размерности, уже не имею- щей такого блока нулей, поэтому можно будет ис- пользовать классические методы решения СЛАУ с разреженной матрицей.
Известно несколько вариантов алгоритма Уза- вы, основанных на различных релаксационных методах решения СЛАУ. В данном программном модуле реализованы варианты, основанные на ме- тоде наискорейшего градиентного спуска, методе сопряженных градиентов (двух- и трехслойная схемы) и трехслойном методе сопряженных невя- зок. Расчетные формулы для этих вариантов мето- да Узавы записываются по аналогии с расчетными формулами соответствующих итерационных ме- тодов [8–12].
Как и классические итерационные методы ре- шения СЛАУ, для разных механических задач мо- гут подходить разные варианты метода (какие-то требуют больше итераций, а какие-то могут не сходиться), но, тем не менее, в большинстве слу- чаев наиболее стабильным и эффективным по ко- личеству итераций и времени расчета оказался ва- риант на основе метода сопряженных невязок. Это можно увидеть из приведенного далее примера.
Сравним полученные результаты для четырех матриц различной размерности:
1) 30 402 строки, из них 20 402 приходятся на главный блок (матрица А);
2) 120 802 строки, из них 80 802 приходятся на главный блок;
3) 246 534 строки, из них 164 738 приходятся на главный блок;
4) 481 602 строки, из них 321 602 приходятся на главный блок.
Рис. 2. Производительность при вычислениях
на графическом (GPU) и центральном (CPU)
процессорах с использованием предобуславливателя
Якоби и без него на примере СЛАУ
с 300 тыс. неизвестных
Fig. 2. Computing productivity for GPU and CPU
using Jacobi precondition and without it, case study
is a system of linear equations in 3000 k unknowns
0 50 100 150 200 250 300 350 400 450
Host
Device
Host + Jacobi
Device + Jacobi
277,3 49,1 38,2 7,2
Время
работы, сек.
AMD Athlon 64 x2 dual Core, GF GTX 260
Intel Pentium Dual Core E6500, GF GTX 295

Программные продукты и системы / Software & Systems
№ 4 (108), 2014
165
Критерием окончания расчета было уменьше- ние нормы невязки в 100 раз по сравнению с нор- мой начальной невязки. Для решения СЛАУ с матрицей А на каждой итерации метода Узавы ис- пользовались прямые методы.
На графике (рис. 3) по вертикальной оси отме- чено время (в секундах), необходимое для реше- ния системы тем или иным вариантом алгоритма, а по горизонтальной оси отмечены матрицы, для которых приведены данные. Разными символами на линиях отмечены разные методы.
Отметим, что количество итераций алгоритма
Узавы почти не зависит от размерности матрицы, но, тем не менее, как можно видеть из графика, время расчета увеличивается с размером матрицы
(матрицы на графике упорядочены по возраста- нию их размерности). Это связано с ростом вре- мени, потраченного на одну итерацию алгоритма и решение СЛАУ с матрицей A, размерность ко- торой тоже растет. Также можно видеть, что наи- более эффективным и стабильным является метод
Узавы, основанный на методе сопряженных невя- зок (CRes).
В заключение отметим, что в статье сравнива- ется эффективность различных методов решения
СЛАУ в задачах теории упругости сжимаемых и несжимаемых материалов при малых и конечных деформациях. Составлен алгоритм выбора наибо- лее эффективного метода с учетом размерности задачи, особенностей ее механической постановки и возможностей компьютера, на котором произво- дится расчет. Алгоритм реализован в виде про- граммного модуля, включающего в себя большой набор методов. Этот модуль может быть исполь- зован в качестве составной части прочностных па- кетов и при решении реальных задач механики в программных комплексах [14].
Литература
1. Левин В.А., Калинин В.В., Зингерман К.М., Верши- нин А.В. Развитие дефектов при конечных деформациях. Ком- пьютерное и физическое моделирование. М.: Физматлит, 2007.
392 с.
2. Zienkiewicz O.C., Taylor R.L. The finite element method.
The basis. Oxford: Butterworth-Heineman, 2000, vol. 1, 691 p.
3. Levin V.A., Vershinin A.V. Non-stationary plane problem of thesuccessive origination of stress concentrators in a loaded bo- dy. Finite deformations and their superposition. Communications in
Numerical Methods in Engineering, 2008, vol. 24, pp. 2229–2239.
4. Бахвалов Н.С., Жидков Н.П., Кобельков Г.М. Числен- ные методы. М.: Наука, 2003. 632 с.
5. Богачев К.Ю. Методы решения линейных систем и на- хождения собственных значений. М.: Изд-во мех.-мат. ф-та
МГУ, 1998.
6. Баландин М.Ю., Шурина Э.П. Методы решения СЛАУ большой размерности. Н.: Изд-во НГТУ, 2000.
7. Saad Y. Iterative Methods for Sparse Linear Systems,
Second edition, SIAM, Philadelphia, 2003.
8. Быченков Ю.В., Чижонков Е.В. Итерационные методы решения седловых задач. М.: БИНОМ. Лаборатория знаний,
2010. 349 с.
9. Чижонков Е.В. Релаксационные методы решения сед- ловых задач. М.: Изд-во ИВМ РАН, 2002. 238 с.
10. Ладыженская О.А. Математические вопросы динамики вязкой несжимаемой жидкости. М.: Наука, 1970.
11. Быченков Ю.В. Оптимизация предобусловленных ме- тодов для седловых задач // Докл. Акад. наук. 2002. Т. 384.
№ 4. С. 439–441.
12. Степин Н.Е., Левин В.А., Зингерман К.М., Верши- нин А.В. Сравнительный анализ различных вариантов алго- ритма Узавы в задачах упругости для несжимаемых материа- лов // Вестн. ТвГУ: Сер. Прикладная математика. № 3 (26).
С. 29–34.
13. Левин В.А., Вершинин А.В., Траченко А.В., Проко- пенко А.С., Степин Н.Е. Некоторые результаты использования технологии CUDA при решении СЛАУ для задач прочности при перераспределении конечных деформаций / В кн.: Совре- менные проблемы математики, механики, информатики: матер.
10-й Междунар. конф. Тула, 2009.
14.
Levin V.A., Zingerman K.M., Vershinin A.V., Nikifo- rov I.V. CAE FIDESYS for strength analysis at large strains and their redistribution. 10th Word Congress on Computational
Mechanics, 8
–
13 July 2012, Sao Paulo, Brazil, 2012, 19579, p. 323.
DOI: 10.15827/0236-235X.108.162-166
Received 11.04.2014
DEVELOPING A SOFTWARE MODULE FOR AN AUTOMATIC CHOICE
OF SPARSE SOLVER IN STRENGTH ANALYSIS
Stepin N.E., Postgraduate Student, N_i_k_i1989@mail.ru
(Lomonosov Moscow State University, Leninskie Gory, Moscow, 119991, Russian Federation)
Аbstract. The paper presents a software module which includes a variety of different algorithms and methods: direct and iterative sparse solvers for symmetric and non-symmetric matrices, various preconditioners, different formats for storing sparse matrices in RAM, HPC technologies based on OpenMP and CUDA. A method for solving incompressible materials problems based on Uzawa algorithm was implemented in this software module. An algorithm for an optimal sparse solver choice depending on mechanical problem, its dimensions and available hardware resources was designed and implemented.
Рис. 3. Время решения системы различными
вариантами алгоритма Узавы для различной
размерности матриц системы уравнений
Fig. 3. System solving time using different
Uzawa's methods and different matrices
0 2
4 6
8 10 12
Матрица 1
Матрица 2
Матрица 3
Матрица 4
Mres
StDes
CG3
Cres
CG2

Программные продукты и системы / Software & Systems
№ 4 (108), 2014
166
A user can limit some options and set parameters that influence a choice of a sparse solver or directly specify which solver to use. In fact, the software module is some kind of a wrapper over individual sparse solvers which takes a matrix and a right hand side as an input as well as some solver’s settings, and then decides which solver to run based on the algorithm, config- ures it and converts a matrix to a corresponding format (different solvers may have different optimal storage formats for a sparse matrix). A series of numerical experiments that confirms the validity of the criteria used in the algorithm was per- formed by the author.
Keywords: theory of elasticity, incompressible materials, finite element method, Uzawa algorithm, iterative methods for sparse systems of linear algebraic equations.
References
1. Levin V.A., Kalinin V.V., Zingerman K.M., Vershinin A.V. Razvitie defektov pri konechnykh deformatsiyakh.
Kompyuternoe i fizicheskoe modelirovanie [Defects Development at Finite Deformations. Computational and Physical Mod- eling]. Moscow, Fyzmatlit Publ., 2007, 392 p.
2. Zienkiewicz O.C., Taylor R.L. The finite element method. Vol. 1. The basis. Oxford, Butterworth-Heineman Publ.,
2000, 691 p.
3. Levin V.A., Vershinin A.V. Non-stationary plane problem of the successive origination of stress concentrators in a loaded body. Finite deformations and their superposition. Communications in Numerical Methods in Engineering. 2008, no. 24, pp. 2229–2239.
4. Bakhvalov N.S., Zhidkov N.P., Kobelkov G.M. Chislennye metody [Numerical Methods]. Moscow, Science Publ.,
2003, 632 p.
5. Bogachov K.U. Methody resheniya lineynykh system i nakhozhdeniya sobstvennykh znacheniy [Linear System Solv- ing Methods and Eigenvalues Finding]. Moscow, MSU Mechanic and Mathematic Department Publ., 1998.
6. Balandin M.Yu., Shurina E.P. Metody resheniya SLAU bolshoy razmernosti [Solving Methods of Large Dimension
Linear Algebraic Equation Systems]. Novosibirsk, NNSTU Publ., 2000.
7. Saad Y. Iterative methods for sparse linear systems. 2nd edition, SIAM, Philadelphia, 2003.
8. Bychenkov Yu.V., Chizhonkov E.V. Iteratsionnye metody resheniya sedlovykh zadach [Iterative Methods for Solv- ing Saddle Point Problems]. Moscow, BINOM Knowledge Laboratory Publ., 2010, 349 p.
9. Chizhonkov E.V. Relaksatsionnye metody resheniya sedlovykh zadach [Relaxation Methods for Solving Saddle Point
Problems]. Moscow, 2002.
10. Ladyzhenskaja O.A. Matematicheskie voprosy dinamiki vyazkoy neszhimaemoy zhidkosti [The Mathematical Theory of Viscous Incompressible Fluid]. Moscow, Science Publ., 1970.
11. Bychenkov Yu.V. Preconditioner methods optimization for saddle point problems. Dokl. Akad. Nauk [Reports of the
Academy of Sciences]. 2002, vol. 384, no. 4. pp. 439–441 (in Russ.).
12. Stepin N.E., Levin V.A., Zingerman K.M., Vershinin A.V. Comparative analysis of Uzava algorithm various options for elasticity problems of incompressible materials. Vestn. TvGU: Ser. Prikladnaya matematika [Bulletin of the Tver State
Univ. Applied Mathematics Series], no. 3 (26), pp. 29–34 (in Russ.).
13. Levin V.A., Vershinin A.V., Trachenko A.V., Prokopenko A.S., Stepin N.E. Some results of using CUDA technolo- gy for linear algebraic equation systems solving for problems of strength in the redistribution of finite deformations. Mater.
10 Mezhdunar. konf. “Sovremennye problemy matematiki, mekhaniki, informatiki” [Proc. 10th Int. Conf. Modern Problems of Mathematics, Mechanics, Informatics]. Tula, 2009.
14.
Levin V.A., Zingerman K.M., Vershinin A.V., Nikiforov I.V. CAE FIDESYS for strength analysis at large strains and their redistribution. Proc. 10th Word Congress on Computational Mechanics. 8–13 July 2012, Sao Paulo, Brazil, 19579, p. 323.
Вниманию авторов!
Редакция международного журнала «Программные продукты и системы» принимает к рассмотре- нию оригинальные, ранее нигде не опубликованные материалы, соответствующие тематике журнала
(специализация 05.13.ХХ – Информатика, вычислительная техника и управление).
Представление чужой работы как авторской, копирование или перефразирование ее существенных частей (без указания авторства и соответствующих библиографических ссылок), заявление собственных прав на результаты чужих исследований неэтичны и неприемлемы. В случае установления редакцией факта плагиата статья не рассматривается, а автор заносится в негласный редакционный список как утра- тивший доверие.
В течение месяца с момента поступления рукопись проходит обязательное слепое рецензирование
(то есть рецензент не знает, кто является автором статьи, а автор – кто рецензент). Эксперты назначаются редакционной коллегией журнала или рабочим советом. Заключение рецензента сообщается автору.
При подаче статьи на рассмотрение просьба учитывать формальные редакционные требования, подробно изложеннные на сайте журнала www.swsys.ru.
Редакция