И. В. Станкевич, канд. техн. наук, Научно-исследовательский институт энергетического машиностроения МГТУ им. Н. Э. Баумана

Хранение и использование разреженных матриц в конечно-элементной технологии

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

В настоящее время в связи с бурным развитием вычислительных средств широкое распространение получили информационные технологии, имеющие разнообразную теоретическую и прикладную направленность. Среди них особое место занимают системы автоматизированного проектирования (САПР), неотъемлемую часть которых составляют подсистемы математического моделирования различных физических процессов. Одним из самых перспективных направлений развития математического моделирования как составной части комплексной автоматизации сквозного цикла: проектирование — конструирование — изготовление, является широкое использование численных методов, таких как метод конечных элементов (МКЭ), метод граничных элементов (МГЭ) и некоторых других. В процессе построения дискретных аналогов краевых задач указанными методами возникают большие системы линейных, а в общем случае нелинейных, алгебраических уравнений. Нелинейные системы уравнений решаются в два этапа: на первом этапе они линеаризуются, а затем полученная система линейных уравнений решается с помощью какого-либо метода линейной алгебры. Если сходимость не достигнута, то процесс повторяется. Таким образом, каждый раз решается система линейных алгебраических уравнений (СЛАУ):

                                                       (1)

где  — вектор узловых неизвестных;  — вектор правой части.

Матрицы A СЛАУ (1), как правило, симметричны с выраженной разреженной структурой. При соответствующей нумерации узлов конечно-элементных моделей (КЭМ) эти матрицы могут иметь ленточную структуру. В связи с этим при использовании конечно-элементной технологии возникает проблема разработки эффективных алгоритмов формирования, хранения и использования разреженных матриц. Вопросы организации, хранения и использования разреженных матриц рассмотрены во многих работах, например [1—6].

Память, используемая для хранения разреженных матриц, состоит из двух частей: основной памяти, в которой содержатся числовые значения элементов матриц, и дополнительной памяти, где хранятся указатели, индексы и другая информация, необходимая для формирования структуры матриц и обеспечивающая доступ к числовым значениям их элементов при выполнении процедур формирования и решении СЛАУ, т. е. так называемые списки связности. Способы хранения и использования данных, хранящихся в основной и дополнительной памяти, весьма разнообразны и определяются, главным образом, выбранным методом решения СЛАУ. Для структурного построения информации, хранящейся в дополнительной памяти — списков связности — широко используется теория графов [2, 4].

Особенно сложные алгоритмы работы с разреженными матрицами возникают при использовании прямых методов решения, например типа фронтальных [1—4]. В первую очередь это связано с тем, что при разложении разреженной матрицы A = LL^T она обычно претерпевает некоторое заполнение, т. е. нижний треугольный множитель L имеет ненулевые элементы в тех позициях, где в исходной матрице A стояли нули (не исключено, что при этом могут появиться и новые нулевые элементы). Кроме того, прямые методы требуют специальной нумерации узлов конечно-элементной сетки, которая обеспечивает минимальную ширину ленты. Подобных проблем при использовании итерационных методов не возникает.

В рамках данной работы применительно к итерационным методам решения СЛАУ (1) был разработан алгоритм хранения и использования разреженных матриц, отличающийся простотой реализации и высокой вычислительной эффективностью.

При реализации итерационного решения системы уравнений (1) весьма частой является ситуация, когда необходимо выполнить умножение матрицы системы A на какой-либо вектор, например, на вектор узловых неизвестных  или на вектор невязки , где . Построение произведений типа  или  является одним из узких мест эффективной реализации всего итерационного процесса решения системы уравнений (1), поскольку требует наибольших затрат процессорного времени. Ниже представлен алгоритм построения произведения  и непосредственно связанный с ним способ плотного хранения матрицы A.

Рис. 1. Конечно-элементная модель

Рассмотрим в качестве примера конечно-элементную модель, представленную на рис. 1. С каждым узлом сетки связано некоторое число узлов, которые могут быть разбиты на два множества: первое множество составляют узлы, номера которых меньше i — номера рассматриваемого узла, а второе множество — узлы, номера которых больше i. На рис. 2 показана структура матрицы СЛАУ, соответствующая узловым связям конечно-элементной модели, представленной на рис. 1. Матрица A хранится в виде двух векторов (рис. 3):  — члены главной диагонали, — ненулевые элементы, расположенные над главной диагональю. Для формирования и использования элементов вектора  создаются дополнительные четыре вектора, содержащие информацию о связях узлов конечно-элементной модели.

Рис. 2. Структура матрицы A системы линейных алгебраических уравнений

