Дж. Деммель - Вычислительная линейная алгебра, страница 22

Разреженные матрицы общего вида Под разреженной матрицей понимают матрицу, содержащую большое число нулевых элементов. С практической точки зрения, это число должно быть настолько большим, чтобы имело смысл применение алгоритма, избегающего хранения нулевых элементов и оперирования с ними. В главе б обсуждаются методы решения разреженных линейных систем, отличающиеся от гауссова исключения и его вариантов. Существует огромное количество методов для разреженных систем и выбор наилучшего из них часто требует значительной информации о матрице системы [24]. В этом разделе мы лишь кратко обрисуем проблематику разреженного гауссова исключения и дадим ссылки на соответствующую литературу и имеющееся программное обеспечение.

94 Глава 2. Решение линейных уравнений В качестве очень простого примера рассмотрим следующую матрицу: 1 .1 1 .1 1 .1 1 .1 .1 .1 .1 .1 1 1 1 1 .1 1 .1 1 .1 1 .1 .96 1 1 .1 .1 .1 .1 1 1 .1 .1 .1 .1 .1 1 .1 1 .1 1 .1 1 А' = = ь'Г 1 .1 1 .1 —.01 .1 —.01 .1 —.01 1 .1 .1 .1 .1 .99 —.01 —.01 —.01 .99 —.01 —.01 .99 —.01 .99 1 —.01 1 —.01 —.01 1 Как видим, теперь 1ь и У' совершенно заполнены и их хранение требует пз слов. В действительности, все нулевые элементы в А' были заполнены уже по- сле первого шага алгоритма, поэтому приходится выполнять такую же работу, как и для плотного гауссова исключения, т.

е. 2~пз операций. Она упорядочена таким образом, что при применении к ней гауссова исключения с частичным выбором главного элемента никаких перестановок строк не происходит, и называется (учитывая характер расположения ее ненулевых элементов) стрслоепдььой хьатрицей. Отметим, что ни один из нулевых элементов в А не был заполнен при проведении гауссова исключения, так что 1 и бь совместно могут быть сохранены на месте, прежде занимавшемся ненулевыми элементами исходной матрицы. Кроме того, если подсчитывать только существенные арифметические операции (исключая умножение на нуль и сложение с нулем), то их будет лишь 12 (4 деления при вычислении последней строки в 1 и 8 умножений и сложений при пересчете элемента (5, 5)) вместо $из 83 операций.

Более общо, для стреловидной матрицы А порядка п хранение требует лишь Зп — 2 машинных слова вместо п, а применение гауссова исключения 2 обходится в Зп — 3 операции с плавающей точкой вместо 2~из. Когда и велико, эти значения для памяти и числа операций становятся очень малыми по сравнению со случаем плотной матрицы. Предположим теперь, что вместо А задана матрица А', получающаяся из А записью строк и столбцов в обратном порядке. Это равносильно обращению порядка уравнений и неизвестных в линейной системе Ах = Ь.

При применении к А' метода СЕРР снова не происходит перестановок строк и с точностью до двух десятичных разрядов будет найдено разложение 2.7. Специальные линейные системы 95 Данный пример свидетельствует, что порядок строк и столбцов крайне важен для зкономии памяти и работы. Выбор оптимальных перестановок для строк и столбцов, минимизирующих память или работу, является чрезвычайно трудной задачей даже в том случае, если не нужно заботиться о выборе главных элементов для обеспечения численной устойчивости [как обстоит дело в алгоритме Холесского). На самом деле, эта задача ХР-полна [111], что означает: все известные алгоритмы для отыскания оптимальной перестановки требуют времени, эксионенциально растущего вместе с и; таким образом, для больших и эти алгоритмы намного менее эффективны, чем даже плотное гауссово исключение.

