Analysis of variants of Cholesky method for solution of big sparse systems of linear algebraic equations

Abstract

Algorithms of solution of big sparse systems of linear algebraic equations based on Cholesky technique is released in program complex of analysis of stress-strain condition under cyclic elastoplastic loading by finite element method. Comparison of memory and computational resource efficiency of algorithms is carried out.

Full Text

Анализ напряженно-деформированного состояния конструкций методом конечных элементов (МКЭ) с математической точки зрения сводится к решению системы линейных алгебраических уравнений (СЛАУ) [1, 2]. Из-за большой размерности получаемой СЛАУ, зависящей от количества степеней свободы решаемой задачи, именно хранение ее матрицы и собственно решение самой СЛАУ является наиболее ресурсоемкой в смысле объема требуе- мой памяти и процессорного времени частью решения задачи МКЭ. Ранее в модуле решателя МКЭ «Cycle2D» [3] использовался ленточный метод хранения и решения системы линейных уравнений, который стал помехой для решения задач большой размерности и большом количестве полуциклов расчета при моделировании циклического нагружения. Если время счета на тестовом компьютере одного полуцикла расчета составляет около пяти минут при количестве уравнений порядка двадцати тысяч, то расчет тысячи по- луциклов - несколько суток, что значительно затрудняет анализ конструкции. Поэтому вста- ла необходимость модернизации модуля решателя и реализации более эффективных методов хранения матрицы жесткости конструкции и решения СЛАУ. При модернизации модуля решателя на основе метода Холесского (см. Приложение) реализованы несколько алгоритмов решения разреженных СЛАУ [4]: профильный; минимальной степени; дерево факторов; параллельных сечений; вложенных сечений. В зависимости от рассчитываемой конструкции сетки конечных элементов имеют различную структуру. Соответственно матрицы получаемой СЛАУ имеют различные характе- ристики (полуширину ленты, структуру заполненности и общую заполненность, т.е. количе- ство ненулевых элементов), в предположении, что нумерация узлов сетки оптимизирована для получения диагональной СЛАУ с минимальной полушириной ленты. На рисунке 1 при- ведены примеры конструкций с различными характеристиками. Рисунок 1. Варианты конструкций: а - сечение диска; б - сектор замкового соединения; в - два диска, соединенные по ободу и ступице; г - два диска, соединенные по полотну; д - три диска, соединенные по ободам и ступицам Конструкция на рисунке 1,а обладает свойствами балочной структуры, и, соответ- ственно, СЛАУ имеет малую полуширину и большую заполненность. Конструкция на ри- сунке 1,б имеет квадратную структуру, и, соответственно, большую полуширину ленты СЛАУ и малую заполненность. Остальные конструкции обладают промежуточными харак- теристиками. Для тестирования алгоритмов хранения матриц СЛАУ и их решения были использова- ния несколько сеток конечных элементов с различной конфигурацией и, соответственно, различной полушириной и степенью заполненности матрицы. Основные характеристики се- ток представлены в таблице 1. Из таблицы видно, что при практически одинаковом количе- стве уравнений матрицы жесткости задачи наименьшей полушириной и наибольшей запол- ненностью ленты обладает балочная структура, и, наоборот, наибольшей полушириной и наименьшей заполненностью - квадратная структура. Остальные структуры имеют проме- жуточные характеристики. Характеристики сеток конечных элементов Таблица 1 В таблицах 2-5 и на рисунках 2-3 представлены результаты решения тестовых задач различными методами. В таблице 2 представлены объемы требуемой памяти для накладных расходов (НР), т.е. массивов индексации систем хранения матрицы жесткости, самой матри- цы жесткости и всей задачи в целом для каждого из методов в единицах ячеек памяти про- граммы (I4 - целое число длиной 32 байта, R8 - число с плавающей точкой двойной точно- сти размером 64 байта), в таблице 3 - в процентах относительно объема всей задачи при ре- шении ленточном методом. На рисунке 2 приведены диаграммы объема требуемой памяти, построенные по таблице 3. Из таблиц и диаграмм видно, что наилучшие результаты по использованию памяти дает метод параллельных сечений. При этом требуемый объем памяти уменьшился от 40% на ба- лочной структуре до 19% на квадратной структуре. В таблице 4 представлены времена счета отдельных этапов расчета и всей задачи в це- лом для каждого из методов в процентах относительно этапов при решении ленточным ме- тодом. В таблице 5 представлены времена счета отдельных этапов итерации относительно времени итерации при решении ленточным методом. На рисунке 3 приведены диаграммы времен счета, построенные по данным таблицы 5. В этап подготовки входит нахождение упорядочения матрицы в соответствии с системой хранения и распределение памяти с созда- нием массивов индексации, в этап формирования глобальной матрицы (ГМ) - вычисление матриц и правых частей элементов и вставка их в глобальную матрицу и правую часть, в этап решения СЛАУ - нахождение разложения матрицы и вычисление вектора неизвестных, в этап итерации - суммарное время одной итерации при решении нелинейной задачи, состо- ящей из формирования матрицы и решения СЛАУ, в этап задачи - суммарное время решения задачи, состоящей из подготовки и нескольких итераций. Требуемый объем памяти (в ячейках памяти) Таблица 2 Требуемый объем памяти (в % от ленточного метода) Таблица 3 Из таблиц видно, что хотя все методы проигрывают ленточному методу при подготовке решения и формировании матрицы жесткости конструкции, время решения самой системы уравнений существенно выше. Это дает большой выигрыш по времени на одной итерации (рисунке 3), а также при решении всей задачи, особенно при нескольких итерациях, не гово- ря уже о нескольких сотнях и тысячах циклов. Наибольший выигрыш времени на одной ите- рации дает метод вложенных сечений - уменьшение времени итерации на сложных структу- рах до 9-29%. а) б) в) г) д) Рисунок 2. Объем памяти для различных задач: а - балочная; б - квадратная; в - кольцевая; г - тавровая; д - замкнутая тавровая. Методы решения: 1 - ленточный; 2 - профильный; 3 - минимальной степени; 4 - дерево факторов; 5 - параллельных сечений; 6 - вложенных сечений Время счета отдельных этапов (в % от ленточного метода) Таблица 4 Время одной итерации (в % от ленточного метода) Таблица 5 а) б) в) г) д) Рисунок 3. Время счета одной итерации разными методами для различных задач (за 100% принято время одной итерации при решении ленточным методом): а - балочная; б - квадратная; в - кольцевая; г - тавровая; д - замкнутая тавровая. Методы решения: 1 - ленточный; 2 - профильный; 3 - минимальной степени; 4 - дерево факторов; 5 - параллельных сечений; 6 - вложенных сечений Реализованные методы показали свою высокую эффективность. В конструкциях, где заведомо хорошо работает ленточный метод, время одной итерации удалось уменьшить до 67%, а объем занимаемой оперативной памяти - до 42%. В остальных же случаях время уда- лось уменьшить до 9%, а использование памяти - до 19%. На рисунке 4 представлены рейтинги для всех методов с учетом всех типов структур. В среднем наибольшую экономию памяти дает метод параллельных сечений, наибольшую экономию времени счета одной итерации - метод вложенных сечений. На рисунке 4,в пред- ставлен обобщенный рейтинг с учетом использования памяти и времени счета одной итера- ции. Видно, что наибольшей эффективностью обладает метод вложенных сечений. Таким образом, реализация приведенных методов позволяет повысить эффективность расчета более чем в 4 раза. При этом можно либо достичь уменьшения времени счета при решении задач от 10 раз, либо увеличить объем решаемых задач до 5 раз при сохранении тех же затрат машинного времени. а) б) в) Рисунок 4. Оценка эффективности реализованных методов с учетом всех типов задач: а - по использованию памяти; б - по времени счета одной итерации; в - обобщенный рейтинг. Методы решения: 1 - ленточный; 2 - профильный; 3 - минимальной степени; 4 - дерево факторов; 5 - параллельных сечений; 6 - вложенных сечений Выводы В модуле решателя МКЭ «Cycle2D» реализованы алгоритмы решения разреженных систем линейных алгебраических уравнений, что, как показано на тестовых задачах, сократило время их решения от 1,5 до 10 раз, уменьшило объем требуемой памяти от 1,5 до 5 раз, по- высило среднюю эффективность решения от 2 до 4 раз по сравнению с базовым ленточным методом. Анализ показывает, что наиболее эффективным алгоритмом для любой структуры яв- ляется алгоритм вложенных сечений. Литература Зенкевич О. Метод конечных элементов в технике. М.: Мир, 1975. 544 с. Бате К.-Ю. Методы конечных элементов / Пер. с англ. В.П. Шидловского под ред. Л.И. Турчака. - М.: ФИЗМАТЛИТ, 2010. - 1024 с. Темис Ю.М., Азметов Х.Х. Математическое моделирование циклического деформирова- ния // Известия МГТУ «МАМИ» 2011, с. 195-202. Джордж А., Лю Дж. Численное решение больших разреженных систем уравнений. Пер. с англ. - М.: Мир, 1984. - 333 с. Метод Холесского Приложение. При анализе конструкции методом конечных элементов задача сводится к решению си- стемы линейных алгебраических уравнений (СЛАУ), занимающей более 90% времени счета. Основными свойствами матрицы данной системы является большая размерность и разре- женность. Предложен ряд алгоритмов решения такой СЛАУ, описанных в работе [4]. Поскольку матрицы получаемой СЛАУ содержат большое количество нулевых элемен- тов, первым шагом каждого алгоритма является переупорядочивание матрицы для того, что- бы собрать все ненулевые элементы вместе, исключив таким образом нулевые элементы из процесса счета. Используя разреженность, можно достигнуть больших сокращений в запро- сах к ресурсам памяти и вычислительной работе. Имеются различные схемы хранения раз- реженных матриц, отличающиеся способом использования нулей. В некоторых случаях до- пускается хранение части нулей в обмен на упрощение схемы хранения, в других - исклю- чаются все нули системы. Полученная пересортированная матрица решается методом Холесского. Общий алгоритм решения можно записать следующим образом: Упорядочение: найти «хорошее» упорядочение (перестановку Р) для данной матрицы А с учетом выбранного метода хранения матрицы. Распределение памяти: определить необходимую информацию о множителе Холесского L матрицы PAPT с тем, чтобы сформировать подходящую схему хранения. Разложение: разложить переупорядоченную матрицу PAPT в произведение LU. Решение треугольных систем: решить системы Ly = b и LTz = y. Получение результата: положить x =PTz. Все используемые методы основаны на численном алгоритме, известном как метод Хо- лесского, - симметричном варианте гауссова исключения для симметричных положительно определенных матриц. Каждый метод использует его при конечном расчете результата, а также некоторые из них пользуются им в процессе сортировки. Предположим, что система уравнений, которую нужно решить, есть: Ax b , (1) где: А - симметричная положительно определенная матрица коэффициентов (матрица жест- кости) размерности N×N; b - вектор длины N, называемый правой частью; х - век- тор-решение длины N, компоненты которого нужно вычислить. Применение к матрице А метода Холесского приводит к треугольному разложению: A LLT , (2) где: L - нижняя треугольная матрица с положительными диагональными элементами. Подставляя данное представление в систему уравнений, имеем: LLT x b . (3) Замена LTx = y показывает, что вектор x можно получить, решая треугольные системы: Ly b и LT x y . (4) Вычислительную схему алгоритма метода Холесского, в так называемой форме внеш- них произведений, можно описать в матричном виде: T d1 0 1 0 1 A H d1 1 T d1 1 1 d 1 0 L LT L A LT H I 0 1 d 1 1 N 1 1 0 H1 d1 1 0 I N 1 1 1 0 H 1 1 1 1 , 1 0 0 1 0 A 0 d T 1 0 0 1 0 0 0 d 0 0 1 0 1 0 0 T 2 0 d 2 2 2 2 L A LT (5) 1 1 0 H 2 2 0 2 H 2 2 d 0 2 I 0 0 H 2 T d 2 2 d N 2 2 0 0 2 2 I N 2 , …, AN 1 N N N L I LT , где: di - положительный скаляр; i - вектор длины N-i; а Hi - положительно определенная симметричная матрица порядка N-i для значений i от 1 до N. После N шагов алгоритма имеем: T T T T A L1L2 ...LN LN ...L2 L1 LL . (6) Можно показать, что L = L1 + L2 +…+ LN - (N -1) IN. Таким образом, i-й столбец L есть в точности 1-й столбец Li. В этой схеме вычисляются один за другим столбцы L. В то же время каждый шаг тре- бует модификации подматрицы Н, посредством внешнего произведения i iT/di. Результатом является матрица Hi, т.е. как раз та матрица, которую остается разложить. Порядок обраще- ния к элементам матрицы А в процессе разложения показан на рисунке П1. Рисунок П1. Порядок вычислений в форме внешних произведений Рисунок П2. Порядок вычислений в методе окаймления Другая формулировка процесса разложения - это метод окаймления. Пусть матрица А представлена в виде: M u A T , (7) u s M причем уже получено симметричное разложение LMLT ведущей главной подматрицы М порядка N-1. Тогда разложением А будет: L 0 LT w где: A M wT t 1 M , (8) 0 t 1 T 2 w LM u и t s w w . (9) Заметим, что разложение LLT подматрицы М также можно получить приемом окаймле- ния. Поэтому схема может быть описана следующим образом: Для i = 1, 2, .... N решить систему: Вычислить l i,i ai,i l1,1 li 1,1 i 2 l 2 i,k ... ... ... . 0 li 1,i 1 li,1 ... li,i 1 ai,1 ... ai,i 1 . (10) k 1 В этой схеме поочередно вычисляются строки L. К еще не разложенной части матрицы обращение не производится до тех пор, пока не будет вычисляться соответствующая часть L. Последовательность вычислений изображена на рисунке П2. При решении задач МКЭ со сложными граничными условиями используются множи- тели Лагранжа [3]. Получаемые матрицы заведомо не обладают положительной определен- ностью из-за нулевых элементов на диагонали, и использование описанного выше алгоритма невозможно. Но имеется возможность использовать усовершенствованный метод окаймле- ния для решения таких матриц. Если A - симметрическая матрица, у которой все главные миноры невырождены. Мож- но доказать, что для нее возможно разложение: A L diag 1 LT , (10) где: L - нижнетреугольная невырожденная матрица. Это делается по следующему алгоритму. Обозначим: ~ LT diag ~ 1 LT , (11) где: LT - верхнетреугольная матрица, отличающаяся от LT только знаками в некоторых строках. Их выбор осуществляется в зависимости от знака элемента на диагонали матрицы при разложении. ~T Если разложение выполнено для главного минора порядка n ( LnLn An ), то подобное разложение для главного минора порядка n+1 осуществляется как показано на рисунке П3. Рисунок П3. Усовершенствованный метод окаймления Чтобы выполнялось равенство, изображенное на рисунке П3, необходимо выполнение условий: LT n ~ b , ~T 2 c . (11) Из данных уравнений находится сначала вектор . После его определения из и с сле- дует найти , выбрав нужный знак. От этого выбора будет зависеть, станет ли n+1 строка Ln матрицы ~T «инвертированной». Если минор An 1 ~T Ln 1Ln 1 невырожден, то не может оказаться нулем.