Джордж, Лю - Численное решение больших разреженных систем уравнений (947498), страница 15
Текст из файла (страница 15)
СИММЕТРНЧНОЕ РАЗЛОЖЕНИЕ ПРООЧЕЛЬНЫМ ° ° ° е С ° ч ° ° ° ° ° * ° ° ИСПОЛЬЗУЕМЫЕ ПОДПРОГРАММЫ " ЕЬБЕУ. Сч ° ° ° ° ° ° С С Сечь ° ° ° ° ° ° ч ь ° С Сь ° ° ° С цикл по стрскыч 2 3,... НЕЯИВ С С С С С С С С С С С С С С С С С С С С Гл. Е. Ленточные и профильные методы ПОДПРОГРАММА ВЫЧИСЛЯЕТ РАЗЛОЖЕНИЕ Ее Е (ТРАНСПОНИРОВАННАЯ) ДЛЯ ПОЛОЖИТЕЛЬНО ОПРЕДЕЛЕННОИ МАТРИЦЫ А хРАИЯшдйся В ПРовильном формате. НспольБУПГСР СТАНДАРТНЫЙ МЕТОД ОКАЙМЛЕНИЯ.
ВХОДНЫЕ ПАРАМЕТРЫИЩИБ - ЧИСЛО УРАВНЕНИЙ. ХЕИУ - ИНДЕКСНЫЙ ВЕКТОР ОБОЛОЧКИ, ИЗМЕНЯЕМЫЕ ПАРАМЕТРЫЕИУ - ОБОЛОЧКА Ь - НА МЕСТО ОБОЛОЧКИ А . ВТАО - ДНАГОНАлы.- нА место ДИАГОНАпй А. 1РЬАО - индикктоР ошивки. ему ПРИСВАИВАйтся ч ЕСЛИ ПРИ РАЗЛОЖЕНИИ ОБНАРУЖЕН О ИЛИ КОРЕНЬ ИЗ ОТРИЦАТЕЛЬНОГО ЧИСЛА, ° ° ° чч ьч ° ч*ч ь* ° ° чь ° чье ° ч БОВЕОВТ1ИЕ ЕБРСТ ( ИЕОИБ, ХЕИУ, ЕИУ, 01АО, 1И.АО ) 000ВЕЕ РПЕС1БТОИ СОПИТ, ОРБ СОИИОИ /БРЕОРБ/ ОРБ ЕЕА(.
01АО(1), ЕЕЧУ(1), Б, ТЕИР 1ИТЕОЕЕ ХЕИУ(1), 1, 1ВАИО, 1Р1ВБТ, 1ИА6, 1ХЕИУ, 1, 1 ВТОР, ИЩИБ 1Р ( 01АО(1) .ЕЕ. О.ОЕО ) 60 ТО 400 01АО(\) Щкт(01А6(1)> 1Р ( ИЗБ .Щ. 1 ) ЕЕТОЕИ ВО 300 1 2, ИЩИБ 1ХЕ>ЧУ ХЕИУ(1) 1В>ОЧО ХЕЬ>У(1+1 > - 1ХРЕ(У ТЕИР 01АО(1) 1Р ( 1ВАИО .Щ. 0 ) 60 ТО 200 1Р1ВБТ 1 - 1ВАИО ВЫЧИСЛИТЬ СТРОКУ 1 ТРЕУГОЛЬНОГО МНОЖИТЕЛЯ САЫ. Е1.Б1 У ( 1 ВАЕА), ХЕИУ ( 1Р1ВБТ ), ИЧУ, 01А6(1РТВБТ>, еиу(1хшчу»' )ВТОР ХЕИУ(1ь1) - 1 00 100 1 1ХЕИУ, 1БТОР Б ЕЕПЕ(1) У 45, Подпрограммы Е5РСТ, Е) 5(.Ч и Е()51Ч 96 ТЕМР ТЕМР - 5 5 СОМТ1М)Е 1Р ( ТЕМР .1.Е.
О.ОЕО ) 60 ТО 400 01АО(1) 5ОИТ(ТЕМР) СООНТ 1ВАМ) ОР5 ОР5 + СООГ(т СОЙ11МЗЕ ИЕТПИН С С о 6 400 УСТАНОВИТЬ 1 ДЛЯ 1РЕАО( МАТРИЦА — НЕ ПОЛОЖИТЕ)ЗЬНО ОПРЕДЕЛЕННАЯ. 1РЕА6 1 ИЕ'ПЗИМ енп где М вЂ” ведущая главная подматрица А с разложением Холесского Тном, то множитель Т, матрицы А можно записать как Т где Тммг = и н 1=(5 — а('тр)'". Тем самым множитель Холесского можно вычислять строка за строкой, начиная с матрицы ам порядка 1 и переходя к матрицам М все большего порядка. Самое интересное в ЕЗРСТ вЂ” зто использование того обстоятельства, что векторы и — «короткие», так как мы работаем с профильной матрицей.
Обращаясь к рис. 4,5.2, предположим, что закончены первые 1 — 1 шагов разложения, так что разложена ведущая главная подматрица А порядка ( — 1. (Тем самым выполнены операторы, предшествующие циклу РО 300 1 = 2, ...; сам цикл был выполнен ( — 2 раза.) Чтобы вычислить ('-ю строку А, нужно решить систему уравнений г.лгп) и, Рис.
4.5.2. Иллюстрация и способу использоваиия разреженности в подпро. грамме ЕБГСТ; в вмчислеиии ш по и участвует лишь г,. 96 Гл. 4. Ленточная и профнльнае метода Однако нз рисунка (а также из лемм 2.2.! н 2.2.4) видно, что в вычислениях участвует лишь часть 1.м, именно та часть, что обозначена г.. Поэтому прн обрашеннн к ЕЬВ).Ч размер треугольной системы, т. е.
порядок и н Е, задает переменная 1ВАИ0; прн этом 1Р1гсЯТ есть индекс в Е первой строг1н 1., Упражнения 4Л.1. Пусть А имеет снмметрвчную структуру, но А чь А'1 в пусть гауссово исключеииц примененное к А, численно устойчиво без выбора главных элементов. Уравнения метода окзймлення для разложения А, аналогичные ура. вненням, использовавшимся в разделе 4.5.2 для подпрограммы ЕЗГСТ, будут М р юзг й хгм л Езгя н, уумы з н г ~а ыгй Здесь Ь вЂ” нижняя треугольная матрица с единицами на дяагонвлн в, разу.
меется, й чь ггг. з) Приняв в качестве основы ЕЕЗЕЧ, составить фортран-подпрограмму Е).!ЗьЧ, решающую ннживе треугольные системы с единичными диагональными злементами с использованием профильной схемы хранения. б) Используя ЕЗРСТ как основу, составить фортран-подпрограмму ЕГСТ, которая бы расхладывалз А в произведение Щ где ь н 0" хранятся по профильной схеме.
в) Какие подпрограммы понадобятся вам для решения системы Ал Ь, где А — матрица описанного выше типа. Указания: 1) В подпрограммах ЕЕЗЕЧ н ЕЗГСТ нужны лишь очень небольшие изменения. 2) Взшз реализация ЕРСТ должна использовать ЕЕ1ЯЕЧ в ЕЕЗЕЧ. 4.$.2. Пусть ь н Ь имеют структуру, изображенную иа рисунке, причем ь хранится массивами ХЕ)чЧ, ЕМЧ н 01АСг, как описано в рзед. 44,1. Сколько арифметических операций проделает подпрограмма Е).З).Ч прн решении системы Ех = Ьз Сколько операций выполнит Е))З1Ч при решении си. стемы х,гл Эз З 4.6. Дополнительные замечания То, что мы не проявляем энтузиазма в отношении ленточных упорядочений, отчасти связано с тем обстоятельством, что й 4о Дололнительные ламе«ание 97 в нашей книге рассматриваются методы, не использующие внешнюю память.
Если же приходится прибегать к внешней памяти, то ленточные упорядочения становятся привлекательными, поскольку очень легко составить подпрограммы разложения и решения, использующие внешнюю память и примерно р(р+ 1)/2 слов оперативной памяти (Еейрра 1970). В работе (%!1зоп е1 а!. 1974) описана ленточная схема, обходящаяся еще меныцим объемом оперативной памяти; программа авторов может работать даже тогда, когда памяти хватает лишь на хранение двух столбцов из ленты 7.. Другая ситуация, в которой важны ленточные упорядочения, — это применение так называемых ленточных методов с минимальной памятью (ЗЬеггпап 1975). Основная вычислительная схема та же, что и для методов с внешней памятью, однако здесь столбцы 7.
вычисляются, используются и затем . «отбрасываются», вместо того чтобы быть записанными во внешнюю память. Части 7., которые понадобятся в дальнейшем, перевычисляются. Было предложено несколько алгоритмов упорядочения для получения малого профиля, отличающихся от приведенных в этой главе. В работе (! ему 1971) описан алгоритм, иумерующий узлы по принципу наименьшего приращения в размере оболочки. Кинг (К!пп' 1970) предложил аналогичную схему с тем отличием, что круг кандидатов на помечивание ограничен такими узлами, у которых есть по крайней мере один нумерованный сосед; следовательно, здесь требуется выбор начального узла.
В последние годы были 'предложены некоторые алгоритмы, более близкие к описанному нами (О)ЬЬз е1 а!. 1976Ь, БВЬЬз 1976с). Некоторые авторы (Ме!озЬ !969, 1гопз 1970) описали «фронтальные» методы (методы «волнового фронта»), использующие вариации в ширине ленты и рассчитанные на привпечение внешней памяти. В этих схемах требуется приблизительно ео(ео+ 1)/2 слов оперативной памяти вместо р(р+ 1)/2 слов в случае ленточных схем, хотя программы, как правило, оказываются значительно более сложными. Идея фронтальных методов была высказана в контексте решения конечноэлементных систем уравнений, и еще одной новой чертой этих методов является то, что в них уравнения генерируются параллельно с процессом решения системы.
Было доказано (СЬап 1979), что при заданном начальном узле алгоритм 'КСМ можно реализовать так, чтобы он выполнялся за время 0(!Е)). Поскольку. каждое ребро графа должно быть исследовано хотя бы раз, этот метод, очевидно, является оптимальным. Набор подпрограмм, аналогичных подпрограммам Е(.$(.У, ЕПЯЛ и ЕБГСТ, содержится в работе ХЕ(зепз1а1 е974). Ю.
Универсальные разреженные методы $5.0. Введение В этой главе будут рассмотрены методы, которые в отличие от методов главы 4 стараются использовать все нулевые элементы треугольного множителя Е матрицы А Метод упорядочения, изучаемый в данной главе, называется алгоритмом минимальной степени (Козе 1972а). Это эвристический алгоритм для отыскания такого упорядочения матрицы А, при котором она в процессе разложения претерпевает лишь небольшое заполнение. Алгоритм широко используется в прикладных расчетах и имеет хорошую репутацию. Машинные реализации подпрограмм распределения памяти и численного решения суть адаптации соответствуюших подпрограмм Йельского пакета для разреженных матриц (Е1зепз1а1 1981).
ф 5.1. Симметричное разложение Пусть А — симметричная разреженная матрица. Структура ненулевых эЛементов А определяетея как Нопг(А)=((1, 1) ~ам ~ О и 1~1). Предположим, что посредством алгоритма Холесского матрица разложена в произведение ЕТ.т,.Матрица заполнения Р(А) для А есть матричная сумма Е + Ет. Когда из контекста ясно, о какой матрице речь, будем писать Р вместо Р(А). Соответствующая ей структура ненулевых элементов— )х)опг(Р)=((ю',!)1)пФО и 1Ф1), Напомним, что во всей книге мы придерживаемся соглашения о том, что точное числовое взаимное уничтожение не имеет места; поэтому прн заданной структуре Хопг(А) соответствуюшее множество Хопг(Р) полностью определено.