Рассмотрим i-е уравнение СЛАУ (1). В общем случае оно имеет вид

                       (2)

где a_ij — элементы матрицы A; f_i — элемент вектора правой части ; x_j — элемент вектора неизвестных ; N — число неизвестных.

В уравнении (2) присутствуют суммы произведений двух типов. Первая сумма состоит из произведений элементов a_ij - матрицы A и элементов дивектора неизвестных , у которых j < i, а для произведений второй суммы имеем j > i. Для построения произведений первой суммы используются векторы  и  пример структуры которых показан на рис. 4. В векторе  хранятся номера тех узлов, которые связаны с данным рассматриваемым узлом, и их номера меньше номера данного узла. Например, рассмотрим узел с номером 5 (см. рис. 1). Этот узел входит в элементы с номерами 2 и 3. В силу этого он имеет связь с узлами 4, 2 и 1, номера которых меньше 5. Номера этих узлов записываются в вектор  в последовательности, определенной номерами элементов и правилом обхода их узлов. На рис. 1 стрелками показано начало и направление обхода узлов элементов. Этим объясняется порядок появления номеров 4, 2 и 1 в векторе. В векторе  хранятся номера индексов, координирующих в векторе  расположение номера первого узла всего массива номеров узлов, связанных с рассматриваемым узлом. Например, как видно на рис. 1, массив узлов, связанных с узлом 5, начинается с индекса 6. Последний индекс массива определяется начальным индексом узла, следующего за рассматриваемым, за вычетом единицы. Например, для узла 5 последний индекс массива в векторе  равен 8 = 9 - 1, где номер 9 определяет начало массива узлов для узла 6 (см. рис. 1).

Рис. З. Структура векторного хранения матрицы A:  — вектор элементов главной диагонали;  — вектор наддиагональных элементов

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

В тех редких случаях, когда тот или иной рассматриваемый узел конечно-элементной модели не имеет связи с узлами, номера которых меньше или больше номера данного узла, тогда в соответствующие ячейки векторов  или  заносятся нули. Например, узел 4 (см. рис. 1) не имеет связи с узлами, номера которых меньше 4, это учитывается записью нуля в ячейку 4 вектора  (рис. 4). Узел 3 не имеет связи с узлами, номера которых больше 3 — этому обстоятельству соответствует запись нуля в ячейку 3 вектора . Исключение составляет запись в 1-ю ячейку вектора , куда заносится номер последней ячейки вектора . Это сделано для того, чтобы зафиксировать параметры оператора цикла "DO" при просмотре списка узловых связей последнего узла сетки конечно-элементной модели. Для всех других узлов эти параметры определяются тривиально. Для каждой конечно-элементной модели информационные векторы строятся один раз — сразу же после того, как закончено построение массива, описывающего конечные элементы сетки.

Рис. 4. Структура векторов-указателей

Структура вектора  (см. рис. 3) в значительной степени повторяет структуру вектора . Ненулевые наддиагональные элементы a_ij матрицы A хранятся в векторе  построчно, а порядок их следования в пределах строки определяется порядком следования индексов в векторе .

Теперь рассмотрим построение сумм произведений, входящих в уравнение (2). Рассмотрение начнем со второй суммы, а именно:

                                 (3)

При фиксированном индексе i определяются номера первой и последней ячеек вектора , в которых размещаются номера узлов, связанных с j-м узлом, и для которых j > i. Например, для узла 5 (см. рис. 1) имеем: номер 1-й ячейки 13, а номер последней ячейки 14 = 15 - 1. Номера первой и последней ячеек определяют в векторе  диапазон размещения ненулевых элементов матрицы A, участвующих в произведениях суммы (3). Кроме того, эти же номера ячеек определяют в векторе  диапазон размещения номеров соответствующих элементов x_i вектора . Таким образом, например, для узла 5 сумма (3) состоит из двух слагаемых:

                                                 (4)

Построение первой суммы из уравнения (2)

              (5)

алгоритмически выполняется сложнее. При фиксированном номере i определяются номера первой n и последней m (n < m) ячеек вектора , в котором хранятся номера j узлов, связанных с узлом i, и для которых j < i. Этим самым задаются компоненты x_i вектора . Теперь необходимо определить элементы a_ij = a_ji (i ≠ j) матрицы A, хранящиеся в векторе . Эта процедура выполняется следующим образом. Последовательно из ячеек с n по m в векторе  берутся номера узлов j (j < i). По этим номерам j из вектора  определяются номера первой к и последней l (l < k) ячеек вектора . Для каждого узла j в диапазо­не ячеек (k — t) вектора  находится номер i, поскольку, во-первых, узел i связан с каждым узлом j, а, во-вторых, i > j. Просматривая ячейки вектора  и сравнивая номера узлов, находящихся в этих ячейках, с номером i, можно определить номер ячейки, в которой для данного узла j и его диапазона номеров ячеек (k — l) находится номер узла i. По номеру ячейки узла i в векторе  можно найти в векторе  элемент a_ij = a_ji (i ≠ j, j < i) т- е- точно так же, как это делалось при построении второй суммы (2).

