Die Grundlehren der mathematischen Wissenschaften in Einzeldarstellungen Band 186 Handbook for Automatic Computation Linear Algebra J.H.Wilkinson, C.Reinsch Heidelberg New York Springer-Verlag Berlin Уилкинсон, Райнш Справочник алгоритмов на языке АЛГОЛ Линейная алгебра Перевод с английского С. П. ЗАБРОДИНА, В. Г. ПОТЕМКИНА, П. И. РУДАКОВА Под редакцией д-ра техн. наук профессора Ю. И. ТОПЧЕЕВА Москва «МАШИНОСТРОЕНИЕ» 1976 6Ф7.3 У36 УДК 681.3,057 i 518.5 Уилкинсон, Райнш У36 Справочник алгоритмов на Линейная алгебра. Перевод с ред. д-ра техн. наук проф. Ю. «Машиностроение», 1976 г. 389 с. с ил. языке АЛГОЛ. английского. Под И. Топчеева. М., В книге приведены алгоритмы решения всех основных задач линей- ной алгебры. Соответствующие алгоритмы реализованы в виде проце- дур на алгоритмическом языке АЛГОЛ-60, сопровождаемых тесто- выми задачами. Особый интерес для специалистов по теории управле- ния представляют алгоритмы решения полной проблемы собственных значений для произвольных матриц. Книга предназначена для широкого круга научных работников, инженеров, использующих вычислительные машины для решения прикладных задач, а также для специалистов по программированию, аспирантов и студентов различных специальностей. U38 (01)-76 30502-013 13—76 6Ф7.3 Издательство «Машиностроение».. 1976 г. Оглавление Предисловие к русскому изданию .................. Предисловие к английскому изданию ................ Часть I Системы линейных уравнений, метод наименьших квадратов и линейное программирование Алгоритмы и их назначение .................... 1. Введение (13). 2. Список процедур (13). 3. Положи- тельно определенные матрицы (14). 4. Неположительные симметрические матрицы (16). 5. Произвольные матрицы (17). 6. Метод наименьших квадратов и псевдообращение матриц (18). 7. Задача линейного программирования (19) Алгоритм 1.1. Треугольное разложение положительно определенных сим- метрических матриц ........................ 1. Теоретические предпосылки (20). 2. Применение алго- ритма (23). 3. Список формальных параметров (24). 4. Про- граммы на языке АЛГОЛ (26). 5. Организация процедур и обозначения (33). 6. Оценка точности решения (35). 7. Примеры использования процедур и результаты их про- верки (36). Алгоритм 1.2. Итерационное уточнение решения системы уравнений с положительно определенной матрицей .............. 1. Теоретические предпосылки (40). 2. Применение алго- ритма (41). 3. Список формальных параметров (42). 4. Программы на языке АЛГОЛ (43). 5. Организация процедур и обозначения (46). 6. Оценка точности реше- ния (47). 7. Примеры использования процедур и резуль- таты их проверки (49). Алгоритм 1.3. Обращение положительно определенных матриц методом Гаусса .............................. 1. Теоретические предпосылки (52). 2. Применение алго- ритма (53). 3. Список формальных параметров (53). 4. Про- граммы на языке АЛГОЛ (53). 5. Организация процедур и обозначения (54). 6. Оценка точности решения (55). 7. Примеры использования процедур и результаты их проверки (55). Алгоритм 1.4. Симметричное треугольное разложение положительно определенных ленточных матриц ................. 56 1. Теоретические предпосылки (56). 2. Применение алго- ритма (57). 3. Список формальных параметров (57). 4. Про- граммы на языке АЛГОЛ (58). 5. Организация процедур и обозначения (59). 6. Оценка точности решения (60). 7. Примеры использования процедур и результаты их про- верки (61). Алгоритм 1.5. Метод сопряженных градиентов ............ 63 1. Теоретические предпосылки (63). 2. Применение алго- ритма (64), 3, Список формальных параметров (65). 4. Программы на языке АЛГОЛ (65). 5. Организация процедур и обозначения (66). 6. Оценка точности реше- ния (67). 7. Примеры использования процедур и резуль- таты их проверки (68). Алгоритм 1.6. Решение систем уравнений с ленточными матрицами и вычисление собственных векторов ленточных матриц ........ 71 1. Теоретические предпосылки (71). 2. Применение алго- ритма (74). 3. Список формальных параметров (75). 4. Про- граммы на языке АЛГОЛ (77). 5. Организация процедур и обозначения (85). 6. Оценка точности решения (86). 7. Примеры использования процедур и результаты их проверки (88). Алгоритм 1.7. Решение действительных уравнений ............ комплексных систем линейных 92 1. Теоретические предпосылки (92). 2. Применение алго- ритма (95). 3. Список формальных параметров (95). 4. Про- граммы на языке АЛГОЛ (97). 5. Организация процедур и обозначения (102). 6, Оценка точности решения (103). 7. Примеры использования процедур и результаты их про- верки (104). Алгоритм 1.8. Применение преобразований Хаусхолдера при реализации метода наименьших квадратов .................. 1. Теоретические предпосылки (107). 2. Применение алго- ритма (108). 3. Список формальных параметров (108). 4. Про- граммы на языке АЛГОЛ (109). 5. Организация процедур и обозначения (111). 6. Оценка точности решения (111). 7. Примеры использования процедур и результаты их про- верки (112). Алгоритм 1.9. Решение систем уравнений и проблема наименьших ква- дратов .............................. 1. Теоретические предпосылки (113). 2. Применение алго- ритма (115). 3. Список формальных параметров (115). 4. Про- граммы на языке АЛГОЛ (116). 5. Организация процедур и обозначения (121). 6. Оценка точности решения (122). 7. Примеры использования процедур и результаты их проверки (123), 107 113 Алгоритм 1.10. Разложение по сингулярным числам и проблема наи- меньших квадратов ........................ 1. Теоретические предпосылки (125). 2. Применение алго- ритма (129). 3. Список формальных параметров (131). 4. Про- граммы на языке АЛГОЛ (131). 5. Организация процедур и обозначения (136). 6. Оценка точности решения (137). 7. Примеры использования процедур и результаты их про- верки (137). Алгоритм 1.11. Реализация симплекс-метода с использованием тре- угольного разложения ...................... 1. Теоретические предпосылки (140). 2. Применение алго- ритма (147). 3. Список формальных параметров (148). 4. Про- граммы на языке АЛГОЛ (156). 5. Организация процедур и обозначения (164). 6. Оценка точности решения (165). 7. Примеры использования процедур и результаты их про- верки (167). Часть II Алгебраическая проблема собственных значений Алгоритмы и их назначение 1. Введение (173). 2. Список процедур (174). 3. Действи- тельные плотные симметрические матрицы (174). 4. Сим- метрические ленточные матрицы (176). 5. Одновременное определение собственных значений и векторов симметри- ческих редких матриц (177). 6. Обобщенная проблема соб- ственных значений для симметрических матриц (177). 7. Эр- митовы матрицы (178). 8. Действительные плотные несим- метрические матрицы (178). 9. Несимметрические ленточ- ные матрицы (180). 10. Плотные несимметрические матрицы с комплексными элементами (180). Алгоритм И.1. Метод Якоби для действительных симметрических матриц 1. Теоретические предпосылки (182). 2. Применение алго- ритма (183). 3. Список формальных параметров (183). 4. Про- граммы на 'языке АЛГОЛ (18)4. 5. Организация процедур и обозначения (185). 6. Оценка точности решения (186). 7. Примеры использования процедуры и результаты про- верки (187). Алгоритм 11.2. Преобразование Хаусхолдера для приведения симметри- ческой матрицы к трехдиагональной форме ............. 1. Теоретические предпосылки (190). 2. Применение алго- ритма (191). 3. Список формальных параметров (192). 4. Про- граммы на языке АЛГОЛ (193). 5. Организация процедур и обозначения (197). 6. Оценка точности решения (199). Г. Примеры использования процедур и результаты их про- верки (200). Алгоритм 11.3. Применение QR и QZ.-методов для определения собствен- ных значений и векторов симметрических матриц .......... 203 1. Теоретические предпосылки (203). 2. Применение алго- ритма (206). 3. Список формальных параметров (207). 4. Про- граммы на языке АЛГОЛ (207). 5. Организация процедур и обозначения (210). 6. Оценка точности решения (211). 7. Примеры использования процедур и результаты их про- верки (211). Алгоритм 11.4. QZ.-алгоритм с неявным сдвигом ........... 216 1. Теоретические предпосылки (216). 2. Применение алго- ритма (218). 3. Список формальных параметров (218). 4. Про- граммы на языке АЛГОЛ (219). 5. Организация процедур и обозначения (221). 6. Оценка точности решения (221). 7. Примеры использования процедур и результаты их про- верки (221). Алгоритм 11.5. Вычисление собственных значений симметрической трсх- диагональной матрицы методом деления отрезка пополам .... 22о 1. Теоретические предпосылки (223). 2. Применение алго- ритма (224). 3. Список формальных параметров (224). 4. Про- граммы на языке АЛГОЛ (224). 5. Организация процедур и обозначения (226). 6. Оценка точности решения (227). 7. Примеры использования процедур и резульгаты их про- верки (228). Алгоритм 11.6. 9/?-алгоритм в сочетании с методом Ньютона для вычи- сления отдельных собственных значений симметрических трехдиаго- нальных матриц ......................... 230 1. Теоретические предпосылки (230). 2. Применение алго- ритма (233). 3. Список формальных параметров (233). 4. Про- граммы на языке АЛГОЛ (234). 5. Организация процедур и обозначения (235). 6. Оценка точности решения (235). 7. Примеры использования процедур и результаты их про- верки (236). Алгоритм 11.7. 9/?-алгоритм для симметрических ленточных матриц . . . 238 1. Теоретические предпосылки (238). 2. Применение алго- ритма (238). 3. Список формальных параметров (238). 4. Про- граммы на языке АЛГОЛ (239). 5. Организация процедур и обозначения (241). 6. Оценка точности решения (242). 7. Примеры использования процедур и результаты их про- верки (242). Алгоритм 11.8. Приведение симметрической ленточной матрицы к трех- диагональной форме ....................... 244 1. Теоретические предпосылки (244). 2. Применение алго- ритма (246). 3. Список формальных параметров (246). 4. Про- граммы на языке АЛГОЛ (246). 5. Организация процедур и обозначения (247). 6. Оценка точности решения (248). 7. Примеры использования процедур и результаты их про- верки (249). Алгоритм 11.9. Метод одновременной итерации для симметрических матриц ............................. 1. Теоретические предпосылки (252). 2. Применение алго- ритма (252). 3. Список формальных параметров (252). 4. Про- граммы на языке АЛГОЛ (255). 5. Организация процедур и обозначения (260). 6. Оценка точности решения (263). 7. Примеры использования процедур и результаты их про- верки (264). Алгоритм 11.10. Решение полной проблемы собственных значений для систем уравнений общего вида с симметрическими матрицами ..... 1. Теоретические предпосылки (267). 2. Применение алго- ритма (268). 3. Список формальных параметров (268). 4. Про- граммы на языке АЛГОЛ (269). 5. Организация процедур и обозначения (272). 6. Оценка точности решения (273). 7. Примеры использования процедур и результаты их про- верки (274). Алгоритм 11.11. Масштабирование матриц при вычислении собствен- ных значений и векторов ..................... 1. Теоретические предпосылки (277). 2. Применение алго- ритма (280). 3. Список формальных параметров (280). 4. Про- граммы на языке АЛГОЛ (281). 5. Организация процедур и обозначения (282). 6. Оценка точности решения (283). 7. Примеры использования процедур и результаты их про- верки (283). Алгоритм 11.12. Решение проблемы собственных значений по методу Якоби с понижением нормы для действительных матриц ....... 1. Теоретические предпосылки (287). 2. Применение алго- ритма (288). 3. Список формальных параметров (289). 4. Про- граммы на языке АЛГОЛ (289). 5. Организация процедур и обозначения (292). 6. Оценка точности решения (293). 7. Примеры использования процедуры и результаты ее про- верки (293). Алгоритм 11.13. Приведение матриц общего вида к форме Хессенберга 1. Теоретические предпосылки (298). 2. Применение алго- ритма (300). 3. Список формальных параметров (300). 4. Про- граммы на языке АЛГОЛ (302). 5. Организация процедур и обозначения (307). 6. Оценка точности решения (308). 7. Примеры использования процедур и результаты их про- верки (309). Алгоритм 11.14. 9/?-алгоритм для вычисления собственных значений действительной матрицы в форме Хессенберга ............ 1. Теоретические предпосылки (316). 2. Применение алго- ритма (319). 3. Список формальных параметров (320). 4. Про- граммы на языке АЛГОЛ (320). о. Организация процедур и обозначения (322). 6. Оценка точности решения (323). 7. Примеры использования процедур и результаты их про- верки (323). Алгоритм 11.15. LR и 9/?-алгоритмы для вычисления собственных век- торов действительных и комплексных матриц ............ 1. Теоретические предпосылки (327). 2. Применение алго- ритма (428). 3. Список формальных параметров (329). 4. Про- граммы на языке АЛГОЛ (330). 5. Организация процедур и обозначения (339). 6. Оценка точности решения (340). 7. Примеры использования процедур и результаты их про- верки (342). Алгоритм 11.16. Модифицированный LR -алгоритм для вычисления соб- ственных значений комплексной матрицы в форме Хессенберга . . . 1. Теоретические предпосылки (348). 2. Применение алго- ритма (350). 3. Список формальных параметров (350). 4. Про- граммы на языке АЛГОЛ (351). 5. Организация процедур и обозначения (353). 6. Оценка точности решения (353). 7. Примеры использования процедур и результаты их про- верки (353). Алгоритм 11.17. Решение проблемы собственных значений по методу Якоби с понижением нормы для комплексных матриц ......... 1. Теоретические предпосылки (355). 2. Применение алго- ритма_(356). 3. Список формальных параметров (356). 4. Про- граммы на языке АЛГОЛ (357). 5. Организация процедур и обозначения (360). 6. Оценка точности решения (362). 7. Примеры использования процедур и результаты их про- верки (362). Алгоритм 11.18. Определение отдельных собственных векторов методом обратной итерации ....................... 1. Теоретические предпосылки (367). 2. Применение алго- ритма (369). 3. Список формальных параметров (370). 4. Про- граммы на языке АЛГОЛ (371). 5. Организация процедур и обозначения (380). 6. Оценка точности решения (381). 7. Примеры использования процедур и результаты их про- верки (382). 348 355 367 Список литературы ............ Литература, добавленная при переводе книги 385 ЗУб Предисловие к русскому изданию Предлагаемый вниманию читателя перевод «Справочника алгоритмов на языке АЛГОЛ. Линейная алгебра» Уилкинсона и Райнша будет с интересом встречен всеми, кто использует в своей научной и практической деятельности цифровые вычислительные машины. Справочник содержит 84 процедуры, напи- санные на языке АЛГОЛ-60, и охватывает широкий класс алгоритмов, связанных с решением систем линейных уравнений, псевдообращением матриц, вычислением собственных значений и собственных векторов. При описании каждого алгоритма исследована область применения, дано подробное описание процедуры на алго- ритмическом языке, а также оценка точности и результаты тестовой проверки. При разработке процедур значительное внимание было'уделено экономии памяти и обеспечению высокой скорости сходимости алгоритмов. Данный справочник существенно дополнит известные справочники алго- ритмов [102—106, 112]. Он позволит значительно повысить эффективность исполь- зования ЦВМ при решении прикладных задач оптимального управления, плани- рования и проектирования автоматизированных систем. Алгоритмы, приведен- ные в справочнике, можно использовать для построения специализированных рабочих программ анализа и синтеза автоматических систем, построения пере- ходных процессов, обработки информации, оценки эффективности автоматизи- рованных систем. С теоретическими основами решения систем линейных уравнений, методов оптимальной обработки информации, линейного программирования и обобщен- ной проблемой собственных значений читатель может ознакомиться в отечествен- ных книгах, приведенных в списке литературы, добавленной при переводе. Для читателей, желающих освоить приемы программирования на языке АЛГОЛ, мы также даем перечень соответствующих книг. Проверка значительной части алгоритмов показала, что без каких-либо модификаций они могут применяться на отечественных вычислительных машинах (БЭСМ-6, БЭСМ-4, М-220; ЕС ЭВМ), оснащенных трансляторами с языка АЛГОЛ-60. Выпуск данного справочника свидетельствует о том, что автоматизация про- граммирования достигла высокого уровня развития и позволяет реализовать пре- емственность алгоритмов, независимо от типа используемых ЦВМ. Книга рассчитана на научных работников и инженеров, использующих вычислительные машины для решения прикладных задач, а также для специа- листов по программированию, аспирантов и студентов различных специаль- ностей. /О. Топчеев Предисловие к английскому изданию Применение в различных странах единого машинного языка, каким является АЛГОЛ, позволяет создать универсальные вычислительные программы, которые могут быть использованы без каких-либо изменений программистами, работаю- щими на различных вычислительных машинах, снабженных транслятором с языка АЛГОЛ. Справочник содержит алгоритмы одного из разделов численного анализа. Выбор темы «Линейная алгебра» авторы считают наиболее целесообразным, поскольку соответствующие алгоритмы широко используются при численном решении задач; кроме того, они образуют некоторый законченный раздел. Приве- денные алгоритмы разделены на две группы: одна связана с решением систем линейных уравнений, вторая — с решением полной проблемы собственных зна- чений. Именно потому книга разделена на две части, каждой из которых предше- ствует введение, в котором проанализированы относительные достоинства алго- ритмов. Несмотря на простоту и ясность формулировок задач линейной алгебры, практика показала, что даже в этом случае подготовка рабочего алгоритма отни- мает достаточно много времени. При выборе рабочих алгоритмов отбирались те, которые обеспечивали, хотя и для ограниченной области применения, наилучшие решения, с высокой ско- ростью сходимости и минимальным объемом используемой памяти машины. Отсут- ствие в тексте того или иного из известных алгоритмов не означает, что его не сле- дует использовать; просто авторы не нашли оптимальной формы его реализации. Само содержание книги указывает на то, что это плод совместных усилий большого числа исследователей. Во многих случаях исходная основа алгоритма претерпевала изменения в течение длительного периода времени и многие мате- матики принимали участие в улучшении первоначального варианта. Мы бла- годарны всем, кто принимал непосредственное участие в этом издании. Публикация материала подобного содержания связана с теми особенно- стями, что выбранные рабочие алгоритмы должны быть полезными в течение длительного периода времени. Авторы убеждены, что создание работоспособных алгоритмов — процесс непрерывный, и они будут благодарны всем, кто укажет на ошибки и трудности, встретившиеся при практической работе с предложен- ными алгоритмами. Уилкинсон, Райнш Часть I Системы линейных уравнений, метод наименьших квадратов и линейное программирование Алгоритмы и их назначение 1. Введение Системы линейных алгебраических уравнений возникают при решении большого числа практических задач. Они порождают три типа вычислительных проблем: 1) решение невырожденных систем линейных уравнений, обращение матриц и вычисление определителей; 2) решение произвольных систем уравнений методом наименьших квадратов и псевдообращение матриц; 3) решение задач линейного программирования. Алгоритмы линейного программирования и метода наименьших квадратов имеют много общего с алгоритмами решения стандартных задач линейной алгебры, и поэтому они объединены вместе. Почти все алгоритмы решения систем уравнений п-го порядка, приведенные в книге, начинаются с разложения исходной матрицы А на треугольные, что позво- ляет очень просто вычислить определитель матрицы А. Задав в качестве правой части системы уравнений единичную матрицу размера яХ п, можно с помощью тех же алгоритмов вычислить матрицу, обратную исходной. При составлении Справочника было решено не включать итерационные методы решения систем линейных алгебраических уравнений, поскольку их в основном используют при решении дифференциальных уравнений в частных производных. Заметим также, что эффективность того или иного алгоритма зависит от свойств конкретной вычислительной машины1, например, объема памяти. В алгоритмах, приведенных в настоящей книге, нет указаний по использованию внешних накопителей — магнитных лент или дисков. Последнее обусловлено тем, что распределение памяти при решении практических задач может оказаться совсем иным, чем это представлялось составителю алгоритма. Например, если требуется итерационное уточнение решения системы линейных уравнений, то исходная матрица должна сохраняться. На практике такое условие легко выполнить, используя внешний накопитель. 2. Список процедур Для упрощения поиска все процедуры перечислены ниже в алфавитном порядке с указанием тех алгоритмов, в которых они описаны. Помимо процедур, связанных с решением задач линейной алгебры, в список включены вспомога- тельные процедуры вычисления скалярного произведения, арифметических опе- рации над комплексными числами и т. д. ' Понятие эффективности алгоритма определяется как областью его применения; так быстродействием и объемом требуемой памяти; поэтому авторы и указывают на воз- ыожную зависимость эффективности от типа машины {Прим. ред. пер-). 13 Список процедур осе inverse, асе solve — 1.2 bande.i I, bansol 1 — 1.6 bandel 2, banwl 2—1.6 cabs, cdiv, csqrt — 11.13 cg — 1.5 chobanddet, chobandsol — 1.4 choldel I, cholsol 1 cholinversion !— 1.1, 1.2 choldet 2, cholsol 2, cholinversion 2—1.1 compdel, compsol, ex ace solve — 1.7 gjdef I, gjdef 2-1.3 innerprod, ex innerprod — 1.2, 1.7 least squares solution — 1.8 Ip- 1.11 minfit, svd— 1.10 wtho 1, orlho 2, oriholin /. orlholin 2— 1.9 svd, minfit — 1.10 symdel, symsol, syminversion— 1.1 unsym ace solve, unsymdet, unsymsol — 1.7 3. Положительно определенные матрицы При построении вычислительных процедур решения систем линейных урав- нений следует различать три типа систем; с положительно определенными симметрическими матрицами; с симметрическими, но неположительными матрицами и, наконец. с произвольными матрицами. Рассмотрение начнем с наиболее простого случая. Как правило, алгоритмы для разложения положительно определенных симметрических матриц связаны с представлением исходной матрицы либо в виде А == LL1, где L — нижняя треугольная матрица (схема Холецкого), либо в виде А == LDL7, где L—нижняя треугольная матрица с единичной диагональю, a D — диагональная матрица с положительными элементами. Эти алгоритмы обе- спечивают устойчивый вычислительный процесс и требуют выполнения малого числа арифметических операций. Существенных преимуществ в представлении исходной матрицы в виде LL7 или LDL7 нет, но в первом случае при решении /'1Л/"Т<1»*1-1 17П01Э1-ЮТТ11Й llQ/^^vi-i-nTi»»" /^»,/ч-> ^..—— >.—— —-—-----— --- -. квадратного корня системы уравнений необходима операция извлечения а во втором она отсутствует. Плотные положительно определенные матрицы Схема Холецкого. Можно выделить две основные группы процедур, связанных с разложением плотных положительно определенных матриц. В первую группу входят процедуры: choldet I, cholsol 1 и cholinversion I; во вторую — choldet 2, cholsol 2 и cholinversion 2. В процедурах второй группы достигнута значительная экономия памяти за счет работы лишь с массивом из -„- п (п 4- 1) элементов ма- трицы А, расположенных на диагонали и ниже. Эти элементы в процессе вычисле- ний заменяются элементами матрицы L. Для реализации процедур второй группы необходимы большие затраты машинного времени, чем для процедур первой группы. Поэтому процедуры вчорой группы используют лишь в тех случаях, когда требуется экономия памяти. * Поэтому в отечественной литературе можно встретить название <метод квадрат- ных корней» (см.; например Демидович Б. П и Марон И. А. «Основы вычислительной математики». М.; »Наука»; 1970). 14 С помощью процедур choldet в обоих вариантах выполняется разложение матрицы А на произведение вида LL7 и одновременно вычисляется ее определи- тель. После выполнения процедуры choldet решение системы уравнений можно получить с помощью процедуры cholsol, а матрицу, обратную А, — используя процедуру cholinversion. Хотя алгоритмы, реализующие схему Холецкого, сами по себе корректны и их точность сравнима с машинной, результаты решения системы уравнении или обращения матрицы могут иметь низкую точность из-за плохой обусловлен- ности исходной матрицы. Правда, точность результатов .и в этом случае может быть увеличена в рамках схемы Холецкого. Для этого надо воспользоваться процедурами асе solve» асе inverse, в которых реализован итерационный процесс уточнения решения системы уравнений или результата обращения до машинной точности, при условии, что исходная матрица А не является вырожденной при задании ее коэффициентов с машинной точностью. Процедура итерационного уточнения в конечном счете определяет решение системы уравнений в виде век- — II х ~ х II со тора х, удовлетворяющего условию ——————^ k-macheps, где macheps— ма- Цх||„ шинная точность1; k—коэффициент порядка 1. Для достижения такой высокой точности в этих алгоритмах использована вспомогательная процедура innerprod, которая позволяет вычислить с двойной точностью скалярное произведение двух векторов, заданных с обычной точностью. Эта процедура должна быть написана в кодах машины. Ради полноты изложения приведены также модифицированные варианты процедур choldet 1 и cholsol 1, в которых все скалярные произведения вычисля- ются с удвоенной точностью. Если реализация накопления скалярного произве- дения на вычислительной машине затруднительна, то применяют обычные варианты процедур choldet 1 и cholsol 1. Однако накопление скалярного произве- дения при вычислении невязок обязательно и является самой существенной осо- бенностью процесса итерационного уточнения. Процедуры асе inverse и асе solve играют важную роль при оценке точности решения. Действительно, если x'r— исходное решение, а х'2'— его первое уточ- нение, то величина х'2'—x'r дает хорошую оценку влияния ошибки в последнем разряде исходных данных (изменение в последнем разряде зависит от длины машин- ного слова) и, таким образом, позволяет получить оценку решения при неточном задании матрицы А. Если операция вычисления скалярного произведения с удвоен- ной точностью на машине возможна, то целесообразно этим воспользоваться, тем более, что такая процедура не требует больших затрат времени. В заключе- ние отметим, что для итерационного уточнения необходимо запоминать матрицу А и векторы правых частей исходной системы уравнений. Представление исходной матрицы в виде LDL1. Разложение исходной ма- трицы А на произведение вида LDL1 реализовано в серии процедур symdet, symsol и syminversion. Требуемое разложение с одновременным вв1числением det А выполняет процедура symdet. После выполнения процедуры symdet можно найти решение исходной системы уравнений и вычислить матрицу, обратную А, с помощью процедур symsol и syminversion соответственно. Следует заметить, что, используя разложение вида LDL^i можно составить алгоритмы, аналогичные choldet 2, cholsol 2 и cholinversion 2, и записать процедуры асе solve и асе inverse с применением алгоритмов symdet и symsol. Существенной разницы между соот- ветствующими алгоритмами, в которых использованы разложения вида LL или LDL7, нет. Вычисление обратной матрицы методом Гаусса. Для обращения положительно определенной матрицы с записью результата на прежнем месте можно применить ' Машинная точность — характеристика вычислительной машины, она определяется числом разрядов машины; так, например, машинное слово, имеющее bj разряда, ооеспе- чивает точность 2-" я» Ю--". 15 процедуры gjdej 1 и gjdef 2. Во второй из них использован объем памяти, равный памяти процедуры choldet 2. Однако алгоритм Гаусса не следует применять для отыскания решения системы уравнений, поскольку он уступает по своей эффек- тивности комбинации процедур choldet—cholsol. При обращении матриц его точ- ность сравнима с точностью алгоритмов cholinversion I, 2. Положительно определенные ленточные матрицы Многие физические задачи описываются системами уравнений с положительно определенными ленточными матрицами. Для таких матриц схема Холецкого осо- бенно эффективна, поскольку в процессе вычислений сохраняется ленточная структура матрицы. Процедура chobanddet выполняет треугольное разложение и вычисляет определитель исходной матрицы А; последующее применение про- цедуры chobandsol позволяет найти решение исходной системы уравнений. Обра- щения таких матриц обычно не требуется (причем отметим, что обратная матрица будет полной, хотя исходная и является ленточной), и поэтому процедура вида choband inverse в справочнике отсутствует. Уточнение решения можно получить с помощью итерационной процедуры, аналогичной асе solve. Другой метод реше- ния уравнений с ленточными матрицами будет рассмотрен ниже. Редкие положительно определенные матрицы специальной формы Многие задачи математической физики приводят к положительно определен- ным матрицам с небольшим числом ненулевых элементов, расположенных в раз- ных местах матрицы. Эти элементы (а^) можно вычислить как функцию пара- метров (', /'. Последнее имеет место при переходе от дифференциальных уравнений в частных производных к разностным. Если исходная матрица является редкой, а функции от переменных г и / достаточно просты, то процедура cg, основанная на методе сопряженного градиента, может оказаться весьма эффективным алго- ритмом для решения таких систем уравнений. Метод сопряженного градиента иногда рассматривают как итерационный, однако при точных вычислениях он дает точное решение за конечное число шагов. Именно поэтому соответствующий алгоритм и включен в Справочник, хотя чисто итерационные методы в нем не рассматриваются. Процедура cg весьма экономно использует память машины, так как запоми- нать матрицу А нет необходимости. Вместо этого нужно лишь вычислять произ- ведение Ах при заданном векторе х. Если матрица А имеет высокий порядок, то при- менение метода сопряженного градиента (или какого-либо другого итерационного метода) оказывается единственно возможным способом решить задачу, поскольку при больших размерах матрицы ее не удается запомнить, даже если она является узкой ленточной матрицей. Это имеет место при аппроксимации уравнений в част- ных производных. Как правило, получаются узкие ленточные матрицы, причем большинство элементов равно нулю, и поэтому процедура вычисления вектора Ах оказывается достаточно простой. Если точность процедуры cg и может оказаться выше, чем в схеме Холец- кого, то по скорости счета она всегда будет уступать последней, поскольку при каждом изменении вектора правой части требуется повторить заново процедуру cg. В алгоритмах, реализующих схему Холецкого, разложение матрицы А выпол- няется только 1 раз независимо от числа правых частей. 4. Неположительные симметрические матрицы В Справочнике нет алгоритмов, в которых было бы использовано свойство симметрии при решении систем уравнений с плотными симметрическими матри- цами. Поэтому последние следует решать с помощью общих процедур для плот- ных матриц. Особо выделим процедуры ba.nd.et I, bansol I и bandet 2, bansol 2, которые можно применить для решения систем с симметрическими ленточными матрицами. Однако в процессе выполнения алгоритма симметрия матрицы нару- 16 шается. Процедуры bandet I или 2 можно использовать для разложения А и вычис- ления del А, а для решения системы уравнений — процедуры bansol I, 2. Первое множество процедур заметно эффективнее второго, но если А — симметрическая матрица, то процедура bandet 2 определяет также все положительные собственные значения. Поэтому ее можно использовать для отыскания всех собственных .вне-- чений, больших заданного р, если в качестве исходной принять матрицу А—р1 5. Произвольные матрицы Алгоритмы для решения систем уравнений с произвольными матрицами сложнее описанных выше для систем уравнений с положительно-определенными симметрическими матрицами. В частности, значительные трудности связаны с про- блемой масштабирования, которая до сих пор не имеет удовлетворительного реше- ния. Поэтому процедуры, приведенные в Справочнике,—эго компромисс между желаемым и достигнутым. Плотные матрицы (действительные и комплексные) Все процедуры для решения систем уравнений с произвольными матрицами основаны на представлении исходной матрицы А в виде произведения LU, где L — нижняя треугольная и U — верхняя треугольная с единичной диагональю. Для разложения матрицы использован алгоритм с частичным выбором главного эле- мента. Он напоминает алгоритм исключения Гаусса, который, правда, приводит к нижней треугольной матрице с единицами по диагонали. Процедура unsymde t позволяет разложить действительную плотную матрицу А в произведение LU и одновременно вычислить det А. После этого с помощью про- цедуры unsymsol можно найти решение для любого числа правых частей. Про- цедуры compdet и co'npsol позволяют решить те же задачи для комплексной матрицы А. В тех случаях, когда матрица А плохо обусловлена, процедуры unsym асе solue и сх асе solve позволяют уточнить решение, используя в процессе итерацион- ного уточнения треугольное разложение исходной матрицы, найденное в про- цедурах unsymdet или compdet соответственно. Если матрица А не является вырожденной при задании ее коэффициентов с машинной точностью, то процедура итерационного уточнения будет в конце концов давать решение х с машинной 1|х-х||„ точностью, т. е. —T~ij——^ k-macheps, гда macheps— машинная точность и k — " "J^ константа, близкая к единице. Отметим, что «малые» компоненты вектора х, как правило, имеют низкую относительную точность. Первый шаг итерационного уточнения позволяет полу- чить оценку чувствительности решения системы уравнений. Если векторы Xе1' и х1-"— первые два уточненных решения, то разность х'"'—х11' дает хорошую оценку точности решения от ошибки задания исходных данных. Несимметрические ленточные матрицы Процедуры bandet i и bansol I предназначены для проведения вычислений с действительными несимметрическими матрицами ленточной формы. Они весьма эффективны при работе с матрицами, имеющими различное число ненулевых эле- ментов по обе стороны от главной диагонали. В этих процедурах использован алгоритм с частичным выбором главного элемента для треугольного разложения исходной матрицы. Данные процедуры применимы и для симметрических лен- точных матриц. Процедура bandet I предназначена для разложения исходной матрицы и вычисления ее определителя; после этого для решения системы урав- нений при произвольном числе правых частей следует применять процедуру bansol I. 17 Процедуры bandet 2 и bansol 2 также можно использовать для решения систем уравнений с несимметрическими матрицами, но с меньшей эффективностью; они- предназначены в основном для работы с симметрическими ленточными матрицами. 6. Метод наименьших квадратов и псевдообращение матриц Метод наименьших квадратов. В Справочнике приведены две процедуры для отыскания решений системы уравнений с матрицей А размера тХп, когда т ^ п и ранг А = я, методом наименьших квадратов. Первая из них — least squares solution связана с определением ортогональной матрицы Q такой, что QA = R, где R — верхняя треугольная матрица размера тХп', А — матрица А с переставленными столбцами. Эта перестановка зависит от выбора главного элемента. Матрица Q представляет собой произведение п т *' элементарных ортогональных матриц вида 1—2ww , где || w |[з = 1. Вектор х гарантирует минимум нормы ||Ь—А х Цз и является решением- уравнения Rx = == Qb. Это уравнение может иметь произвольное число правых частей. Если в каче- стве правой части взять единичную матрицу размера тХт, то решением будет матрица (А^)" А" размера пХт. Ее называют псевдообратной матрицей, и она, по предположению, имеет ранг, равный п. Процедура включает итерационный процесс для уточнения полученного реше- ния. Этот процесс аналогичен описанному ранее для систем линейных уравнений п-го порядка. Для совместной системы уравнений, т. е. системы, для которой min [\ Ь—Ax |[g = 0, процесс уточнения приводит к решению, определенному с машинной точностью, но если система уравнений несовместна, то окончательная точность может оказаться недостаточной. Метод наименьших квадратов часто применяют для отыскания решений несовместных систем уравнений. В таких случаях процедура уточнения не играет существенной роли. Процедура least squares solution применима и в том случае, когда т == п; при этом определяют решение невырожденной системы линейных уравнений. Именно поэтому ее можно использовать при обращении матриц. Хотя ортого- нальное разложение имеет некоторые преимущества перед разложением вида LU. однако, что касается численной устойчивости, то этого преимущества вряд ли достаточно, чтобы оправдать применение процедуры least squares solution вместо процедур unsymdet и unsymsol. Наряду с алгоритмом least squares solution в Справочнике приведены алго- ритмы, которые используют ортогональные преобразования исходной матрицы А к треугольной форме, опираясь на процедуру ортогонализации Грама—Шмидта. При этом вычисляют ортогональный базис пространства, порожденного столб- цами матрицы А, и каждый вектор правых частей Ь выражают через этот базис и вектор невязки г, ортогональный ему. Такая факторизация реализована в четырех процедурах. Процедура ortho- tin 1 определяет решение по минимуму наименьших квадратов для матричного уравнения АХ = В, где В состоит из k векторов правой части. Процедуру ortholln 2 используют лишь тогда, когда имеется единственная правая часть. Процедура ortho 1 позволяет осуществить обращение матрицы размера пХ п. В этих процедурах предусмотрено уточнение полученного решения. В процедуре ortho 2 реализовано обращение матрицы размера пХп без уточнения решения, что позволило сэкономить память за счет размещения обратной матрицы на месте исходной. Основные характеристики этих процедур сравнимы со свойствами процедуры least squares solution. В частности, итерационное уточнение не всегда дает точное решение для несовместных систем уравнений яг-го порядка с п неиз- вестными. Обобщенный метод наименьших квадратов и псевдообращение матриц. Если матрица А размера тХ п имеет ранг, меньший п, то решение по методу наи- меньших квадратов не является единственным. В этом случае очень часто ищут решение в виде вектора х, для которого норма || х ||з имеет минимум. Это решение можно записать следующим образом: х = А*Ь, где А* — матрица, псевдо- 18 обратная А. Для решения этой задачи существует два алгоритма, в которых использовано разложение исходной матрицы А по сингулярным числам, т. е. осуществлено приведение матрицы А к диагональной форме. Если т ^ п, то это разложение можно записать в виде А =s USV7. где L)1!! = V1 V == VV^^ == In, S = tliag ((TI, . . ., Од) и о, — сингулярные числа матрицы A (U — ма- трица размера тХ п со взаимно ортогональными столбцами и V — ортогональная матрица). Соответствующее выражение для псевдообрашой матрицы имеет вид VS"^ , где S" = diag (of, . . ., о^) и ^ f о7'1, если о, > 0; о, — •J ( 0, если 0{ == 0. Разложение по сингулярным числам является самым надежным методом определения ранга матрицы А, элементы которой подвержены ошибкам. Процедура svd позволяет вычислить множество сингулярных чисел, а также матрицы LI и V, а процедура minfit — решение А* Ь для каждой из р заданных правых частей. Эти две процедуры можно широко использовать для отыскания решений по методу наименьших квадратов и вычислений псевдообратных матриц. 7. Задача линейного программирования Исследования в области линейного программирования до недавнего времени слабо увязывались с проблемами'линейной алгебры. Это отрицательно сказыва- лось на исследованиях, поскольку оба направления имеют много общего. Поэтому включение в настоящий Справочник процедуры линейного программирования обратит внимание читателей на эту общность и позволит раскрыть взаимосвязь между данными областями. Процедура 1р предназначается для решения следующей задачи: минимизи- ровать выражение Со + с.щХ-т + • • • + c-i.c.1 -г- ел + • • • + Wn при наличии следующих ограничений: /.' л-' + S aikxk =bi' '= 1,2, .... от; k=-l где /+ и /° — непересекающиеся множества из множества {—т, . . ., —1, 1, . ... п]. Процедура является одной из разновидностей симплекс-метода, основанного на приведении к треугольной форме текущего базиса А./ согласно уравнению LA.J = R, где R — невырожденная верхняя треугольная матрица. Замечено, что эффективность решения задач линейного программирования существенно зависит от типа машины. Для этих задач характерны большие мас- сивы перерабатываемой информации, что может привести к необходимости исполь- зовать внешние накопители. Для преодоления этих трудностей предложены две вспомогательные процедуры ар и р, первая из которых обрабатывает матрицу А, а вторая выполняет операции с матрицами R и I- и различными вариантами пра- вых частей В. Приведены упрощенные варианты процедур ар и р для тех случаев, когда матрицы A, R, L, В размещены в стандартных прямоугольных массивах. 19 Алгоритм I.I Треугольное разложение положительно определенных симметрических матриц 1. Теоретические предпосылки Все приведенные ниже варианты треугольного разложения основаны на сле- дующей теореме Холецкого [25]. Если А — симметрическая положительно опре- деленная матрица, то существует действительная невырожденная нижняя тре- угольная матрица L, такая, что т » LL' (1) Более того, если диагональные элементы матрицы L выбраны положитель- ными, то разложение единственно. Элементы матрицы L можно определять по строкам или по столбцам, прирав- нивая соответствующие элементы матриц в выражении (1). Если матрицу вычис- ляют по строкам, то для элементов г-ой строки справедливы следующие соотно- шения: k=\ или /-1 i/ — ^ liklft k=\ (2) Ч, = is likhk = = "It, k=\ или , 1-1 \1/2 hi- = "u- -s ^ . \ t=l (3) Таким образом, для реализации метода необходимо га действий извлечения 113 квадратного корня и приблизительно -„- умножении. Существует другой вариант разложения, в котором удается избежать извле- чения квадратного корня. Для этого определим новую матрицу L следующим образом: L=Ldiag (/„•), (4) где L — матрица, вычисленная ранее по схеме Холецкого. Матрица L существует (поскольку все элементы 1ц положительны) и является нижней треугольной матрицей с единицами по диагонали. При этом справедливо следующее соотно- шение: А= LL^Ldiag^^diag^^L^Ldiag (l^U^LDU, (5) .где n — положительно определенная диагональная матрица. 20 Такое разложение исходной матрицы можно выполнить за n шагов, причем на 1-ом шаге определяют 1-ую строку матрицы L и г-ый элемент d, матрицы D. Соответствующие выражения для нахождения этих элементов имеют вид k=i или (6) ^/^у ^ о»/— ]^ ~likdii~ljk, == 1, t=i S '^kdSik^- йц. k^.1 или t'—i di == ац — ^ 'likdkhk- (7) Алгоритм, представленный в таком виде, требует выполнения вдвое большего числа умножений, чем в разложении Холецкого. Однако если ввести вспомогатель- ные переменные а„, удовлетворяющие условию a,,=~li,dj, (8) то выражения (6) и (7) примут .вид /-1 V^ •" ', . , . , (9) •i, = "I/ —— ^J Cliklik, I •= 1, • . • , I —— 1; (10) t--=l i -1 ri( ^ ац— ^ afklik; k-:l поэчому сначала вычисляют вспомогательные величины ас,, а затем их используют для определения величин d, и 1ц. Заметим также, что для определения элемен- тов а;, в 1-ой строке не требуется знания элементов предыдущих строк. Коли- п3 чество операции умножения в этом варианте составляет приблизительно -,—» но не требуется операций извлечения квадратного корня. При любом варианте разложения матрицы А можно одновременно вычислить определитель матрицы в соответствии с выражениями то П la '(•=1 del А == detL.det L1 det A =-- dei L.det D det U = П dr. i==\ Можно также найти решение системы уравнений вида Ах .--= Ь. Действительно, если ввести два новых уравнения Ly-b; \ !Л=у,| Ах = Ll/x = Ly = b. (") (12) (13) (14) (15) 21 Соотношения (14) позволяют записать решение уравнения (13) в следующем виде: «•-1 ьi— S ^kVk k=i . , Vi = ———T~.————• t = 1, . . . , n; (16) У1 Xi =• • S lklxk k=i+l la , i = n, (17) что связано с выполнением n2 операций умножения и 2/г операций деления. Ана- логично, если то (18) (19) Ly == Ь; Lтx=D-ly, Ax = LDL ^x = LDD- 'у = b. Уравнения (18) можно решить, последовательно вычисляя величины 1-1 Vi = bi — ^ likVk, t = 1, . . . . n; ' (20) k==l Vl Xi. == • S ~1^' '=". .. • . 1- (21) ft=t+l При этом необходимо выполнить n2 операций умножения и лишь п операций деления. Обращение матрицы можно заменить решением п уравнений вида Ах=ег, t= 1, . . . , п; (22*) причем решению с t-ой правой частью соответствует (-ый столбец матрицы А~1. Однако алгоритм упростится, если в полной мере использовать свойство симме- трии А и особый вид правых частей в уравнении (22). Заметим, что А-^^-^-^^^-^-Ч-1, (23) где L~1 — нижняя треугольная матрица. Обозначая матрицу L"1 через Р, запишем 1 Р»=7~! 4i /-i ^J IjkPkt pi, = — •—————, / == f + 1, . . . , n, При вычислении элементов матрицы А как элементов матрицы Р Р (сим- метричная матрица) необходимо определить лишь элементы, расположенные ниже главной диагонали, и для их вычисления воспользоваться формулой (24) (25) (А-1)у, = (Р7?),, = ^ pkjpki. •}=- i, . . . , п. k=j (26)* * е^ — единичный вектор; 1-я компонента которого равна 1, а остальные равны О (Прим. ред. пер.). •* Запись X..,. используют для обозначения элемента матрицы Х (Прим. ред. пер,). 22 Заметим, что после того, как вычислен элемент (А"'1),,, величины р,, далее не используются. Аналогично, если нижнюю треугольную матрицу L"1 с единицами на диаго- нали обозначить через Qi, то будут справедливы следующие соотношения: ?«=1; (27) /-1 9/i == — ^ likqki, !==i+1' k=i ^•'-ST-' '=i- ••••"; <•=/ (28) (29) и здесь после вычисления элемента (А'1)/; знание величины дц более не требуется. Что касается числа арифметических операций, то при использовании разло- жения по схеме Холецкого их значительно больше, но проще организованы сами вычисления. Следует отметить, что на практике при обращении матриц исполь- зуют оба способа треугольного разложения исходной матрицы. 2. Применение алгоритма Разложение Холецкого, дающее представление исходной матрицы в виде А = LL , можно использовать только при расчетах с положительно определен- ными матрицами. Если ненулевые элементы матрицы А располагаются вблизи главной диагонали, то следует использовать процедуры chobanddet и chobandsol [47]. Если матрица А симметрическая, но не является положительно определен- ной, то разложений вида (1) или (5) может не существовать. Примером такой матрицы является матрица 1: ;]• Заметим, что невозможность разложения не означает, что матрица А плохо обусловлена. Если матрица А неотрицательная, то можно найти действительную матрицу L, удовлетворяющую соотношению (1), но она будет вырожденной и, вообще говоря, неединственной. Например, справедливо разложение ГО 01 f0 01 ГО х\ [о М. J'[o ,J (э" при условии )? + у2 = 2. Аналогично, но уже для всех х, справедливо разло- жение ;][° 2][; •О О^П 01ГО 1Г1 х1 О 2 \х l|| 2||0 l| (32) Если матрица А имеет хотя бы одно отрицательное собственное значение, то не существует действительной, но может существовать комплексная матрица L, удовлетворяющая соотношению (1). Треугольному разложению с комплексными матрицами можно поставить в соответствие треугольное разложение с действи- тельными матрицами L и D при условии, что диагональные элементы матрицы D отрицательны. При этом справедливо представление (33) 23 ч поскольку матрица D имеет отрицательные компоненты, то матрица LD- == L имеет столбцы, состоящие из чисто мнимых элементов. Существование разло- жения вида LDL1 над полем вещественных чисел является полезной особенностью для рассматриваемых матриц. Однако следует отметить, что численная устойчи- вость треугольного разложения гарантируется только для случая положительно определенной матрицы. Ниже приведены две группы алгоритмов, в которых использована схема Холецкого. Каждая группа состоит из трех процедур. Для процедур первой группы характерно сохранение исходной матрицы, в то время как в процедурах второй группы вычисляемую матрицу L записывают на место А. Именно поэтому требуемый объем памяти для процедур второй группы вдвое меньше, чем для первой, но отсутствие исходной матрицы не позволяет вычислить невязки и выпол- нить итерационное уточнение полученного решения системы уравнений или результата обращения матрицы. Разложение вида LDl^ реализовано в третьей группе процедур, причем здесь предусмотрено запоминание исходной матрицы А. Остальные особенности следуют из анализа самих процедур. Первые две группы алгоритмов реализуют схему Холецкого и включают следующие процедуры: choldet — процедура вычисления матрицы L и определителя det A; cholsol—процедура определения множества решений уравнения Ах= Ь при задании г правых частей. Процедуре cholsol всегда должна предшествовать процедура choldet. Процедуру cholsol можно применять произвольное число раз после единственного обращения к процедуре choldet; cholinversion — процедура обращения матрицы. Ее можно использовать независимо от процедур choldet и cholsol и следует применять в тех случаях, когда необходима подробная информация о деталях процесса обращения. 3. Список формальных параметров Входные параметры процедуры choldet 1: п — порядок матрицы А; а — массив элементов положительно определенной матрицы А. Его следует задать, как описано в п. 5. Выходные параметры процедуры choldeU: а — массив элементов нижней треугольной матрицы L в разложении Холец- кого; его размещение в памяти вычислительной машины списано в п. 5; р — одномерный массив величин, обратных диагональным элемента»] матрицы L; d\, d2 —два числа для записи вычисленного значения det А; fail — выходная метка, если матрица А не удовлетворяет условию положи- тельной определенности, возможно, из-за ошибок округления. Входные параметры процедуры cholsoll: п — порядок матрицы А; г — число, характеризующее количество правых частей, для которых необ- ходимо решить систему Ах == Ь; а — массив элементов нижней треугольной матрицы L, являющейся резуль- татом разложения исходной матрицы А с помощью процедуры choldetl; р — массив величин, обратных диагональным элементам матрицы L, опре- деленный с помощью процедуры choldet]; Ь — массив элементов, состоящий из т векторов правых частей исходной) уравнения (матрица размера пХг). Выходные параметры процедуры cholsoll: х — массив элементов, состоящий из г векторов решений системы уравнении Ах = Ь (матрица размера пХг}. Входные параметры процедуры cholinversion}: п — порядок матрицы А; 24 а — массив элементов положительно определенной матрицы А. Его следует задать, как описано в п. 5. Выходные параметры процедуры cholin version/: а — массив элементов матрицы А. размещение которого описано в п. 5; fail — выходная метка, если матрица А не удовлетворяет условию положи- тельной определенности, возможно, из-за ошибок округления. Входные параметры процедуры choldet2: п — порядок матрицы 'А; а — массив элементов положительно определенной матрицы А. Его следует вадагь в соответствии с рекомендациями п. 5. Выходные параметры процедуры choldet2: а — массив элементов нижней треугольной матрицы L, являющейся резуль- татом разложения матрицы А. Его размещение описано в п. 5; d/, d2 — два числа, позволяющие вычислить det А; mil — выходная метка, если матрица А не удовлетворяет условию положи- тельной определенности, возможно, из-за ошибок округления. Входные параметры процедуры cholsol2: п — порядок матрицы А; г — число правых частей, для которых необходимо решить систему Ах == Ь; а — массив элементов нижней треугольной матрицы L, являющийся резуль- татом разложения исходной матрицы А с помощью процедуры choldet2: Ь — массив элементов, состоящий из г векторов правых частей исходной системы уравнений (матрица размера пХг). Выходные параметры процедуры cholsol2: Ь •— массив элементов, состоящий из г векторов решений системы уравнений Ах =- b (матрица размера пХг). Входные параметры процедуры, cholin version2: п — порядок матрицы А; а — массив элементов положительно определенной матрицы А. Его следует задать в соответствии с рекомендациями п. 5. Выходные параметры процедуры cholinversion2: а — массив элементов матрицы А"1, размещение которого описано в п. 5; fail — выходная метка, если матрица не удовлетворяет условию положи- тельной определенности, возможно, из-за ошибок округления. Входные параметры процедуры symdet: п — порядок матрицы А; а — массив элементов положительно определенной матрицы А. Его следует задать в соответствии с рекомендациями п. 5. Выходные параметры процедуры symdet: а — массив элементов нижней треугольной матрицы L в разложении вида LDL. Его размещение описано в п. 5; р — одномерный массив величин, обратных диагональным элементам ма- трицы D; dl, d2—два числа, позволяющие вычислить det А; fail — выходная метка, если матрица А оказывается вырожденной, возможно, из-за ошибок округления. Входные параметры процедуры symsol: п — порядок матрицы А; г — число, характеризующее количество правых частей, для которых необ- ходимо решить систему Ах = b; а — элементы нижней треугольной матрицы L, являющейся результатом процедуры symdet; р — массив величин, обратных диагональным элементам матрицы D, опре- деленной с помощью процедуры symdet; Ь — массив элементов, состоящий из г векторов правых частей исходной системы уравнений (матрица размера пХг}. Выходные параметры процедуры symsol: < — массив элементов, состоящий из г векторов решения 1;исгемы уравнений Ах =г b. 25 Входные параметры процедуры symin version: п — порядок матрицы А; а — массив элементов положительно определенной матрицы А. Его разме- щение описано в п. 5. Выходные параметры процедуры syminversion'. а — массив элементов матрицы А~1; его размещение описано в п. 5; fail — выходная метка, если матрица А оказывается вырожденной, возможно, из-за ошибок округления. 4. Программы на языке АЛГОЛ Первая группа процедур procedure choldet I (n) irons: {a} result: (p, dl, d2) exit: (fail); value n; integer d2, n; real dl; array a, p; label fail; comment — элементы положительно определенной симметрической матрицы А размещены в массиве a [i, j], i = 1, . . , n, / = 1, . . ., n (см. п. 5). Процедура реализует разложение по схеме Холецкого вида А = LU, где U — матрица, транспонированная L. Элементы нижней треугольной матрицы L, за исключением диагональных, накапливаются в части массива а (см. п. 5). Величины, обратные диагональным элементам, размещены в массиве p [i], i == 1, . . ., п. Матрица А сохранена, поэтому в дальнейшем возможно осуществить итерационное уточне- ние. При вычислении определителя матрицы А использовано представление diX 2 f d2. Выход из процедуры по метке fail происходит в том случае, если из-за ошибок округления вычисленные значения квадратов диагональных элементов матрицы L оказываются неположительными; begin integer i, j, k; real x, dl-.^l; d2: = 0; fon : == 1 step 1 until n do for; : = i step 1 until n do begin x : = a [i, /]; for k: = i — 1 step— 1 until 1 do x : = x—a[j, k] x n[t, k]; if j =; i then begin dl: = dl x x; Пл:=0 then begin d2 : = 0; go to fail end; L}: if abs(dl)^s 1 then begin dl: == dl X 0. 0625; d2 : = d2 + 4; go to LI end; L2i if abs (dl) < 0. 0625 then begin dl: =dl X 16; d2 : = d2 — 4; go to L2 end; if x < 0 then go to fail; /?[(]:= \lsqrt (x) end else • c[/, i]:=x Xp [i] end ij end choldet 1; procedure cholsol 1 (n) data: (r, a, p, b) result: (x); value r, n; integer r, n; array a, b, p, x; comment — процедура решения системы уравнений Ах — b. 26 где А — положительно определенная симметрическая матрица и В — матрица г правых частей размера пХг. Процедуре cholsoll должна предшествовать про- цедура choldet 1 для вычисления матрицы L в виде двух массивов а (i, /'] и p [i\. Решение системы АХ = В найдено как результат последовательного решения систем LY == В и U X = Y. Матрица В сохранена для выполнения итерацион- ного уточнения. Матрица Х размера пХг размещена на месте V. При обращении к процедуре числовой массив х необходимо задать. Например, его можно отожде- ствить с массивом b, поскольку они имеют одинаковую размерность. В этом случае обращение к процедуре будет таким: cholsol I (n, r, a, p, Ь, b); begin integer i, j, k; real г; for / i == 1 step 1 until r do begin comment — решение системы Ly == b; for i s == I step 1 until n do begin г •. = b[i, j]; for k !===(' — 1 step — 1 until 1 do г : =z —a[t, k] x K[k, j]; x[i, /]i==z x p[i] end i; comment — решение системы Ux == у; for i i == n step — I until 1 do begin zs =x[i, /•]; for k i == t + 1 step 1 until n do г : = z — a[k, t] X x \k, j}; x[i, /]i=z X p[f] end (" end / end cholsol 1; procedure cholinversion 1 (я) trans: (a) exit: (fail); value n; integer n; array a; label fail; comment — элементы положительно определенной симметрической матрицы А размещены в массиве а [i, /'], i == 1, . . ., n + 1, / = 1, . . ., n (см. п. 5). Процедура реализует разложение по схеме Холецкого А == UJ, где U — матрица, транспонированная L. Нижняя треугольная матрица L накапливается в свободной части массива а (см. п. 5). На месте диагональных элементов матрицы Ь размещены величины, им обратные. В процессе обращения на место L записы- вают ей обратную матрицу, а на это место, в свою очередь, записывают элементы матрицы, обратной А. Исходная матрица сохраняется, так что возможно итера- ционное уточнение обратной матрицы. Выход из процедуры по метке fail проис- ходит так же, как в процедуре choldet /; begin integer (, ;', k, il, jl; real x, y, s; comment—формирование матрицы L; for t s = 1 step 1 until n do begin il: = ( + 1; for j ;== i step 1 until n do begin jl i = /' + 1; x; =a[f, /j; for k : == i — 1 step — 1 until 1 do x : =x—a[jl, k] X a[il, k]; if ; == i then begin if л:г$0 then go to fail; a[il, i]: = у : = l/sqrt (x) end else a[jl, i]i=xxy end ;• endi; comment — обращение матрицы L; for ( i = 1 step 1 until n do for / : = i -t- 1 step I until n do begin г : s= 0; 27 il :=?+!; for ft : = /—1 step— 1 until i do г s =г—а[/7, ft] x a ]ft + 1, l then begin dl: = d/ X 0. 0625; d2 : = d2 + 4; go to LI end; L2: if aus (dl) < 0. 0625 then begin dl: =dl x 16; d2 : = ri2 — 4; go to L2 end; if A: < 0 then go to fail; a[p}:= l/sqrf (x) 28 end; else a\p]:=xxa\r] end / end t. end choldet 2; procedure cholsol 2 (n) data: (r, a) trans : (b}; value .;, л; integer n, r; array a, b; comment — решение системы уравнений AX == В, где А — положительно опре- деленная симметрическая матрица и В—матрица размера пХг, сформирован- ная из г векторов правых частей. Процедуре cholsol 2 должна предшествовать процедура choldet 2 для вычисления и размещения матрицы L в массиве а [t]. Решение системы AX =s В найдено как результат последовательного решения систем LV = В и UX = V. Для размещения матрицы V, а затем и Х исполь- зован массив Ь; begin integer t, ;', ft, , q, •>; real x; ior /':==! step 1 until r do begin comment — решение системы LY == В;] p:=l for ( : == 1 step 1 until n do , begin q : == t — 1; x:==b[i, j}; tor ft : = 1 step 1 until q do begin x : == x — a [p\ X Ь [ft, /']; p : = p+ 1 end ft; b[i, j}:^x X a\p\; p: =p+1 end (; comment — решение системы UX === Y; for ( : == n step — 1 unti 1 1 do begin s i == p : = p — 1; 9:==i+l; x:^b[i, П; for ft : == n step — 1 until q do begin x: = x — a [s] X Ь [ft, /']; s:=s—ft+1 end ft; b[i, j}:== x X a[s] end ( end ; end dwisol 2; procedure cholinversion 2 (n) trans: (a) exit: ([ail); value n; integer n; array a; label fail; comment — элементы положительно определенной симметрической матрицы А n(re+ 1) размещены по строкам в массиве a[f], (==1, . . ., . (см. п. 5). Процедура реализует разложение по схеме Холецкого А == LU, где U — матрица, транс- понированная L. Матрица L накапливается в массиве а. На месте диагональных элементов матрицы L размещены им обратные. В процессе обращения на место L записывают матрицу ей обратную, а на это место, в свою очередь, записывают элементы матрицы, обратной А. Выход из процедуры по метке fail происходит так же, как в процедуре choldet 2; begin integer i, /', ft, p, q, r, s, t, u; real x, y; comment—формирование матрицы L; p:=0; for t : = 1 step 1 until n do begin q: = p+ 1; .- = 0; 29 for /: = 1 step 1 until f do begin x : == a [p + 1]; for k: = q step 1 until p do begin r : = r 4- 1; x : = x — a [k] X а [г] end k; r:=r+l; p : =p+ 1; if t = / then begin if ;<г$0 then go to fail', a[p}:=\isqrt(x) end else a [p] : = x x a [r] end / end t; comment — обращение матрицы L; for f : == 2 step 1 until n do begin p =p+ 1; r =r+i; у =a[r+l], for /': == 2 step 1 until г do begin p =p+ 1; s =f.=t+l; и =t—2; x =0; for k : = r step — 1 until p do begin x: = x — a [k\ x a [s], s : == s — u; и : = и— 1 end k; a [p] : = x x у; end /•; end f; comment—вычисление матрицы, обратной А; р:=0; for i : = 1 step 1 until n do begin r : ==• p + t; for /': = 1 step 1 until t do begin s : == f; <:==P:=JO+1; .<•: = 0; for k: == i step 1 until n do begin x =x+a [s] X а Ш; s = s 4- k; t =i+k end ^; a [p] : = ж end /' end i end cholinversion 2:. Третья группа процедур procedure symdet (n) trans: (a) result: (p, dl, d2) exit: (fail); value ?;; integer d2, n; array a, p; real dl; label fail; comment — элементы положительно определенной симметрической матрицы А размещены в массиве а [i, /'], i = 1, . . ., n; j =• 1, . . ., n (см. п. о). Процедура 30 реализует разложение вида А == LDU, где LI — матрица, транспонированная Ь. Здесь L — нижняя треугольная матрица с единичной диагональю; О — диаго- нальная матрица. Элементы матрицы, за исключением единичной диагонали, размещены в массиве а (см. п. 5). Величины, обратные элементам матрицы D, размещены в массиве р [i], < = 1, . . ., п. Матрица А сохранена для возможного уточнения решения. Определитель матрицы вычислен по формуле dIX2^d2. Выход из процедуры по метке fail происходят в том случае, если из-за ошибок округления элемент диагональной матрицы D оказался равным нулю; begin integer t, /', k; real x, у, г; dl:= 1; d2: = 0; for с = 1 step 1 until n do for ;': == 1 step 1 until (' do begin x: = a [/, i}\ if i = / then begin for k: = i—1 step—1 until 1 do begin y: = a [i, k], z: ==a[t, k]: == yxp[til; x: = x—ухг end k; dl:==dlxx: if x==0 then begin d2: = 0; go to fail end; L}: if abs(d!)^ 1 then begin dl; =d/x 0.0625;. d2: == d2 + 4; go to LI end; L2: if abs (dl) < 0.0625 then begin dl:=dl X 16; d2: == d2 — 4; go to L2 end; p [f]: == l/x end else begin for k: = j—1 step—1 until 1 do x: = x—a[i, k] X a[/, k]; a [(. /•]: = x end end // end symdet; procedure symsol (n) data: (r, a, p, b) result: (x); value n, r, integer n, r; array a, p. b, x; comment—решение системы уравнений AX == В, где А—положи- тельно определенная симметрическая матрица и В — матрица г правых частей размера пХг. Процедуре symsol должна предшествовать процедура symdet для вычисления и размещения матрицы L в массиве а [(', /'], а матрицы D — в массиве Р [i]. Решение системы АХ = В найдено как результат последовательного реше- ния систем LY = В и DUX == Y. Матрица В сохранена для выполнения итера- ционного уточнения решения X. Матрица Х размера nX r размещена на месте Y. При обращении к процедуре необходимо задать числовой массив x. Например, его можно отождествить с массивом Ь, поскольку они имеют одинаковую размер- ность. В этом случае обращение к процедуре будет таким: symsol (n, г, а, р, Ь, b); begin integer t, /, k; real y; for /•: = l step 1 until r do begin comment — решение системы LY = B; for t; = ! step ! until n do 31 begin y: ==b[i, /•]; for k: == i — 1 step — 1 until 1 do y: = y—a[i, k] x x[k, j]; x[i, j]: =(/ end i; comment—решение системы DUX == Y; for f: = n step—1 until 1 do begin y: =x[i, j] Xp [i]; for k: =• i -)- 1 step 1 until n do y:=y—a[k, i] X x[k, ;•]; •r[t. J}'-==V end i end / end symsol; procedure syminversion (n) trans: (a) exit: (fail); value n; integer n; array a- la- bel fail; comment—элементы положительно определенной симметрической матрицы А размещены в массиве a [i, j}, i == 1, . . ., n + 1, /' == 1, . . ., n (см. п. 5). Процедура реализует разложение вида А == LDU, где U — матрица, транспонированная L, L — нижняя треугольная матрица с единичной диаго- налью, D — диагональная матрица. Элементы матриц L и D за исключением единичной диагонали матрицы L накапливаются в свободной части массива а (см. п. 5). На месте элементов матрицы D размещены величины, им обратные. В процессе обращения на место матрицы L записывают ей обратную, а на место обратной L и матрицы D, в свою очередь, записывают элеменчы матрицы, обрат- ной А. Исходная матрица сохраняется, так что возможно итерационное уточне- ние обратной матрицы. Выход из процедуры по метке fail происходит так же, как в процедуре symdet; begin integer i, /', k, il, jl; real x, у, г; comment — формирование матрицы L; for i: =- 1 step 1 until n do begin il: == t + 1; for j: = 1 step 1 until (' do begin jl: == /Ч- 1; x: = a [j, i]; if i== j then begin for k: = j—1 step—1 until 1 do begin y: =a[il, k\; г: = a [il, k]: == y x a [k -+- 1. k]; x: == x—y X z end k; if x=0 then go to fail; a [il, i]: = l/x end else begin for k:=j—1 step—1 until 1 do x: = x —a [il, k] X a [jl, k}; a[il, t\:=-x end end j end i comment — обращение матрицы L; for i : = 2 step 1 until n do begin il: = t + 1; for ;: = 2 step 1 until (' do begin jl: = j — 1; x: = —a[il, j\]; for k: == i — 1 step — 1 until j do x: ==x—a[il. k] X a[k-[- 1, f/1; a[il. jl]:=x end / 32 end i; comment—вычисление матрицы, обратной А; for j: == j step 1 until n do for i: == / step 1 until n do begin il: == t + li x: =a(il, j}; if i^=/ then begin for k: = il step 1 until n do x:==x+a[k+l, i] Xa[k+ 1, /] end else for k: = t'l step 1 until n do begin y: =a[A+ 1, /]; z:==a[k+l, /]:=a[fe+l, k]Xy; X: = Х+уX 2 end k', a[il, j}:=x end ij end syminversion;. 5. Организация процедур и обозначения Первая группа процедур В процедуре choldet 1 матрица А определена как массив размера пХ n', однако правильная запись элементов необходима лишь в верхнем треугольнике, поскольку поддиагональные элементы в процедуре не используются. Матрица L, исключая ее диагональные элементы, накапливается на месте поддиагональных элемен- тов матрицы А. Необходимые для последующих вычислений величины —— , 41 i = 1, . . ., n сформированы отдельно в виде вектора р. Все это означает, что при выходе из процедуры массивы а и р имеют вид "и °i2 ^з "и ^21 °22 °23 °24 Ач1 Лчй ^.1-4 ^' 32 "33 "34 '42 ^43 a•^^ (34) При обращении к процедуре choldet 1 обязательно вычисляется определи- тель матрицы А. Поскольку величина определителя может принимать значения в широком диапазоне, для его представления введены два параметра dl и d2, ра- зная которые, можно вычислить определитель по формуле det A = 2 ~dl, где d2—целое число, а модуль dl удовлетворяет неравенству -ic-sS I ri/ | < 1. Процедура cholsol 1 определяет решение матричной системы уравнений вида АХ == В для г различных правых частей, формирующих матрицу В. Матрица решений Х размера пХ г накоплена отдельно от В, так что матрицу В можно исполь- зовать для вычисления невязок. В процедуре cholinversion t память распределена по-иному, нежели в про- цедуре choldet 1, для того, чтобы иметь возможность записать элементы вычислен- ной матрицы А'1 в естественной форме. Если ввести обозначения Р = L-i, Q == А-1, - ^'илкинсон (35) 33 то размещение матриц L, P, Q в массиве а в процессе работы алгоритма имеет вид "11 "12 ''IS "14 Pll "22 "23 "24 ""11 "12 "13 "U" "11 "12 "13 "14 9ll "22 "33 "24 'П1 "22 "23 "24 ^21 ^22 "33 °34 hi ^32 ^33 a44 _/41 /42 /43 /44 . 136) ?21 ?22 "33 "34 |, | 921 922 "33 "34 931 932 933 "44 .9,1 942 943 944 РЗ! /'32 РЗЗ "44 .Pll P42 Р43 ?44. Вторая группа процедур В процедуре choldet 2 поддиагональные и диагональные элементы симме- трической матрицы А описаны как одномерный массив а, состоящий из —л (га-(- 1) элементов. Для матрицы 4-го порядка это соответствие имеет вид (37) "11 Oi "21 022 "2 "3 "31 "32 "33 а, "в Яв а^ 042 "43 "44 a, as ад ащ / /, и элементы ац переходят в элементы о^, где и == -п- (г — 1) + /'. Элементы /,/ матрицы L вычисляют последовательно один за другим и размещают в одномер- ном массиве о на месте элементов aij матрицы А. На месте элементов 1ц располо- жены им обратные. При каждом обращении будет вычислен определитель матрицы, для представления которого использована формула det А = 2d2X dl. В процедуре cholsol 2 матрица Х решения уравнения АХ = В размещена на месте матрицы В. Это связано с тем, что при применении процедур этой группы нельзя вычислить невязки. В процедуре cholinversion 2 одномерный массив а, соответствующий матрице А, сначала заполняют элементами матрицы L по аналогии с процедурой choldet 2, а затем последовательно элементами матриц L"1 и А~1. Таким образом, матрица А"1 будет размещена на месте исходной матрицы А. Третья группа процедур Три процедуры symdet, symsol, syminversion во многом сходны с процедурами группы chol 1 и отличаются лишь тем, что вместо величин 1^ вычисляют d^1. Кроме того, в процедуре symdet на промежуточных стадиях вычисления исполь- (W Ч эуются величины а,/. Элементы /(•;; вычисляются по строкам. В процессе вычисле- ния четвертой отроки матрицы пятого порядка соответствующие двумерные массивы принимают вид "11 °12 "13 "14 "15 /21 "22 "IS "24 "25 /31 /32 "33 "34 "36 ' X X X "44 "45 |_Х X X X a„_ "11 "n "13 "и ^б ^1 "22 ^3 ^ "25 fv IV /31 /32 "33 ^'M "Зб ' "41 "43 "43 "44 "45 X X X X 0^ "ll "12 "13 "14 "15 /21 "22 "23 "24 "26 /31 /32 "33 "34 "35 /41 /48 /43 "44 "45 X X X X Ogg ,(38) 84 где элементы, отмеченные знаком X, могут принимать произвольное значение. Дополнительно происходит накопление вектора величин, обратных d,, так же как в процедуре choldet 1. 6. Оценка точности решения Для схемы Холецкого характерны высокая точность и устойчивость числен- ного процесса. Это следует из условия i^-0"^1- •••")• (39) /=i которое гарантирует ограниченность величин lij, если ограничены величины а;/. На самом деле, если все ац удовлетворяют неравенству | а,/1 <: 1,то все Iff удовлетворяют тому же неравенству. Именно поэтому реализация схемы Холец- кого даже на вычислительной машине с фиксированной запятой не вызывает затруднений. Далее при точных вычислениях справедливы соотношения || L ]|, = || А ||У2 = ^2, 11L-1 II, = Ц А-1 ||У2 ^ ^2, (4С) если Ki — собственные значения матрицы А, расположенные выбывающей после- довательности. Вторая схема разложения не столь удобна при вычислениях с фик- сированной запятой, поскольку элементы матрицы L могут оказаться существенно большими единицы, несмотря на выполнение условия | a:j I •< 1. Что же касается вычислений с плавающей запятой, то верхние границы ошибок для обеих схем в основном совпадают. Границы ошибок вычисления определителя, отыскания решения системы уравнений, обращения матрицы всегда могут быть выражены с помощью числа обусловленности я: h. (41) Если вычисления выполняются с /-разрядными числами в арифметике с пла- вающей запятой, то вычисленная матрица удовлетворяет равенству LL1 = А + F, (42) причем для матрицы F справедливы следующие оценки: || F Из ^ k^l-1 Ц А На; II F ||., ss k^l-1 max | а ц I. (43) Здесь и всюду в дальнейшем величины ki порядка единицы и зависят от спо- соба округления. При статистическом подходе к анализу ошибок || F |] g вряд ли будет превышать величины га3^"' ИАНг илип^^тах | а.ц \ при условии, что га значительно. Если собственные значения матрицы А + F обозначить через К; то из первой оценки (43) следует - k^2-t\ ^ \, - К', ^ k^lh-\. Поскольку det (LL1') = П^с == Ш'{, то П (К, - k^l^-t ^) ^ defdl^) ^ П (К( + ^22-' Xi), откуда (44) (45) П (, kn3^-^}-^^^1^ О+^З-' h}' 1Ц1— *i" i ~^)'^ detA ^IIV 1 1 A,/ (46) Таким образом, относительная ошибка при вычислении определителя зависит • , ,.3/2n—f„, яовном от величины fei/i ^ "•• 35 г , ,.3/2п—1.. основном от величины fei/i ^ "•• Из неравенства (44) следует, что при х > З'/*,/!^2 (47)' :обственное значение ?-„ может стать отрицательным, т. е. ошибки округления могут привести к потере свойства положительной определенности, если вели- чина У. достаточно велика. Существует три источника ошибок при вычислении обратной матрицы: процесс разложения, обращение матрицы L и вычисление матрицы (L^^L-1. Ошибки последних двух составляющих зависят от величины к172, так что, если матрица А плохо обусловлена, результирующая ошибка в основном зависит от ошибок разложения. Если матрица Q— результат обращения, то справедливо черавенство "M"^"2 ^з»372^. (48) правая часть которого на практике редко превосходит величину п3^"^. При решении системы уравнений Ax = b преобладающими оказываются ошибки решения систем уравнений с треугольными матрицами. Результирующая ошибка удовлетворяет неравенству L^^!^ 3/2,-f -у.—s^k^n i к, ||А-Ч)|| Следует отметить, что во всех процедурах точность решения в основном зависит от точности вычисления скалярного произведения. Если имеется воз- можность накапливать скалярное произведение с двойной точностью, то границы сшибок можно существенно уменьшить. Так, например, справедливо условие LL1' = А + F, И F ||, sg й„п1/2 2-' [I A Ib. (50) причем и эта оценка завышена. Вывод формул для границ этих ошибок дан в рабо- тах [94, 95]. Модифицированные процедуры choldet I, cholsol 1 и cholinversion 1, в которых реализовано накопление скалярного произведения, приведены в алго- ритме 1.2. Причина выхода из процедуры зависит от текущей стадии вычислений. Если в момент останова параметр i <^ п и величина (—1\,) существенно более 2~' || А |]а, то полученное разложение просто не имеет смысла. Это связано с тем, что матрица А, которая должна быть положительно определенной, или неточно вычислена, или неверно задана из-за ограниченной точности вычислителя. В любом случае схему Холецкого нельзя завершить без дополнительных моди- фикаций. Однако разложение вида LDL7 можно было бы закончить даже при отрицательных значениях d,, хотя устойчивость процедуры гарантируется только для положительно определенных исходных матриц. Если процедуру choldeti иди symdet используют для определения методом обратной итерации собственного вектора матрицы В, то матрицу А формируют равной В—U, где 'k выбирают близким к наименьшему собственному значению матрицы В, но меньше его. В этом случае ошибки округления как раз и могут нарушить положительность матрицы В—М, и останов чаще всего происходит при i = п. Тогда допустимо положить /,,„ = 2"'!] А И^72. Что касается вычисления определителя, то обычно происходит останов при определении элемента <„„, а поэтому множитель Гц вычисляют до проверки диагонального элемента на поло- жительность. 7. Примеры использования процедур и результаты их проверки Разложение LDL имеет незначительные преимущества по сравнению со схе- мой Холецкого. Разложение LL1 в сочетании с L R-алгоритмом Рутисхаузера [70] находит применение для вычисления собственных значений и векторов сим- 36 (49) метрических матриц, хотя обычно этот алгоритм применяют лишь для ленточных матриц. Разложение Холецкого используют также для определения корней урав- нений: det (А — ?iB) == 0; det (АВ — U) == 0, (51) где В — положительно определенная и А — симметрическая матрицы. Если В == Ll^, то эти уравнения принимают вид det [L-1 А (L-1 У1 — XI] = 0, det [l^ AL — U] == 0, (52) где матрицы L^A (Ь"^ и I^AL симметрические. Приведенные выше процедуры были применены к матрице 7-го порядка, являющейся главным минором неограниченной гильбертовой матрицы. Для того чтобы избежать ошибок округления, матрица была масштабирована с помо- щью множителя 360 360 так, что все элементы оказались целыми числами. Таким образом, матрица А имела вид ~ 360 360 180 180 120120 180 180 120120 90090 120120 90090 72072 90090 72072 60060 72072 60060 51480 60060 51 480 45 045 51 480 45 045 40040 72072 60060 51 480- 60060 51 480 45045 51480 45045 40040 45045 40040 36036 40040 36036 32760 36036 32760 30030 32760 30030 27 720 "ЗбОЗбО 180180 120120 90090 180180 120120 90090 72072 60060 51480 45045 120120 90090 72072 60060 51480 45045 40040 А= 90090 72072 60060 51480 45045 40040 36036 (53). 45045 40040 36036 Выбор матрицы объясняется тем, что, во-первых, она плохо обусловлена ( к ==-„-• 10е ) и, во-вторых, представляет интерес сравнить полученные здесь результаты обращения с результатами применения алгоритма 1.2, в котором реализовано накопление скалярного произведения. Обращение матриц выполнялось с помощью процедур cholinversion 1, 2 и syminversion, а решение системы АХ = I было получено с использованием про- цедур cholsol I, 2 и symsol. Таким образом, было исследовано 6 способов обраще- ния матрицы. Вычисления проводились на вычислительной машине KDF9 с двоичной системой счисления, плавающей запятой и длиной машинного слова, равной 39 двоичным разрядам. Принимая во внимание большое значение обусло- вленности матрицы А, можно сказать, что все 6 результатов были достаточно точными. Поскольку полученные результаты существенно не отличались по точности один от другого, то ниже приводятся лишь два результата обращений, получен- ных с помощью процедур cholsol I (табл. 1) и cholinversion 1 (табл. 2). Обе про- Таблица I Результат вычисления матрицы А-' с помощью совокупности процедур choldet 1 и cholsol I I столбец II столбец III столбец +1,35971314283 —3,26325113738 +2,44740603601 —8,15793917875 +1,34604954738 —1.07683295031 +3,33303738252 10-10-10-10-10-10-10" —3,26325113753.10-3' +1,04422850066-Ю-1 —8,81060620282.10-' +3,13264074692.10° —5,38420043662.10» +4,43041069139.Ю» —1,39988227242.Ю» +2,44740603295. Ю-2 —8,81060620311.10-' +7,92950934930.10» —2,93684512786.10' +5,1919084898l.l0' —4,36119354458.101 +1,39988675007.10' 37 Продолжение табл. I IV столбец V столбец VI столбец —8,16793917156 +3,13264074649 —2,93684512777 4-1,11879649588 —2,01907567906 +1,72294292795 —5,59956005459 10-' 10" 10» 10' 10' 102 10* +1,34604954657 —5,38420043660 +5,19190848985 —2,01907567907 +3,70163885040 —3,19821597914 +1,04991937763 10-» 10» 101 10' 10' 102 102 —1.07683294804 +4,43041069058 —4,36119354451 +1,72294292794 —3.19821597912 +2,79117246791 —9,23930348819 10-' 109 10*ю2102 102 10» VII столбец +3,33303738787 —1,39988227277 +1,39988675014 —5,59956005457 +1,04991937762 —9,23930348819 +3,07977132849 10-2 10» 101 10« 10' 101 101 Таблица f Результат вычисления матрицы А"' с помощью процедуры cboUncersion J I строка II строка III строка +1,36971314362.10-< —3,26325113762.10-' + 1,04422850063-Ю-1 +2,44740603601.10-' —8,81060620290-Ю-» +7,92950934908.10° IV строка V строка VI строка —8,15793917877-Ю-2 +3,13264074697.10» —2,93684512785 •101 +1,11879649600.102 + 1,34604954738-Ю-1 —5,38420043665.10" +5,19190848982.101 —2,01907567907. 102 —3,70163885041 .108 —1,07683295031 +4,43041069138 —4,36119354457 +1,72294292795 —3,19821597913 +2,79117246790 10-» 10» 10» 10' 102 10' VII строка +3,33303738252 —1,39988227242 +1,39988675007 —5,59956005459 +1,04991937763 —9,23930348819 +3,07977132849 10-•• Ю» 10> 101 10' 101 101 цеДУры дают одинаковый результат для матрицы L, но получение обратных матриц существенно различается в этих процедурах. Процедура cholsol 1 осуществляет обращение матрицы путем последовательного выполнения прямой и обратной под- становок *. В процедуре cholinversion 1 вычисляется L""1, а А"1 определяется как произведение матриц. Приведенные результаты иллюстрируют сказанное выше в п. 6. Элементы патриц при двух методах расчета точно совпадают в девяти значащих цифрах и хорошо согласуются до 11 десятичных разрядов мантиссы, несмотря на то, что оба результата, вообще говоря, верны лишь до 5 значащих цифр. Это обусло- влено тем, что основной источник ошибок — процесс приведения исходной матрицы к треугольной форме, но данная процедура одинакова при обоих спо- собах расчета. По той же причине результат обращения с помощью процедуры cholsol I дает практически симметрическую матрицу, ведь ошибки от прямой и обратной подстановок несущественны, и потому вычисленное обращение очень близко к матрице (L'1)7 L-1, но с ошибками, обусловленными определением L, причем матрица (L'^L'1 симметрическая, хотя L и вычислена с ошибками. Нельзя утверждать, что приведенные выше ошибки обращения порядка 10"* обусловлены ошибками от неточного задания элементов матрицы в /-разрядной арифметике (в данном случае t = 39). Вычисленное с помощью различных процедур значение определителя матрицы А равно: 8,47418829503 Х Ю-2 Х 252 для процедуры choldet 1; 8,47371300872 х 10~2 X 252 для процедуры choldet 2; 8,47317605976 х Ю-2 X 262 для процедуры symdef. Точное значение определителя равно 381614277072600 == 8,47353913 X X 10*2 X 252. Таким образом, точность его вычисления, имея в виду значение числа обусловленности матрицы А, достаточно высокая. Второй просчет той же самой матрицы был выполнен на вычислительной машине TR4. Машинное слово в этом вычислителе представляется 38-разрядной двоичной мантиссой, однако средняя точность результатов была приблизительно на один двоичный порядок выше, чем для вычислителя KDF9. * Прямой и обратной подстановками называют метод решения систем уравнений соответственно с нижней и верхней треугольными матрицами {Прим. ред пер.). Были вычислены два наименьших собственных значения. Они очень плохо обусловлены и, хотя удовлетворяют некоторой близкой возмущенной матрице, не очень точны. Тем не менее каждый собственный вектор был определен за не- полную первую итерацию и таким образом с максимально возможной точностью (табл. 2). Таблица 2 Собственные векторы для двух наименьших собственных значений матрицы Н Собственный вектор для ?,i ,2,77553987291-Ю-2 Собственный вектор для А,, =4,41194423277-Ю-1 +4,27925794055 —1,39361241100 10-' 10- 1 -2,27225077213 +5,52600937792 10-10- 0 +3,02780088463 —5,08333303903 +6,89752171663 —7,67806578297 10-10-10-10- —5,82538031647 —1,15841897038 +1 57134904714 —3,16813976535 10-10-10-10- +7,01463733572 —5,20726524397 +3,07670997672 —1,39678066480 +4,58752082989 —9,72244601275 +1,00000000000 10-10-10-10-10-10-10- 8 2 +4,01474709225 .-3,67167116269 +2,49511499438 —1,24479436263 +4,34794099107 —9,55880557676 +1,00000000000 10-10-10-10-10-10-10° Процедура cxinvit. В качестве теста для процедуры cxinvit была выбрана комплексная матрица А + iB, где 1 3 2 1 4 3 2 3 1 ~1 Г 3 1 2 124 4 2 1 1 3 5 I I 3 5 1 1315 Были вычислены два наименьших собственных значения и соответствующие собственные векторы (табл. 3). Таблица 3 Собственные векторы, вычисленные процедурой cxinvit Собственный вектор для К, = 2,22168234771- 10»+t • 1,84899335974- 10» Re Im —7,96620817435-Ю-' —1,78834394738-10-1 —2,52814311194-10-' +1,00000000000-10° +3,04980785941•10-» +4,29724147742-10-1 +3,81734049494-10-2 +0,00000000000 Собственный вектор для >., = 1,36566930995.10»-t • 1,40105359747.10° —8,46598714563-10-* —8,79417015790-10-2 +1,00000000000-10° —4,32089374804-10-' +7,30203551264-10-' —3,87901799614.10-» +0,00000000000 —4,33426369343.10-» 384 Список литературы 1. Bartels R. H. A Numerical Investigation of the Simplex Method. Tech. Report. Computer Science Department, Stanford University, California, 1968. 2. Bartels R. H., Golub G. H. The Simplex Method of Linear Programming Using LU-decomposition. Comm. ACM. 12, pp. 266—268, 1969. 3. Barth W., Martin R. S., Wilkinson J. H. Calculation of the Eigenvectors of a Symmetric Tridiagonal Matrix by the Method of Bisection. Numer. Math. 9, pp. 386—393, 1967. 4. Bauer F. L. Das Verfahren der Treppeniteration und verwandte Verfah renzur Losung algebraischer Eigenwertprobleme. Z. Angew, Math. Phys. 8, pp. 214— 235, 1957. 5. Bauer F. L. Optimally Scaled Matrices. Numer. Math, 5, pp. 73—87, 1963. 6. Bauer P. L. Elimination with Weighted Row Combinations for Solving Linear Equations and Least Squares Problems. Numer. 7, pp. 338—352, 1965. 7. Bauer F. L., Heinhold J., Samelson K., Sauer R. Moderne Rechenanla- gen. Stuttgart, Teubner, 1965. 8. Bauer F. L. Geriauigkeitsfragen bei der Losung Linearer Gleichungssysteme. Z. Angew. Math. Mech. 46, pp. 409—421, 1966. 9. Bauer F. L., QD-method with Newton Shift. Tech. Rep. 56 Computer Science Department, Stanford University, 1967. 10. Bowdler H. J., Martin R. S., Peters G., Wilkinson J. H. Solution of Real and Complex Systems of Linear Equations. Numer. Math. 8, pp. 217—234, 1966. 11. Bowdler H., Martin R. S., Reinsch C., Wilkinson J. H. The Q/?-and QL- algorithms for Symmetric Matrices. Numer. Math. 11, pp. 293—306, 1968. 12. Businger P., Golub G. Linear Least Squares solutions by Householder Transformations. Numer. Math. 7, pp. 269—276, 1965. 13. Dantzing G. В. Linear Programming and Extensions. Princeton, Princeton University Press, 1963. 14. Dekker Т. 1., Hoffman W. Algol-60 Procedures in Linear Algebra. Mathe- matical Centre Tracts, 23. Mathematisch Centrum, Amsterdam 1968. 15. Dubrulle A., Martin R. S., Wilkinson J. H. The Implicit QD-algorithm. Numer. Math. 12, pp. 377—383, 1968. 16. Eberlein P. 1. A Jacobi-like Method for the Automatic Computation of Eigenvalues and Eigenvectors of an Arbitrary Matrix. J. Soc. Indust. Appl. Math. 10, pp. 74—88, 1962. 385 17. Eberlcin P. 1., Boothroyd I. Solution to the Eigenproblem by a Norm- reducing Jacobi-type Method. Numer. Math. 11, pp. 1—12, 1968. 18. Eberlein P. 1. Solution to the Complex Eingenproblem by a Norm-redu- cing Jacobi-type Method. Numer. Math. 14, pp. 232—245, 1970. 19. Engeli M., Ginsburg Th., Rutishauser H., Stiefel E. Refined Iterative Methods for Computation of the Solution and the Eigenvalues of Self-adjoint Boun- dary Value Problems. Mitt. Inst. angew. Math. E. I. H. Zurich, No. 8, Basel, Birk- hau'ser, 1959. 20. Forsythe G. E. Crout with Pivoting. Comm. ACM 3, pp.507—508, 1960. 21. Forsythe G. E., Henrici P. The Cyclic Jacob! Method for Computing the Principal Values of a Complex Matrix. Proc. Amer. Math. Soc. 94, pp. 1—23, 1960. 22. Forsythe G. E., Golub G. On the Stationary Values of a Second-degree Polynomial on the Unit Sphere. J. Soc. Indust. Appl. Math. 13, pp. 1050—1068, 1965. 23. Forsythe G. E.,'Moler С. В. Computer Solution of Linear Algebraic Systems. Englewood Cliffs, New' Jersey, Prentice Hall, 1967 (имеется русский перевод; Форсайт Дж., Молер К, Численное решение систем линейных алгебраических уравнений. M., «Мир», 1969). 24. Forsythe G. H., Straus E. G. On best Conditioned matrices. Proc. Amer. Math. Soc. 6, pp. 340—345, 1955. 25. Fox L. Practical Solution of Linear Equations and Inversion of Matrices. Nat. Bur. Standards Appl. Math., Ser. 39, pp. 1—54* 1954. 26. Francis J. G. F. The Q^-transformation — a Unitary Analoque to the L^-transformation, Parts I, II. Comput. J. 4, pp. 265—271, 1961, pp. 332—345, 1962. 27. Cinsburg T. The Conjugate Gradient Method. Numer. Math. 5, pp. 191— 200, 1963. 28. Givens 1. W. A Method for Computing Eigenvalues and Eigenvectors Sugges- ted by Classical Results on Symmetric Matrices. Nat. Bur. Standards Appl. Math. Ser. 29, pp. 117—122, 1953. 29. Givens I. W. Numerical Computation of the Characteristic Values of a Real Symmetric Matrix. Oak Ridge National Laboratory, ORNL-1574, 1954. 30. Golub С., Kahan W. Calculating the Singular Values and Pseudo-inverse of a Matrix. J. SIAM. Numer. Anal., Ser. B2, pp. 205—224, 1965. 31. Colub G., Kahan W. Least Squares, Singular Values, and Matrix Appro- ximations. Aplikace Matematiky 13, pp. 44—51, 1968. 32. Gregory R. T. Computing Eigenvalues and Eigenvectors of a Symmetric Matrix on the ILLIAC. Math. Tab. and other Aids to Сотр., 7, pp. 215—220, 1953. 33. Golub G. H. Numerical Methods for Solving Linear Least Squares Problems. Numer. Math. 7, pp. 206—216, 1965. 34. Golub G. H., Reinsch C. Singular Values Decomposition and Least Squares Solutions. Numer. Math., 14, pp. 403—420, 1970. 35. Henrici P. On the Speed of Convergence of Cyclic and QuasicycHc Jacob. Methods for Computing Eigenvalues of Hermitian Matrices. J. Soc. Indust. Appli Math. 6, pp. 144—162, 1958. 36. Hestenes M. R., Stiefel E. Methods of Conjugate Gradients for Solving Linear Systems. Nat. Bur. Standards J. of Res. 49, pp. 409—436, 1952. 37. Hestenes M. R. Inversion of matrices by Biorthogonalization and Related Results. J. Soc. Indust. Appl. Math. 6, pp. 51—90, 1958. 38. Householder A. S. Unitary Triangularization oi a nonsymmetric Matrix. J. Assoc. Comput. Mach. 5, pp. 339—342, 1958. 39. Householder A. S. Principles of Numerical Analysis, N. Y., 1953 (есть русский перевод: Хаусхолдер А. С. Основы численного анализа, M,, изд-во »Ин. лит.», 1956). 386 40. Householder A. S., Bauer F. L. On Certain Methods for Expanding the Characteristic Polynomial. Numer. Math. 1, pp. 29—37, 1959. 41. Jacob! С. С. J. Uber ein leichtes Verfahren. die in der Theorie der Saku- larstorungen Vorkommenden Gleichungen numerisch aufzulosen. Crelle's Journal 30, pp. 51—94, 1846. 42. Jennings A. A direct Method for Computing Latent Roots and Vectors of a Symmetric Matrix. Proc. Cambridge Philos. Soc. 63, pp. 755—765, 1967. 43. Kahan W. Accurate Eigenvalues of a Symmetric Tridiagonal matrix. Tech. Rep. No. Cs.— 41, Computer Science Department, Stanford University, 1966. 44. Kogbetliantz E. G. Solution of Linear Equations by Diagonalization of Coefficients Matrix. Quart. Appl. Math., 13, pp. 123—132, 1955. 45. Кублановская В. H. О некоторых алгорифмах для решения полной проблемы собственных значений.— «Журнал вычислительной математики и ма- тематической физики», т. 1, № 4, стр. 555—570, 1961. 46. Lauchli P., Waldburger H. Berechnung von Beulwerten mit Hilfe von Mehrstell^noperatoren, Z. Angew. Math. Phys. 11, pp. 445—454, 1960. 47. Martin R. S., Wilkinson J. H. Symmetric Decomposition of Positive Definite Band Matrices. Numer Math., 7, pp. 355—361, 1965. 48. Martin R. S., Peters G., Wilkinson J. H. Symmetric Decomposition of a positive definite Matrix. Numer. Math., 7, pp. 362—383, 1965. 49. Martin R. S., Peters G., Wilkinson J. H. Iterative Refinement of the Solution of a Positive Definite System of Equations. Numer. Math. 8, pp. 203— 216, 1966. 50. Martin R. S., Wilkinson J. H. Solution of Symmetric and Unsymmetrio Band Equations and the Calculation of Eigenvectors of Band Matrices. Numer. Math. 9, pp. 279—301, 1967. 51. Martin R. S., Wilkinson J. H. Reduction of the Symmetric Eigenproblera Ax-^Bx and Related Problem to Standard Form. Numer. Math. II, pp. 99—110, 1968. 52. Martin R. S., Reinsch C., Wilkinson J. H. Householder's Tridiagonali- zation of a Symmetric Matrix. Numer. Math. II, pp. 181—195, 1968. 53. Martin R. S., Wilkinson J. H. Similarity Reduction of a General Matrix to Hessenberg Form. Numer. Math. 12, pp. 349—368, 1968. 54. Martin R. S., Wilkinson J. H. The Modified L/?-algorithm for Complex Hessenberg Matrices. Numer. Math. 12, pp. 369—376, 1968. 55. Martin R. S., Wilkinsin J. H. The Implicit QL-algorithm. Numer. Math. 12, pp. 377—383, 1968. 56. Martin R. S., Peters G., Wilkinson J. H. The Q^-algorithm for Real Hessenberg Matrices. Numer. Math. 14, pp. 219—231, 1970. 57. Martin R. S., Peter G., Wilkinson J. H. The Q^-algorithm for Band Sym- metric Matrices. Numer. Math. 16, pp. 85—92, 1970. 58. McKeeman W. M. Crout with Equilibration and Iteration. Comm. ACM 5, pp. 553—555, 1962. 59. Ortega J. M. An Error Analysis of Householder's Method for the Symmetric Eigenvalue Problem. Numer. Math. 5, pp. 211—225, 1963. 60. Ortega J. M., Kaiser H. F. The LZA — and Q^-methods for Symmetric Tridiagonal Matrices. Comput. J. 6, pp. 99—101, 1963. 61. Osborne E. E. On Pre-conditioning of Matrices. J. Assoc. Comput. Math. 7, pp. 338—345, 1960. 62. Parlett В. N. The Development and Use of Methods of LR-tyw. STAM. Rev., 6, pp. 275—295, 1964. 63. Parlett В. N. Global Convergence of the Basis Q^-algorithm on Hessen- berg Matrices. Math. Сотр. 22, pp. 803—817, 1968. 64. Parlett В. N.. Reinsch C. Balancing a Matrix for Calculation of Eigen- values and Eigenvectors. Numer. Math, 13, pp. 293—304, 1969. 387 63. Peters G., Wilkinson J. H. Eigenvectors of Real and Complex Matrices by LR and 0^-Triangularizations. Numer. Math. 16, pp. 181—204, 1970. 66. Petrie G. W. Ill Matrix Inversion and Solution of Simultaneous Linear Algebraic Equations with the IBM-604 Electronic Calculating Punch. Simulta- neous linear Equations and Eigenvalues. Nat. Bur. Standards Appl. Math. Ser. 29, pp. 107—112, 1953. 67. Pope D. A., Tompkins C. Maximizing Functions of Rotations-Experiments Concerning Speed of Diagonal isation of Symmetric Matrices Using Jacobi's Method. J. Assoc. Comput. Mach. 4, pp. 459—466, 1957. 68. Reinsch С., Bauer F. L. Rational QR -transformation with Newton Shift for Symmetric Tridiagonal Matrices. Numer. Math. II, pp. 264—-272, 1968. 69. Rutishauser H. Automatische Rechenplanfertigung fur programmgesteuerte Rechenanlagen. Z. Angew. Math. Phys. 3, pp. 312—316, 1952. 70. Rutishauser H. Solution of Eigenvalue Problems with the L^-transfor- mation. Nat. Bur. Standards Apll. Math. Ser. 49, pp. 47—81, 1958. 71. Rutishauser H. Stabile Sonderfalle des Quotienten-Differenzen-Algo- rithmus. Numer. Math. 5, pp. 95—112, 1963. 72. Rulishauser H., Schwarz H. R. The LR-transformation Method for Sym- metric Matrices. Numer. Math. 5, pp. 273—289, 1963. 73. Rutishauser H. On Jacobi Rotation Patterns. Proceedings of Symposia in Applied Mathematics, vol. 15, pp. 219—239 Experimental Arithmetic, High Speed Computing and Mathematics, 1963. 74. Rutishauser H. The Jacobi Method for Real Symmetric Matrices. Numer. Math. 9, pp. 1—10, 1966. 75. Rutishauser H. Lineares Gleichungssystem mit Symmetrischer Positive Definitpr Bandmatrix nach Cholesky. Computing I, pp. 77—-78, 1966. 76. Rutishauser H. Description of ALGOL—60. Handbook of Automatic Computation, vol. la. Berlin-Heidelberg-New York, Springer, 1967. 77. Rutishauser H. Computational aspects of F. L. Bauer's Simultaneous Iteration Method. Numer. Math. 13, pp. 4—13, 1969. 78. Rutishauser H. Simultaneous Iteration Method for Symmetric Matrices. Numer. Math. 16, pp. 205—223, 1970. 79. Schoenhage A. Zur Konvergenz des Jacob i-Verfahrens. Numer. Math. 3, pp. 374—380, 1961. 80. Schwarz H. R. Reduction of a Symmetric Bandmatrix to Triple Diagonal Form. Comm. ACM. 6, pp. 315—316, 1963. 81. Schwarz H. R. Die Reduktion einer symmetrischen Bandmatrix auf tri- diagonale Form. Z. Angew. Mech. (Sonderheft) 45, T-76, T-77, 1965. 82. Schwarz H. R. Tridiagonalization of a Symmetric Band Matrix. Numer. Math. 12, pp. 231—241, 1968. 83. Schwarz H. R. Matrizen-Numeric. Stuttgart, Teubner, 1968. 84. Stiefel E. Einige Methoden der Relaxations rechnung. Z. Angew. Math. Phys. 3, 1—33, 1952. 85. Stiefel E. Relaxationsmethoden bester Strategie zur Losung linearer Glei- chungssysteme. Comment. Math. Helv, 29, pp. 157—179, 1955. 86. Stiefel E. Kernel Polynomials in Linear Algebra and their Numerical Applications, Nat. Bur. Standards. Appl. Math. Ser. 49, pp. 1—22, 1958. 87. Stuart G. W. Accelerating the Orthogonal Iteration for the Eigenvalues of a Hermitian Matrix. Numer. Math. 13, pp. 362—376, 1969. 88. Varah J. Ph. D. Thesis. Stanford University, 1967. 89. Wilkinson J. H. Error Analysis of Floating-point Computation. Numer. Math., 2, pp. 319—340, 1960. 90. Wilkinson J. H. Note on the Quadratic Convergence of the Cyclic Jacobi Process. Numer. Math. 4, pp. 296—300, 1962. 388 91. Wilkinson J. H. Householder's Method for Symmetric Matrices. Numer. Math. 4, pp. 354—361, 1962. 92. Wilkinson J. H. Calculation of the Eigenvectors of a Symmetric Tridiagona! Matrix by the Method of Bisection. Numer. Math. 4, pp. 362—367, 1962. 93. Wilkinson J. H. Calculation of the Eigenvectors of a Symmetric Tridiagonal Matrix Inverse Iteration. Numer. Math. 4, pp. 368—376, 1962. 94. Wilkinson J. H. Rounding Error in Algebraic Processes. London, Her Majesty's Stationery Office; Englewood Cliffs, N. J„ Prentice-Hall, 1963. 95. Wilkinson J. H. The Algebraic Eigenvalue Problem. London, Oxford University Press, 1965. (Имеется русский перевод. Уилкинсон Дж. X. Алгебраи- ческая проблема собственных значений. М., «Наука», 1970). 96. Wilkinson J. H. Error Analysis of Trasformations Based on the Use of Matrices of the Form I—2 ww". Error in Digital Computation, vol. II, L. B. Rail. ed., pp. 77—101. New York, John Wiley Sons, Inc., 1965. 97. Wilkinson J. H. Convergence of the LR; QR- and Related Algorithms. Comput. J. 8, pp. 77—84, 1965. 98. Wilkinson J. H. The QR-Algorithm for Real Symmetric Matrices with Multiple Eigenvalues. Comput. J. 8, pp. 85—87, 1965. 99. Wilkinson J. H. Global Convergence of Tridiagonal QR -algorithm with Origin Shifts. Lin. Alg. and its Appl. I, pp. 409—420, 1968. Литература, добавленная при переводе книги 100. Автоматическое программирование, численные методы в функциональ- ный анализ. Сб. работ под ред. Фадеевой В. H. Л., «Наука», 1968, 263 с. 101. Агеев М. И. Основы алгоритмического языка АЛГОЛ-60. Изд. ВЦ АН СССР, 1965, 123 с. 102. Агеев М. И., Алик В. П., Гали с Р. М. Алгоритмы (1—50). М., ВЦ АН СССР, 1966, 106 с. 103. Агеев М. И., Алик В. П., Марков Ю. И. Алгоритмы (51—100). М., ВЦ АН СССР, 1966, 107 с. 104. Агеев М. И., Кривонос Л. С., Марков Ю. И. Алгоритмы (101—150). М., ВЦ АН СССР, 1967, 116 с. 105. Агеев и др. Алгоритмы (151—200). М., Инст. проблем упр.— ВЦ АН СССР, 1970, 218 с. 106. Агеев М. И., Марков Ю. И., Швакова Г. М. Алгоритмы (201—250). М., Инст. проблем упр.— ВЦ АН СССР, 1971, 210 с. 107. Айнберг В. Д. Основы программирования на алгоритмическом языке АЛГОЛ-60. М., «Машиностроением, 1973, 150 с. 108. Балуев А. П. и др. Сборник упражнений по АЛГОЛ-60. Изд. Ленинград- ского университета, 1967, 86 с. 109. Барсов А. С. Линейное программирование в технико-экономических задачах. М., «Наука», 1964, 278 с. 110. Белявский E. И., Деген А. Б., Этин Ю. Б. АЛГОЛ-60. Л„ Гидрометео- издат, 1966, 208 с. 111. Березин И. С., Жидков H. П. Методы вычислений. М., «Наука», 1966, 632 в. 112. Библиотека алгоритмов 16—506. М., «Сов. радио», 1975, 17& в, 113. Брудно А. Л. АЛГОЛ. М., «Наука», 1971, 80 е. , ..', -:. 114. Брудно А. Л. Программирование в содержательных обозначениях. М., «Наука», 1968, 144 с. - , 115. Бухтияров А. М., Фролов Г. Д. Сборник задач-по программированию на алгоритмических языках. М., «Наука», 1974, 240 с. .- - .г,, 389 116. Воеводин В. В. Линейная алгебра. М., «Наука», 1974, 336 с. 117. Воеводин В. В. Численные методы алгебры. Теория и алгорифмы. М., «Наука», 1966, 248 с. 118. Вычислительные методы и программирование. Под ред. Воеводина В. В. и Климова Г. П. М., изд. Московского университета, 1972, 205 с. 119. Гантмахер Ф. Р. Теория матриц. М., «Наука», 1967, 575 с. 120. Глушков В. М. Введение в АСУ. Киев, «Техника», 1972, 310 с. 121. Демидович Б. П., Марон И. А. Основы вычислительной математики. М.. «Наука», 1970, 664 с. 122. Лавров С. С. Универсальный язык программирования. М., «Наука», 1972, 183 с. 123. Лавров С. С. Введение в программирование. М., «Наука», 1973, 352 о. 124. Линник Ю. В. Метод наименьших квадратов и основы математико- статистической теории обработки наблюдений. М., «Физматгиз», 1962, 349 о. 125. Ляшенко В. Ф. Программирование для ЦВМ с системой команд типа М-20. М., «Советское радио», 1974, 414 с. 126. Сборник научных программ на ФОРТРАНе, вып. 1. Статистика, 316 с., вып. 2. Матричная и линейная алгебра, 224 с., пер. с англ. М., «Статистика», 1974. 127. Солодовников А. С. Введение в линейную алгебру и линейное програм- мирование. М., «Просвещение», 1966, 183 с. 128. Фаддеев Д. Кч Фаддеева В. Н. Вычислительные методы линейной алгебры. М.—Л., Физматгиз, 1963, 734 с. 129. Фаддеева В. Н. Вычислительные методы линейной алгебры. М.—Л., Гостехиздат, 1950, 240 с. 130. Юдин Д. Б., Гольдштейн Е. Г. Линейное программирование. М., «Наука», 1969, 424 G. Уилкинсон, Райнш СПРАВОЧНИК АЛГОРИТМОВ НА ЯЗЫКЕ АЛГОЛ Линейная алгебра Редактор издательства инж. Л. П. Строганов Технический редактор А. И. Захарова Корректоры И. М. Борейша и Н. И. Шару чина Переплет художника А. Я. Михайлова Сдано в набор 8/1 1976 г. Подписано к печати 13/VIII 1976 г. Формат бОхеО1/!,!. Бумага типографская № 1. Усл. печ. л. 24,5. Уч.-изд. л. 28,25 Тираж 42 000. Заказ 844. Цена 1 р. 74 к. Издательство «Машиностроение», 107885, Москва, Б-78, 1-й Басманный пер.,- 3 Ленинградская типография № 6 Союзполиграфпрома при Государственном комитете Совета Министров СССР по делам издательств, полиграфии и книжной торговли 193144, Ленинград, С-144, ул. Моисеенко, 10 НОВЫЕ КНИГИ 1. ДРОЗДОВ Е. А., КОМАРНИЦКИЙ В. А., ПЯТИБРА- ТОВ А. П. Электронные вычислительные машины единой системы. Объем 45 л. В книге изложены принципы построения и функционирования электронных вычислительных машин единой системы, а также основы выполнения операций на них. Материал базируется на схемах, организации и возможностях использования машин моделей ЕС-1050, ЕС-1030, ЕС-1020. Рассмотрены не только общие схемы машин, но и их элементная база, кодирование ин- формации и машинные алгоритмы выполнения арифметических, логических и других операций на различных моделях машин; описаны особенности построения, основные функциональные и структурные схемы процессоров, запоминающих устройств, ка- налов, устройств ввода-вывода, средств телеобработки; показаны аппаратные и программные средства контроля и диагностики, структура и возможности математического обеспечения, способы организации совместной работы нескольких машин. 2. ЗОЗУЛЕВИЧ Д. М. Программные средства машинной гра- фики в автоматизированном проектировании. Объем 15 л. В книге изложены основы теории алгоритмизации процессов отображения графической информации в системах автоматизиро- ванного проектирования; описаны методы построения математи- ческих моделей изделий, конструкторских документов ЕСКД и ЕСТД, а также процессов автоматического отображения изде- лий в графических конструкторских документах; рассмотрены особенности алгоритмизации и программирования задач отобра- жения графической информации, основанные на системном ана- лизе объектов и процессов; приведен ряд практических примеров применительно к конструированию и технической подготовки производства. Эти книги издательство «Машиностроение» выпустит в 1976 г.