Форсайт Дж., Малькольм М., Моулер К. - Машинные методы математических вычислений (1040536), страница 13
Текст из файла (страница 13)
Приведенная ниже основная программа демонстрирует ис- пользование подпрограмм ОЕСОМР и 301.'уЕ. Заметьте, что заявленная размерность массива Л вЂ” переменная Х01М— имеет значение 10, в то время как Л' — действительный порядок матрицы — равно 3. В нашей практике обнаружилось, что люди, впервые пользующиеся этими программами, очень часто ошиба- ются, неправильно присваивая значение переменной Л)О(М. Значения для правых частей устанавливаются лишь после того, как выяснилось, что матрица достаточно хорошо обуслов- лена, чтобы стоило решать линейную систему.
При этом мы уста- З.З. ПОДПРОГРАММЫ ОЕСОМР И ЗОЬЧЕ ИЛЛЮСТРИРУЮЩАЯ ПРОГРАММА ДЛЯ ПОДПРОГРАММ РЕСОМР И ЗОЬУЕ С С [[ЕАЬ А ([О, 10), В (10), %0[[К (10), СОХР, СОХРР1 1ХТЕОЕЙ [РЧТ(10), 1, Л Х, Х0[М Х01М=10 Х=З А(1, !)=!О А(2, 1) = — 3 А(3, 1) =5 А(1, 2)=- — 7 А (2, 2) =-2 А[З, 2)=- — 1 А(1, 3) = — О А(2, З)=6 А (3, 3) =.5 00 1 1==1, Х %[[[ТЕ (б, 2) (А(!, 2), 2=1, Х) СОХТ[Х[[Е 2 РО[[МАТ (1Н, [ОР5.0) 40[[[ТЕ (б, 8)' САЬЬ ОЕСОМР(ХР!М, Х, А, СОХО, [РЧТ, %О[[К) [и[[!ТЕ (б, 3) СОХР 3 ГО[[МАТ(6Н СОХО=, Е15.5) 4УД[ТЕ (б, 8) СОХРР[ =.СОХР-[-1 1Г (СОХРР! .ЕО, СОКР) %К[ТЕ (6, 4) 4 РОНМАТ (40Н МАТК[Х 15 5!ХО[[САК ТО \ЧОКК[ХО РКЕС[5[ОХ) !Г (СОХРР! .Е().
СОХО) ЗТОР В (1) == 7 В (2) = 4 В (3) =-. 6 0051= — 1, Х 40[[[ТЕ(6, 2) В (!) 5 СОХТ[Х[7Е навливаем значения злементов матрицы и правой части именно посредством операторов присваивания лишь для того, чтобы не беспокоиться о формате данных. Для примера, использованного в 2 3.1, на выходе будет получена такая информация: 10, -7. О. -3. 2. 6. 5. — 1. 5. СОХО= 0.12608Е+02 7. 4.
б. 0.00000 -1.00000 1.00000 3.3. ПОДПРОГРАММЫ РЕООМР И 5ОЕИЕ 67 ОПРЕДЕЛИТЕЛЬ МАТРИЦЫ А МОЖЕТ БЫТЬ ПОЛУЧЕН НА ВЫХОДЕ ПО ФОРМУЛЕ ПЕТ (А)=!РУТ(К) 'А (1, 1) 'А(2, 2) в...вА (Н, Н). ЯЕАЬ ЕК, Т, АНОЯМ, УНОЯМ ХВОЯМ 1НТЕОЕй ХМ1, 1,), К, КР!, КВ, КАП, М 1РУТ (Н)=1 1Р (Н.ЕО.!) 60 ТО 80 НМ! =)( — 1 ВЫЧИСЛИТЬ 1-НОРМУ МАТРИЦЫ А') АНОЯМ=О.О ВО 10Д=1, Н Т=О.О !)О 5 1=1, Н т=т+АВ5<А <1, .))) СОНТ15 ЦЕ 1Р (Т .ОТ.
АНОЕМ) АНОЕМ=Т СОНТ1Н() Е ГАУССОВО ИСКЛЮЧЕНИЕ С ЧАСТИЧНЫМ ВЫБОРОМ ВЕ. ДУЩЕГО ЭЛЕМЕНТА ОО ЗБ К=1, ХМ! КР1=К+1 НАЙТИ ВЕДУЩИЙ ЭЛЕМЕНТ 5 10 М=.К 1)О 15 1= КР1, Н !Р <АВ5<А<1, К)) .От. АВ5(А(М, К))) М=1 сонт!М<)е !РУТ (К) =- М 1Р (м .не. к) 1Рут (к) = — !рут (к) Т= А (М, К) Л(М, К)=А(К, К) А(К, К)=Т ПРОПУСТИТЬ ЭТОТ ШАГ, ЕСЛИ ВЕДУЩИЙ ЭЛЕМЕНТ РАВЕН НУЛЮ 15 1Р (Т .ЕО. 0.0) 60 ТО 05 ВЫс1ИСЛИТЬ МНОЖИТЕЛИ С С С ОО 20 1= КР!, К А (1, К) =- — А (1, к))т СО!!Т(КСГ 20 ПЕРЕСТАВЛЯТЬ И ИСКЛЮЧАТЬ ПО СТОЛБЦАМ 3) Для векторной нормы, позирую используют авторы, соответствуюизую катри ГНУю нормУ, определеннув в 2 Ги2, принято помечать индексом !.— Лупан перев. з линейные системы уРАВнениЙ 70 с 1Г(М .Е().
!) СО ТО 50 'Т(М! =М вЂ” ! 00 20 К=(, (4М! КР! =К+1 М = !РОТ(К) Т =- В(М) В(М) = В(К) В(К) = Т 00 10 1=КР1, М В(1)=В(1)-1-А(1, К) 'Т !О СОХТ1)4ОЕ 20 СО(4Т1)(ОЕ С С ОБРАТНАЯ ПОЛСТА?10ВКА С ))О 40 КВ= 1, )чм! КМ! =)4 — КВ К=КМ!+ ! В(К).=в(К))А(К, К) Т = — В(К) ОО Зо 1=1, КМ! В (1) =-В(!)+А(1, К)'Т Зо СО)(Т))ЛЬЕ 40 СОМ Т1)4()Е во в(н=в(Н(л(1, Н нет(!Еы Е(4Р 3.4. Бопъшие разреженные системы Методы исключения Многие линейные алгебраические системы имеют столь боль, шое число и уравнений и неизвестных, что хранить в оперативной памяти полную квадратную матрицу из и' элементов невозможно.
Обычно такие системы получаются при дискретизации дифференциальных уравнений, обыкновенных илн с частными производными, а также при расчете различных конструкций или цепей. Зачастую матрицы таких задач настолько разрежены, что быстродействующей памяти хватает для хранения всех ненулевых элементов вместе с информацией об их расположении. Как же нужно решать соответствующую систему линейных уравнений? Если применение гауссова исключения возможно, то оно остается очень экономичным, точным и полезным алгоритмом.
Исключение возможно, если удается разместить в памяти ненулевые элементы треугольных матриц, строящихся в ходе исключения, вместе с информацией о том, где онн помещены. Пусть Е(7 — массив, нижний треугольник которого есть матрица множителей, а верхний треугольник — триангулированиая матри- З.ь БОльшие РАЭРеженные системы 7! ца". Тогда г' (/обычно значительно более заполнен, чем А, хотя и остается все еще разре>кенным. Говорят, чта позиции, соответствующие нулевым элементам А и ненулевым элементам Е(7, заполнились в гауссовом исключении.
В принципе возможно предсказать точную степень заполнения, вызываемого исключением, однако детали могут быть столь сложны, что программист отступится от такой возможности. Для некоторых матриц размер заполнения оценить легко, Примероат является ленточная матрица. Договоримся пока рассматривать ленточную матрицу как заполненную„т.
е. игнорировать все нули внутри ленты. Если гауссово исключение можно проводить без выбора ведущих элементов, что действительно возможно и безопасно для одного класса матриц, называемых положительно определенньсим (Форсайт, Моулер, (1967)), то заполнения нет вообще: т' У имеет ту же ширину ленты, что и А. Если же выбор ведущих элементов необходим (для невырожденной матрицы общего вида), то заполнение ограничено лентой полуторной ширины по сравнению с лентой матрицы А. Ленточную матрицу А удобно хранить прямоугольным массивом длины и и ширины 2>п+1; более широкий ленточный массив ЕУ также удобно хранить и обрабатывать.
Отсюда заключаем: линейные системы с ленточными матрицами легко решаются посредством исключения при условии, что ленточный массив Е(7 можно разместить в быстродействующей памяти. Существует несколько программ для обработки ленточных матриц; они являются сравнительно простыми модификациями таких программ исключения, как ОЕСОМР и ВОКАЛЕ. Если необходимо принимать в расчет нули внутри ленточной структуры матрицы для экономии памяти либо сокращения числа арифметических операций, то структура данных и ее обработка могут значительно усложниться. Такая детальная структура может возникать при применении конечио-разностных методов к краевым задачам для эллиптических уравнений в частных производных с переменными коэффициентами.
Решение подобных систем посредством гауссова исключения все еще является исследовательской задачей, и универсальные программы этого сорта только начинают создаваться. Иногда полная или даже ленточная матрица слишком велика, чтобы разместить ее в быстродействующей памяти. В этом случае необходимо часть матрицы хранить во внешней памяти — на дисках или магнитных лентах. Если задача столь велика, то хотя гауссово исключение и возможно, один только объем вычислений делает его очень дорогостоящим. Время выполнения ') То есть верхняя треугольная матрица, получаемая прямым ходом гауссова исключения.— Прим, перев. х линейные системы углвнений 72 арифметических операций обычно значительно больше, время, нужное для операций обмена частями матрицы с внеш памятью.
Поэтому в целях экономии важно так организов ход вычислений либо работу операционной системы, чт процессор не простаивал, ожидая завершения обмена. Это м но сделать различными способамн. Однако готовых унив сальных программ пока не существует, хотя большинство м ных вычислительных установок имели опыт решения таких бо ших систем. В последние годы огромные усилия системных программистов~ были сосредоточены на том, чтобы создать у пользователя иллюзию, будто он располагает очень большой (так называемой вирту-' альной) памятью для своих данных, хотя в действительности эти данные разбиты на страницы нли сегменты, которыми постоянно идет обмен с внешней памятью. Таковы, например, системы машин Впггоцяйз В5500 или некоторых машин фирмы 1ВМ серий 350 и 370.
Наличие виртуальной памяти освобождает программиста от забот, связанных с обменом данными. Однако, цена этой свободы может для некоторых стратегий разбиения на страницы оказаться очень высокой: если программа вынуждена ждать, пока механизм замещения отыскивает каждую очередную строку матрицы, как это происходит на машине Внггоцййз В5500 в пакетном режиме, то время выполнения для большой матрицы может вырасти в недопустимой степени.
Однако, виртуальная память обычно объединена с режимом мультипрограммирования, н во время поиска очередной страницы процессор обычно принимает к исполнению другую программу. Во многих операционных системах время прерывания ие записывают на счет пользователя даже в том случае, когда другой программы нет. Следовательно, «стоимость>-выполнения матричной программы остается приблизительно той же, бьи или нет обмен страницами. Однако в любом случае обмен удлиняет время ожидания завершения работы программы, Итерационные методы Имеется важный класс систем линейных уравнений, для которых элементы матрицы А задаются какой-нибудь простой формулой и, следовательно, могут вычисляться по мере необходимости.
Например, при моделировании уравнения Лапласа конечно-разностными методами элементы А обычно равны одному из чисел О, ! или — 4; какое именно значение принимается, легко определить, исходя из геометрии сетки. В таком случае элементы можно вовсе не хранить, а генерировать их, когда они понадобятся. Кроме того, часто порядки и столь велики, что было бы невозможно хранить заполненный массив т'.17. э.к Большие РАЗРеженные системы 73 Желательно решать такие линейные слстемы Ах==Ь методами, которые вообще не меняют матрицу Л и требуют хранения лишь нескольких векторов длины п.
(Заметим, что вектор Ь обычно должен храниться, так же, как и вектор х.) Методы, отвечающие этим требованиям, существуют и называются итерационными. В них начинают с какого-нибудь пробного приближения хкв и выполняют некоторый процесс, исполь. зующий А, Ь и хьи и приводящий к новому вектору х"'. Затем процесс повторяют. На Ь-м шаге итерационного процесса по А, Ь и х'»-" получают х'»'. При соответствующих предположениях векторы хив сходятся к пределу прп Ь в . Имеется множество подобных итерационных процессов.
Наиболее успешные нз них обязаны своим успехом тесной связи с решаемой задачей. Поэтому редко можно встретить библиотечные подпрограммы итерационных методов. Хотя итерационный процесс математически может быть несложен, структура матрицы А, скорей всего, сложна и специальным образом связана с данной конкретной задачей.
Один простой итерационный метод обсуждается в 5 24 книги Форсайт, Моулер ((967); этот метод известен под различными названиями: процесс Либ»4ана, метод Гаусса — Зейделя или метод пюследовательные лакеи(ений. Здесь основной итерационный шаг состоит в том, что Ье уравнение решается отпосительяо гхй компоненты х; нового вектора х, причем для всех прочих компонент х берутся прежде вычисленные значения. Можно доказать, что метод сходится для некоторых типов матриц, например, для любой симметричной положительно определенной матрицы А. Однако сходимость обычяо медленная, У многих итерационных процессов численного анализа сходимость столь медленная, что очень важная задача — найти способ ускорения сходимости, например, векторов х'»' к решению. И действительно, алгоритмы, ускоряющие сходил4ость последовательностей, составляют важную часть вычислительной математики.