Рассмотрим пример. Построим сумму (5) для узла 5 конечно-элементной модели, показанной на рис. 1. С этим узлом, как уже отмечалось выше, связаны узлы 4, 2 и 1, Номера этих узлов хранятся в векторе  в ячейках с 4 по 6, причем их номера меньше 5. Таким образом, компоненты дивектора  известны: х₄, х₂ и x₁. Определим элементы матрицы A. Узлу 4 в векторе  соответствуют ячейки с 8 по 12. В этих ячейках хранятся номера узлов, связанных с номером 4, причем все эти номера больше 4. Один из этих номеров — 5. Этот номер (узла) находится в ячейке вектора  с номером 12. В ячейке с этим же номером 12 в векторе  хранится элемент a₄₅ = a₅₄ (j = 4 < i = 5) матрицы А. Аналогично отыскиваются элементы a₂₅ (j = 2 < i = 5) и «a₁₅ (j = 1 < i = 5). Таким образом, получаем выражение для :

                                           (6)

Окончательно с учетом (4) вся сумма S⁵ имеет следующий вид:

        (7)

Из представленного выше алгоритма видно, что вычислительные затраты на построение сумм произведений  и  разные. Вторая сумма () строится значительно быстрее. Это связано с тем, что при построении первой суммы (^) необходимо проводить дополнительный сравнительный анализ номеров узлов, размещенных в векторе . При рассмотрении процедур алгоритма предполагалось, что все необходимые данные располагаются в оперативной памяти компьютера. Длина одномерных массивов, в которых размещаются векторы ,  и , зависит от структуры рассматриваемой конечно-элементной модели. В таблице в качестве примера представлены данные по размерности этих векторов для трех типов конеч­но-элементных моделей, построенных в результате разбиения единичного куба на 20-узловые квадратичные элементы. Куб разбивался в следующих пропорциях: 5:5:5 (КЭМ № 1), 10:10:10 (КЭМ № 2) и 20:20:20 (КЭМ № 3).

Сравнение данного алгоритма с алгоритмами, рассмотренными в работах [1—6], показывает следующее. Для реализации предложенного алгоритма требуется больше оперативной памяти, поскольку для координирования ненулевых элементов используются не два массива (векторы номеров строк и столбцов), а четыре. Дополнительные векторы  и , размерности которых не превышают числа диагональных элементов, позволяют локализовать зоны поиска элементов-сомножителей при построении сумм (3) и (5). За счет этого достигается экономия временных затрат на решение СЛАУ (1).

* * *

В заключение отметим, что изложенный выше алгоритм формирования сумм в уравнении (2) может быть полностью применен и для построения других аналогичных матричных произведений, необходимых для реализации итерационных процессов.

Список литературы

1. Тьюарсон Р. Разреженные матрицы / Пер. с англ. М.: Мир, 1977. 191 с.

2. Брамеллер А., Аллан Р., Хэмэм Я. Слабозаполненные матрицы / Пер. с англ. М.: Энергия, 1979. 192 с.

3. Эстербю О., Златев 3. Прямые методы для разреженных матриц / Пер. с англ. М.: Мир, 1987. 120 с.

4. Писсанецки С. Технология разреженных матриц / Пер. с англ. М.: Мир, 1988. 410 с.

5. Сабоннадьер Ж.-К., Кулон Ж.-Л. Метод конечных элементов и САПР / Пер. с англ. М.: Мир, 1989. 192 с.

6. Станкевич И. В. Численные методы линейной алгебры: Учеб. пособие. М.: МГТУ им. Н. Э. Баумана, 1991. 44 с.

7. Самарский А. А., Николаев Е. С. Методы  решения сеточных уравнений. М.: Наука, 1978. 592 с.

8. Хсйгеман Л., Янг Д. Прикладные итерационные методы / Пер. с англ. М.: Мир, 1986. 448 с.

МАТЕМАТИЧЕСКОЕ МОДЕЛИРОВАНИЕ

ИНФОРМАЦИОННЫЕ ТЕХНОЛОГИИ, № 12, 1998

Ключевые слова: Метод конечных элементов, разреженные матрицы, СЛАУ, итерационные методы, векторы указателей, умножение матриц.