Джордж, Лю - Численное решение больших разреженных систем уравнений (947498), страница 34
Текст из файла (страница 34)
1, 1, 1 1$ТОР, 1$ТВТ, 1.АБТ, МВЬКБ, МС01., МВОИ, ВОИ, 1 ВОИ1, ВОИ2 С» ° ° ° » ' ° ° ° ° ° » ° »» ° ° » ° » ° ° » ° ° ° ° ° ° ° с ° ° ° ° ° ° ° ° ° тввьч ... Решение треУГОльных систем неЯВнОЙ С С С С С С С С ЗОВ ПРЯМОЙ ХОД ......, .. 00 БОО 1 1, МВЬКБ ВОИ1 ХВ1.К(1) ВОИа - ХВЬК<1+1) - 1 ЬАБТ ХМОМ2(ВОИ2+!) 1Р ( 1 .Щ.
1 .Ой. 1АБТ .Щ. ХМОМ2(ВОИ1) ) 60 ТО ЗОО ВЫЧЕСТЬ НЗ ЯНБ ПРОИЗВЕДЕНИЕ ВНЕДИАГОНАЛЬНОГО БЛОКА НА СООТВЕТСТВУЮЩУЮ ЧАСТЬ ЭТОГО ВЕКТОРА. 00 200 ВОИ ВОИ1. ВОИ2 1БТВТ ХМОМ2(йОИ) 1Р ( 1$ТВТ .Щ. 1.АБТ ) 00 ТО 300 1$тОР х)а)я(ВОЯ+1) - 1 1Р ( 1$ТОР .ЬТ. 1$ТВТ ) 60 ТО 200 $0. ОЕО СООМТ 1$ТОР - 1$ТВТ + 1 ОРБ ОРБ + СОВЕТ 00 100 1 )БТВТ, )БТОР СО!. М2$03$(1) Б Б + ВНБ(СОЬ) МОНЕ(1) СОМТ1 МОЕ ВНЕ(ВОИ) ВВБ(йОИ) 3 6 6.Б. Подпрограмм« ТЗРСТ и ТЗБ1У 229 сомтгмпе МКО« 80«2 - КО«1 + 1 САЫ.
ЕЬЗЬУ ( !ЧКО«, ХЕНЧ(КО«1), ЕМУ. 91АО(ЯО«1), КН5(ВО«!) ) САЬЬ ЕЦЗЬУ ( МЯО«, ХНЧУ(КО«!>, ЕМЧ, 91АО(ЯО«!>, КНБ(ЯО«!» гоо зоо 1 чоо с с с сомттмие ОБРАТНАЯ ПОДСТАНОВКА... 1Р ( МВГКБ .Щ, ! > КЕПЛКМ ЬАБТ ХВЬК(МВЬКБ) - 1 по Зоо 1 - 1, ЬАБТ ТЕМР(1) О.ОЕО СОМТ1МОЕ 1 МВЬКБ "Ош хвьк(1) СОЬ2 ХВЬК(1+!) — 1 1Р ( 1 .Щ. ! ) ЯЕП>КМ (.АБТ ХМОМЕ(С01.2+1) 1Р ( ЬАБТ .Ео. %ЧОХЕ(СОЬ!) ) 60 ТО 900 УМНО>НИТЬ ВНЕДИАГ.
БЛОК НА СООТВ. ЧАСТЬ РЕШЕНИЯ. ПРИБАВИТЬ РЕЗУЛЬТАТ К ТЕМР. 00 800 СОЬ С01.1, СОЬ2 5 ЯНБ(С01.) 1Р ( 5 .Щ. О.ОЕО ) 60 ТО 800 ЛБТЯТ ХМОЯ(СОЬ> 1Р ( Л5ТЯТ .Щ. 1.АБТ ) 60 ТО 900 Л5ТОР ХМОМХ(СОЬ+!) - 1 1Р ( Л5ТОР .!.Т.'ЛЗТЯТ ) 60 ТО 800 СОВЧГ ЛБТОР - ЛБТЯТ + 1 ОРЗ ОР5 + СОФТ ПО 700 Л ЛБ'ГЯТ, Л5ТОР КО« М25085(Л> тЕМР(ЯО«> - ТЕМЯ(КО«> + БтМОМЕ(Л) СОМТ1МИЕ СОМТ 1 М(Е 1 1 - ! СО!.1 ХВЬК(1) соьг х81л!(1+1) - 1 МСОЬ СОЬ2 - С01.1 + 1 СМЛ. ЕЬБЬУ ( МСОЬ, ХЕНЧ(СОЬ!), ЕМУ, ВТАО (сош >, тевчР (сох( > > САЬЬ ЕПБЬУ ( МСОЬ, ХБЧУ(СО(„1), ЕМУ, 91АО(СОЬ1) ° ТЕМР (СО!. 1 ) ) 00 1000 Л С01,1, СОЬ2 КНБ(Л) ЯНБ(Л) - ТЕМР(Л) СОМТ1МОЕ 60 ТО 800 ЕМО 800 С с С С 100 800 900 !ооо ') Это не относится к ! = !. — Лрам.
перев. стей !). При обратной подстановке подготовка к вызову Е(017 н Е(ло1'тг предусматривает накапливание в вещественном рабочем векторе ТЕМР произведений внедиагональных блоков и частей решения. По окончании работы подпрограммы в массиве ((Но будет находиться вектор решения, 230 Гл. б. Методы фактор-деревьев Упражнение азк<.
Пусть Ь и «' — матрицы, описанные в упр. 4.2.5, причем «' имеет только три ненулевых элемента в каждом столбце. Сравните число операций при вычислении проиэведения )ггХ-гЬ-'Р как )гг(Ь-г(Ь-')г)) и каи (Ргй-г)(Е-~)г). Считайте и и р большими числами, так что младшие члены можно игнорировать 2 6.6.
Дополнительные замечания Идея «отбрасывания» внедиагональных блоков множителя Ь матрицы А, обсуждавшаяся в данной главе, может применяться рекурсивно (беогде 1978с). Чтобы пояснить эту стратегию, предположим, что А — матрица блочного порядка р, а векторы х и Ь разбиты на части соответственно А. Через А<в< обозначим ведущую главную подматрицу блочного порядка Й, через х<,) и Ь<м — соответствующие части х и Ь. Наконец, определим подматрицы в А так же, как в (6А.1); аналогичным образом разобьем Ь, Это разбиение показано на рисунке для р = 5. В этих обозначениях систему Ах =Ь можно записать как (.") .";)( )=(") а разложение матрицы А — в виде <43 н! а т где Ам =Аэв — рвА<4) <'в.
Формально систему Ах = Ь можно решать следующим образом: а) Разложение: вычислить матрицу Авв и разложить ее г т в произведение ЬввЬм (Заметим, что поправку )ГьА<4<й'в можно вычислять по столбцам, затирая каждый столбец после того, как он использован,) э" <> ><оаоэнитееьные эамечании За! б) Реп<ение: б.1) Решить систему А<4>у<4> = Ьмь т 6.2) Решить систему Аммэ = Ьэ — !<ау<э>. б.З) Решить систему Ам>я<4> = Уэхэ. 6.4) Вычислить х<4> = у<э> — х<п. Отметим, что мы всего лишь использовали идеи 5 6.1 и упр. 6.1.2, позволяющие избежать хранения %<э, требуется толь- ко Уэ.
Главное здесь заключается в следующем: все, что нужно для решения системы блочного порядка 6 без хранения В'а,— это возможность решать системы блочного порядка 4. Очевид- но, что той же стратегией можно воспользоваться для решения систем блочного порядка 4 без хранения %'4 и т. д. Таким образом, мы получили метод, который, по видимости, решает систему блочного порядка р, требуя хранения лишь диагональных блоков Е и внедиагональных блоков У, исходной матрицы. Заметим, однако, что каждый уровень рекурсии тре- бует рабочего вектора х<,> (на шаге б.
3); поэтому наступит момент, когда более мелкое разбиение уже не даст уменьше- ния в запросах к памяти, С этой процедурой связано много интересных и неисследованных вопросов. Изучение степени при- менимости идей разбиения и отбрасывания ядляется, по-види- мому, плодотворной областью исследований. Блочные методы успешно применялись и в варианте с ис- пользованием внешней памяти (Ъоп Рцс(>з 1972). Значение р выбирается таким образом, чтобы наличный объем оператив- ной памяти выражался каким-либо подходящим кратным числа (Ж/р)э. Так как А разрежена, то некоторые из блоков будут нулевыми. В оперативной памяти хранится массив указателей; каждая его компонента либо указывает текущий адрес соот- ветствующего блока, если в нем есть ненулевые элементы, либо равна нулю.
Если (рХ р)-матрица указателей слишком велика для того, чтобы храниться в оперативной памяти, то эту матрицу также можно разбить на клетки и рекурсивно ис- пользовать ту же идею. Эта схема управления памятью, оче- видно, связана с определенными издержками. Практика пока- зывает тем не менее, что она является жизнеспособной альтер- нативой другим схемам, использующим внешнюю память, на- пример ленточным или фронтальным методам. Одно из до- стоинств этой схемы в том, что для реальных матричных опе- раций используются простые структуры данных, а именно только квадратные и прямоугольные массивы.
В работе (ЯЫег 1976) рассмотрено применение древовид- ных разбиений к задаче явного обращения матрицы и предло- жен алгоритм вычисления древовидного разбиения графа, от- личающийся от нашего, 7. Методы параллельных сечений для хонечноалементных задач $7.0. Введение В этой главе мы рассмотрим стратегию упорядочения, предназначенную главным образом для задач, возникающих в приложениях метода конечных элементов. Эта стратегия сходна с методом главы 6 в том отношении, что отыскивается древовидное разбиение и используются вычислительные идеи неявного решения и асимметричного разложения. Основное преимущество описываемого в этой главе алгоритма параллельных сечений — гораздо меньшие, как правило, запросы к памяти, чем в ленточных или древовидных схемах, обсуждавшихся в предыдущих главах. Для кпнечноэлементных задач, за исключением систем сверхвысокого порядка, методы данной главы с точки зрения требуемой памяти часто оказываются наилучшими среди всех методов, рассматриваемых в этой книге.
Время решения также будет очень мало, хотя время разложения может быть больше, чем у некоторых других методов. Так как упорядочения, изучаемые в этой главе, являются древовидными, то применимы методы хранения и вычисления из главы 6; по этой причине названные вопросы не будут здесь рассматриваться. Однако схема параллельных сечений требует несколько более сложной процедуры распределения памяти, чем та, что описана в разделе 6А,З. Эта более общая процедура распределения будет предметом 5 7.3. ф 7.1.
Пример — задача на (т Х 1)-сетке 7.!. 1, Упорядочение посредством параллельных сечений В этом разделе мы обсудим задачу на (т Х!)-сетке, которая мотивирует алгоритм из ф 7.2. Рассмотрим (тХ 1)-сетку, показанную на рис. 7.1,1; она имеет У = т1 узлов и т < 1, Ассоциированная конечноэлементная матричная задача Ах=й, которую мы будем изучать, имеет следующее свойство: при некотором упорядочении уравнений (узлов) от 1 до А! неравенство ач чь 6 означает, что узел 1 и узел ! принадлежат одиомч и тому же малому квадрату. К 7П.
Пример — задача на (лс Х1)-сетке 233 Пусть теперь о — целое число, такое, что 1 ( о (( 1. Выберем о вертикальных ливий (называемых в дальнейшем разделителями), которые рассекут сетку на о+ 1 независимых блоков приблизительно одинакового рззмера.
Это показано на Рис. 7.!Л. !тн Х 1)-сетка дли пз = 6 н 1 11. рис. 7.!.2 для о= 4. Пронумеруем вначале о+ 1 независимых блоков (построчио), а затем разделители, как указано стрелками на рис. 7.1.2. Матричная структура, которую такое упо- Š— а фи чт+! Рис. 7.1.2. Упорядочение (нсХ1)-сетки посредством параллельных сечений указано стрелкзмн Здесь а = 4, что приводит к разоиенню с 2а+ 1 = 9 членами, пронумерованными взятыми в кружок числами.
рядочение индуцирует на треугольном множителе Ь, показана на рис. 7.!.3, где заштрихованы внедиагональные блоки, подвергшиеся заполнению. Положим 1 — а Ь--; — р. Смысл этой величины — расстояние между секущими. 234 Гл. 7.
Методы лараллельныл сечений Рассматривая А и Ь как блочные матрицы, разбитые на дв подматриц (а = 2о+ 1), отметим прежде всего размерности различных клеток: 1)Ааа — тбХтб, если !~ай<:о+1. 2) Ааг — т Х тб, если й > о + 1 и ! ц~: о + 1. (7.1.1) 3)Ааг — тХт, если />о+1 и й>о+1. Разумеется, реально о должно быть целым числом, и если б не является целым, то ведущие о + 1 диагональных блоков не Рис. 7.!.3. Структура матрицы У., индуцированнан посредством параллельных сечений сетки рис 712 Заштрихованные области указывают позиции внедиагональных блоков, подвергшиеся наполнению будут в точности одинаковой величины.