Другие журналы

научное издание МГТУ им. Н.Э. Баумана

НАУКА и ОБРАЗОВАНИЕ

Издатель ФГБОУ ВПО "МГТУ им. Н.Э. Баумана". Эл № ФС 77 - 48211.  ISSN 1994-0408

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

#10 октябрь 2005
И

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

 

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

 

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

 

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

                                                       (1)

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

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

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

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

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

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

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

 

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

 

 

1

2

3

4

5

6

7

8

9

1

a11

a12

a13

0

a15

a16

0

0

0

2

a21

a22

a23

0

a25

a26

0

0

0

3

a31

a32

a33

0

0

0

0

0

0

4

0

0

0

a44

a45

a46

a47

a48

a49

5

a51

a52

0

a54

a55

a56

0

0

a59

6

a61

a62

0

a64

a65

a66

0

0

a69

7

0

0

0

a74

0

0

a77

a78

a79

8

0

0

0

a84

0

0

a87

a88

a89

9

0

0

0

a94

a95

a96

a97

a98

a99

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

 

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

                       (2)

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

В уравнении (2) присутствуют суммы произведений двух типов. Первая сумма состоит из произведений элементов aij- матрицы 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).

 

 

 

 

1

a11

1

a15

2

a22

2

a16

3

a33

3

a12

4

a44

4

a13

5

a55

5

a25

6

a66

6

a26

7

a77

7

a23

8

a88

8

a47

9

a99

9

10

a48

 

a49

11

a46

12

a45

13

a59

14

a56

15

a69

16

a78

17

a79

18

a89

 

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

 

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

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

 

 

 

 

 

 

 

1

17

 

1

1

 

1

1

 

1

5

2

1

 

2

1

 

2

5

 

2

6

3

2

 

3

2

 

3

0

 

3

2

4

0

 

4

4

 

4

8

 

4

3

5

4

 

5

2

 

5

13

 

5

5

6

7

 

6

1

 

6

15

 

6

6

7

11

 

7

4

 

7

1

 

7

3

8

12

 

8

5

 

8

18

 

8

7

9

13

 

9

2

 

 

 

 

9

8

 

 

 

10

1

 

 

 

 

10

9

 

 

 

11

4

 

 

 

 

11

6

 

 

 

12

7

 

 

 

 

12

5

 

 

 

13

7

 

 

 

 

13

9

 

 

 

14

8

 

 

 

 

14

6

 

 

 

15

4

 

 

 

 

15

9

 

 

 

16

6

 

 

 

 

16

8

 

 

 

17

5

 

 

 

 

17

9

 

 

 

 

 

 

 

 

 

18

9

 

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

 

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

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

                                 (3)

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

                                                 (4)

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

              (5)

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

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

                                           (6)

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

        (7)

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

 

Номер КЭМ

Число узлов

Число элементов

1

756

125

16079

16071

16079

2

4961

1000

121709

121691

121709

3

35721

8000

946619

946619

946619

 

Сравнение данного алгоритма с алгоритмами, рассмотренными в работах [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

 

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

Поделиться:
 
ПОИСК
 
elibrary crossref ulrichsweb neicon rusycon
 
ЮБИЛЕИ
ФОТОРЕПОРТАЖИ
 
СОБЫТИЯ
 
НОВОСТНАЯ ЛЕНТА



Авторы
Пресс-релизы
Библиотека
Конференции
Выставки
О проекте
Rambler's Top100
Телефон: +7 (915) 336-07-65 (строго: среда; пятница c 11-00 до 17-00)
  RSS
© 2003-2019 «Наука и образование»
Перепечатка материалов журнала без согласования с редакцией запрещена
 Тел.: +7 (915) 336-07-65 (строго: среда; пятница c 11-00 до 17-00)