Поэтому нам придется обратиться к эвристическим подходам, среди которых есть несколько удачных. Некоторые из них будут проиллюстрированы ниже. Помимо трудности выбора хороших перестановок для строк и столбцов, имеются и другие причины, почему разреженные реализации гауссова исключения или алгоритма Холесского намного более сложны, чем соответствующие методы для плотных матриц. Во-первых, нужно организовать структуру данных для хранения только ненулевых элементов матрицы А; существует несколько таких общеупотребительных структур [93[. Далее, нужна структура данных для хранения новых ненулевых элементов, возникающих в 7 и П в ходе исключения.

Это означает, что либо структура данных должна динамически изменяться в алгоритме, либо мы должны найти дешевый способ ее предварительного вычисления без реального проведения исключения. Наконец, структура данных должна быть использована для минимизации числа операций с плавающей точкой, а также выполнения целочисленных и логических операций в количестве, самое большее пропорционапьном числу флопов. Иными шювами, мы не можем позволить себе тратить 0[из) целочисленных и логических операций для отыскания тех немногих операций с плавающей точкой, которые мы согласны выполнять. Более полное описание этих вопросов выходит за рамки данной книги (см. [114, 93)), однако мы дадим справку об имеющемся программном обеспечении.

Пример 2.9. Проиллюстрируем разреженный алгоритм Холесского на реалистическом примере, возникающем при моделировании задачи о смещении механической структуры под воздействием внешних сил. На рис. 2.9 показана простая сетка, нанесенная на двумерное сечение механической структуры с двумя внутренними полостями. Математически задача состоит в том, чтобы определить смещения всех узлов сетки (внутренних по отношению к структуре) под влиянием некоторых сил, приложенных к границе структуры. Узлы сетки пронумерованы числами от 1 до и = 483; в практических приложениях значения и были бы гораздо ббльшими. Уравнения, связывающие смещения и силы, приводят к системе линейных уравнений Ат = Ь, где каждому из 483 узлов соответствуют одна строка и один столбец, причем ам ~ 0 в том и только том случае, если узел 1 связан отрезком прямой с узлом ~. Отсюда можно вывести, что А — симметричная матрица; она оказывается также положительно определенной, поэтому для решения системы Ат = Ь можно использовать алгоритм Холесского.

Заметим, что А имеет лишь пг = 3971 ненулевых элементов из возможного их числа 483т = 233289; таким образом, А заполнена лишь на 3971/233289 = 1.7%. (Сходные задачи моделирования механических структур рассматриваются в примерах 4.1 и 5.1, где приведен подробный вывод матриц А.) 96 Глава 2. Решение линейных уравнений Механическая структура и сетка 0.$ -2.5 -3 1 -3 -2 -1 0 1 2 3 4 3 Рис.

2.9. Сетка для механической структуры. В верхней части рис. 2.10 показана та же сетка, что и на рис. 2.9, а в нижней части этого рисунка — расположение ненулевых элементов в матрице А. При этом 483 узла сетки пронумерованы в некотором «естественном» порядке: узлы логически прямоугольных подструктур нумеруются построчно и одна подструктура нумеруется вслед за другой. Ребра каждой подструктуры окрашены одним н тем же цветом; этим цветам соответствует окраска ненулевых элементов матрицы.

Каждая подструктура имеет метку «(т: «)», указывающую, что ей соответствуют в А строки и столбцы с номерами от 1 до ~. При этом подматрица А(1: у, 1: Я является ленточной матрицей с узкой лентой. (В примере 2.8 и разд. 6.3 описаны другие ситуации, где сетка порождает ленточную матрицу.) Ребра, связывающие различные подструктуры, окрашены красным цветом и соответствуют красным элементам матрицы А, наиболее удаленным от ее главной диагонали. На верхней паре картинок рис. 2.11 снова показана структура разреженности матрицы А при естественном порядке узлов, а также структура разреженности ее множителя Холесского Ь. Ненулевые элементы в Ь, отвечающие ненулевым элементам в А, окрашены черным цветом; новые ненулевые элементы в Ь, составляющие заполнение, окрашены красным цветом.

Всего в Ь 11533 ненулевых элемента, что более чем в пять раз превышает число ненулевых элементов в нижнем треугольнике матрицы А. На вычисление Е алгоритмом Холесского затрачивается лишь 296923 флопа, что составляет всего 0.8% от -'пз = 3.76. 101 флопов, которых потребовал бы алгоритм для плотных матриц. ~ Цветной вариант рнс. 2ЛО. см. нв обложке книги. — Прим. иерее.

2.7. Специальные линейные системы Узлы сетки пронумерованы в естественном порядке г Н5 0.5 -2.5 '-5 -4 -З 2 3 4 5 -г -4 о 50 100 150 200 250 300 350 400 450 100 200 300 400 пг 255 Рис. 2.10. Ребра сетки в верхней части рисунка окрашены и пронумерованы;они находятся в соответствии с ненулевыми элементами разреженной матрицы А в нижней части рисунка. Например, первые 49 узлов сетки (самые левые узлы, окрашенные зеленым цветом) соответствуют строкам и столбцам матрицы А с номерами от 1 до 49. Глава 2. Рептенне линейных уравнений 100 100 200 300 400 100 200 400 200 200 300 Я при естественном упорядочении О 0 100 200 300 400 пх Э971 А после применения обратного алгоритма Катхилла-Макки 0 0 100 200 ЭОО 400 пв Э971 А после применения алгоритма минимальной степени 0 0 100 ЭчО 300 400 пх 397! Рис.

2.11. Разреженность и число матрицы А. Множитель Холесского, число флопов =296923 0 100 200 300 400 пх 11533, МножительХолесского,числофлопов =181525 0 100 200 ЭОО 400 пх 907Э, Множитель Холесского, число флопов =198236 0 100 200 300 400 пх 8440, Жирно — новые ненулевые элементы флопов при различных упорядочениях 99 2.7, Специальные линейные системы Конечно-элементная сетка возле крыла (НАЯА) 1 о.т а1 о ол оа оа ол ое ае оз ол оо 4253 узла Структура разреженности для крыла ((чАЯА) 0 0 600 1000 1600 6000 6600 ШО ИОО 4000 лг = 2883! Рис.

2.12. Сетка возле крыла (ХАБА). Число ненулевых элементов матрицы Е и количество флопов при ее вычислении можно существенно изменить, переупорядочнвая строки и столбцы в А. На средней паре картинок рис. 2.11 показаны результаты, полученные с помощью одной популярной схемы переупорядочения, называемой обратным алгоритмом Катлхилл — Макки (114, 93]; ее назначение — попытаться превратить А в ленточную матрицу с узкой лентой. Как видим, в данном случае это Глава 2. Решение линейных уравнений вполне удалось: заполнение в Л уменьшено на 21% (с 11533 элементов до 9073), а число флопов — почти на 39% (с 296923 до 181525). Другая популярная схема упорядочения называется алгоритмом минимальной стпепеии [114, 93]; ее принцип — минимизировать заполнение на каждом шаге алгоритма Холесского.