CC BY
108
Программные продукты и системы /Software & Systems
№ 4 (108), 2014
УДК 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) используется наиболее эффективный (быстрый, но необязательно самый «сильный») пре-добуславливатель, имеющийся для выбранного решателя;
162
Программные продукты и системы /Software & Systems
№ 4 (108), 2014
f \
Ядро
/----------------\
Решатель #N
Рис. 1. Схема передачи СЛАУ решателям и получения результата
Fig. 1. A system of linear equations transfer circuit and results obtaining
4) если компьютер позволяет использовать технологию CUDA, выбирается решатель, использующий эту технологию;
5) если в задаче используется несжимаемый материал, выбирается метод на основе алгоритма Узавы.
Проведен ряд численных экспериментов, подтверждающих обоснованность приведенных выше критериев.
Если по какой-то причине после выбора решателя задача не была решена (например, если выбранным методом и с выбранным предобуславли-вателем итерационный процесс не сошелся), может быть запущен другой вариант (например, с более «сильным» предобуславливателем).
Программный модуль уже был использован в программном комплексе прочностного анализа CAE FIDESYS [14]. Пользовательский интерфейс для задания настроек программного модуля создан так, что некоторые параметры (например, указывающий на решение задачи с несжимаемым материалом или на симметричность задачи) скрыты от пользователя и определяются программно в зависимости от постановки задачи или параметров, заданных пользователем для материала.
Выбор решателя 1. Всегда, если есть возможность, стоит использовать прямые методы. Их преимущество в том, что они дают точное решение (итерационные методы будут давать приближенное), конечно, с учетом машинной точности. Другим преимуществом прямых методов является то, что они находят решение за один проход, поэтому очень быстры. Но у них есть и серьезный недостаток: требуется большой объем памяти для проведения численной факторизации, поэтому при больших размерах матрицы их использование становится невозможным. Если, помимо прочего, известно, что матрица системы симметрична, то стоит использовать это для оптимизации (напри-
мер, можно экономить память, храня только верхнюю часть матрицы) [4, 5].
Выбор решателя 2. В большинстве реальных SD-задач матрица имеет слишком большую размерность и оперативной памяти не хватает для того, чтобы можно было пользоваться прямыми методами. В таких случаях используются итерационные методы (проекционные методы на сопряженных подпространствах пространства Крылова, основанные на минимизации невязки).
Здесь симметричность может служить уже не просто для оптимизации. Существуют специальные методы только для симметричных СЛАУ, которые являются более стабильными и эффективными, нежели методы для несимметричных систем. В симметричном случае хорошо зарекомендовал себя метод сопряженных градиентов (Conjugate Gradient). Для несимметричных задач можно использовать его модификацию (Biconjugate Gradient Stabilized) или модификацию метода обобщенных минимальных невязок (Flexible Generalized Minimal Residual) [4-7].
Выбор решателя 3. При сложной геометрии задачи и использовании неструктурированных расчетных сеток портрет матрицы может быть нерегулярным (портрет - множество пар индексов, соответствующих ненулевым элементам), а число ее обусловленности - очень большим (от этой величины зависит, насколько хорошо сходится решение задачи итерационными методами; чем больше эта величина, тем решение сходится хуже).
Для уменьшения числа обусловленности в программном модуле реализованы различные возможности использования предобуславливате-лей [4-7]:
- без предобуславливателя;
- предобуславливатель Якоби (диагональный);
- SSOR (symmetric successive overrelaxation -симметричный метод последовательной верхней релаксации);
- AINV - параметризуемый, поэтому реализовано несколько видов под разные наборы параметров;
- неполное LU-разложение и его разновидности: ILU(0), ILU(2), ILUT;
- AMG (алгебраический многосеточный).
Стоит отметить, что некоторые из них подходят только для симметричных случаев, так как использование такого предобуславливателя может нарушить симметрию СЛАУ. Некоторые предобу-славливатели реализованы только для CPU, а другие эффективны только для GPU.
Обратим также внимание на то, что использование предобуславливателя вынуждает нас выполнять некоторые дополнительные действия, связанные с созданием предобуславливателя и домножением на него. Чтобы убедиться в оправ-
163
Программные продукты и системы /Software & Systems
№ 4 (108), 2014
данности этих действий и в том, что уменьшение количества итераций перекрывает все дополнительные затраты на использование предобуслав-ливателя, были проведены некоторые численные эксперименты. Как правило, чем сложнее в построении предобуславливатель, тем больше времени требуют эти дополнительные действия. Но есть и обратная сторона: чем сложнее предобу-славливатель, тем сильнее он уменьшает число обусловленности матрицы, а значит, дает лучшую сходимость. Необходимо искать компромисс
Выбор решателя 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). Расчет делался для задачи теории упругости с симметричной
Время работы, сек.
400 ■ 350 ■ 300 ■ 250 ■ 200 ■ 150 ■ 100 ■ 50 ■ 0 ■
277,3
Г
49,1 38,2 7,2
Host Device Host + Jacobi Device + Jacobi ■ AMD Athlon 64 x2 dual Core, GF GTX 260 □ Intel Pentium Dual Core E6500, GF GTX 295
Рис. 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
матрицей системы линейных алгебраических уравнений, содержащей 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 приходятся на главный блок.
164
Программные продукты и системы /Software & Systems
№ 4 (108), 2014
Критерием окончания расчета было уменьшение нормы невязки в 100 раз по сравнению с нормой начальной невязки. Для решения СЛАУ с матрицей А на каждой итерации метода Узавы использовались прямые методы.
На графике (рис. 3) по вертикальной оси отмечено время (в секундах), необходимое для решения системы тем или иным вариантом алгоритма, а по горизонтальной оси отмечены матрицы, для которых приведены данные. Разными символами на линиях отмечены разные методы.
Рис. 3. Время решения системы различными вариантами алгоритма Узавы для различной размерности матриц системы уравнений
Fig. 3. System solving time using different Uzawa's methods and different matrices
Отметим, что количество итераций алгоритма Узавы почти не зависит от размерности матрицы, но, тем не менее, как можно видеть из графика, время расчета увеличивается с размером матрицы (матрицы на графике упорядочены по возрастанию их размерности). Это связано с ростом времени, потраченного на одну итерацию алгоритма и решение СЛАУ с матрицей 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 body. 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. Степин Н.Е., Левин В.А., Зингерман К.М., Вершинин А.В. Сравнительный анализ различных вариантов алгоритма Узавы в задачах упругости для несжимаемых материалов // Вестн. ТвГУ: Сер. Прикладная математика. N° 3 (26). С. 29-34.
13. Левин В.А., Вершинин А.В., Траченко А.В., Прокопенко А.С., Степин Н.Е. Некоторые результаты использования технологии CUDA при решении СЛАУ для задач прочности при перераспределении конечных деформаций / В кн.: Современные проблемы математики, механики, информатики: матер. 10-й Междунар. конф. Тула, 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. 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)
Abstract. 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.
165
Программные продукты и системы /Software & Systems
№ 4 (108), 2014
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, configures 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 performed 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 Modeling]. 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 Solving 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 Solving 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 technology 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.
Редакция
166
