Методы решения симметричной проблемы собственных значений и проблемы определения сингулярного разложения с оцениваемой точностью тема автореферата и диссертации по математике, 01.01.07 ВАК РФ

На правах рукописи

Мацех Анна Михайловна

МЕТОДЫ РЕШЕНИЯ СИММЕТРИЧНОЙ ПРОБЛЕМЫ

СОБСТВЕННЫХ ЗНАЧЕНИЙ И ПРОБЛЕМЫ ОПРЕДЕЛЕНИЯ СИНГУЛЯРНОГО РАЗЛОЖЕНИЯ С ОЦЕНИВАЕМОЙ ТОЧНОСТЬЮ

01 01 07 — вычислительная математика

АВТОРЕФЕРАТ диссертации на соискание ученой степени кандидата физико-математических наук

Л. Г

Новосибирск — /иш

003174222

Работа выполнена в Институте вычислительных технологий Сибирскою отделения РАН

Научный руководитель доктор технических наук,

профессор

Шурина Элла Петровна

Официальные оппоненты докюр физико-математических наук

профессор

Воскобойников Юрий Евхеньевич

Защита сосюихся 31 октября 2007 1 в 15 ч 00 мин на заседании диссертационного совета Д003 061 01 при Институте вычислиаельной математики и математической хеофизики СО РАН но адресу 630090, I Новосибирск, нр ак М А Лаврентьева 6

С диссертацией можно ознакомиться в чи!альном зале библиотеки ИВМиМГ СО РАН

Автореферат разослан 28 сентября 2007 1

Ученый секретарь диссертационно! о совета доктор физико-математических наук, профессор

кандидат физико-математических наук, доцент

Костин Виктор Иванович

Ведущая ор1анизация Институт математики им С Л Соболева

Сибирского отделения РАН

Кузнецов Ю И

ОБЩАЯ ХАРАКТЕРИСТИКА РАБОТЫ

Актуальность работы. Проблема разработки алгоритмов, имеющих низкую вычислительную сложность и гарантирующих точность при определении спектральных и сингулярных характеристик матриц конечно-разностных и конечно-элементных аппроксимаций дифференциальных операторов, является центральной проблемой современной вычислительной алгебры

Знание спектральных характеристик необходимо при решении задач моделирования электромагнитных полей, в задачах структурной механики, ядерной физики, квантовой химии, при решении задач, описывающих свойства диэлектрических волноводов, при анализе океанографических моделей и во многих других приложениях математической физики Знание нулевых сингулярных чисел произвольной матрицы позволяет определить ее ранг, а также необходимо для построения обобщенного нормального решения систем уравнений, у которых число уравнений не совпадает с числом неизвестных Сингулярное разложение (БУЮ) используется для распознавания, сжатия и восстановления изображений, а в генетике и биологии позволяет делать выводы о эволюции генов и связях между протеинами

Сложность определения спектральных и сингулярных характеристик матриц конечно-разностных и конечно-элементных аппроксимаций дифференциальных операторов обусловлена сложными спектральными свойствами, большим размером и разреженной структурой этих матриц Не менее важной является проблема реализации методов определения спектральных и сингулярных характеристик матриц дискретных аппроксимаций дифференциальных операторов в виде библиотек подпрограмм и прикладных пакетов для современных вычислительных систем

Цель работы Целью диссертации является разработка методов определения спектральных характеристик симметричных вещественных матриц и сингулярных характеристик несимметричных вещественных матриц с оцениваемой точностью и реализация этих методов в виде комплекса программ, предназначенного для решения задач моделирования физических процессов современными вычислительными средствами Под реализацией методов решения задач на собственные значения с оцениваемой точностью подразумевается реализация этих методов в стандартной IEEE-754 модели арифметики, основанная на использовании собственных интервалов и, при необходимости, двусторонних последовательностей Штурма

Защищаемые положения. На защиту выносятся

• метод обратной итерации с оцениваемой точностью, предназначенный для расчета заданного числа собственных векторов симметричных вещественных трехдиагональных и трехдиагонализован-ных матриц и сингулярных векторов несимметричных вещественных двухдиагональных и двухдиагонализованных матриц,

в метод Ланцоша с оцениваемой точностью, предназначенный для расчета заданного числа различных собственных значений и собственных векторов симметричных вещественных разреженных матриц большого размера, а также для расчета заданного числа различных сингулярных чисел и сингулярных векторов несимметричных вещественных разреженных матриц большого размера,

• библиотека ANSI С-программ Sturm, модули которой реализуют методы обратной итерации, Ланцоша и бисекций с оцениваемой точностью, предназначенные для определения спектральных и сингулярных характеристик матриц, возникающих в дискретных аппроксимациях задач моделирования физических процессов

Научная новизна. В методе обратной итерации с оцениваемой точностью при выборе сдвига используются собственные интервалы - наименьшие машинно-представимые интервалы, гарантированно содержащие собственные значения рассматриваемой матрицы Формулируется и доказывается теорема о выборе сдвига в методе обратной итерации с оцениваемой точностью Расчет собственных интервалов производится методом бисекций с оцениваемой точностью, представляющем собой новую реализацию метода бисекций в сочетании с QR-методом При определении неколлинеарной системы векторов начального приближения, соответствующих почти кратным и численно кратным собственным значениям, используется новая вычислительная схема неполного спектрального исчерпывания Вычислительная схема спектральной рекурсии используется для расчета вектора начального приближения, соответствующего изолированному собственному значению

В методе Ланцоша с оцениваемой точностью собственные интервалы используются для корректной идентификации истинных собственных значений Метод обратной итерации с оцениваемой точностью применяется для оценки точности собственных и сингулярных чисел, а также при расчете собственных и сингулярных векторов Частичная двухдиагонализация Голуба-Кахана-Ланцоша без реортогона-лизации применяется при расчете сингулярных характеристик несимметричных матриц

Практическая значимость работы. Алгоритмы, описывающие методы обратной итерации, Ланцоша и бисекций с оцениваемой точностью, были реализованы в виде программных модулей ANSI С библиотеки Sturm, реализованной в стандартной модели арифметики IEEE-754 с двойной точностью для операционной системы Linux Библиотека Sturm является отчуждаемым программным комплексом, предназначенным для определения с оцениваемой точностью спектральных ха-

рактеристик симметричных вещественных матриц и сингулярных характеристик несимметричных вещественных матриц дискретных аппроксимаций задач моделирования физических процессов Вычислительный эксперимент демонстрирует высокую эффективность и точность программных модулей библиотеки Sturm

Достоверность полученных результатов проверялась в ходе сравнительного анализа результатов решения задач моделирования физических процессов в библиотеке Sturm с результатами решения этих же задач в библиотеке LAPACK и пакете Matlab Тестирование проводилось на матрицах из наборов данных Matrix Market, Harwell-Boemg, NEP, SPARSKIT и TOKAMAK, возникающих в конечно-разностных и конечно-элементных аппроксимациях задач моделирования конвективно-диффузионных процессов, задач моделирования электромагнитных полей в ускорителях элементарных частиц и в диэлектрических волноводах, а также на серии матриц большой размерности, генерируемых при решении задачи моделирования электрического поля в нефтяной скважине векторным методом конечных элементов

Апробация работы. Результаты исследований докладывались автором на семинаре С К Годунова 'Математика в приложениях' в ИМ СО РАН, на семинарах ИВМиМГ СО РАН и ИВТ СО РАН, на Всероссийской конференции по вычислительной математике 'КВМ-07' (Новосибирск, Россия, 2007г), на минисимпозиуме 'Последние достижения в плотной линейной алгебре' в рамках конференции 'Современный уровень развития научных и параллельных вычислений' (Умеа, Швеция, 2006г) по приглашению организаторов минисимпозиума, на 13-ой Конференции международного общества линейной алгебры (Амстердам, Нидерланды, 2006г), в форме доклада по запросу группы разработчиков библиотеки LAPACK, на Восьмой конференции по ите-

рационным методам (Коппер Маунтейн, Колорадо, США, 2004г), на Восьмой конференции общества индустриальной и прикладной математики (Уильямсбург, Вирджиния, США, 2003г), на Шестом международном симпозиуме по итерационным методам в научных вычислениях (третье призовое место в конкурсе статей аспирантов, Денвер, Колорадо, США, 2003г), на Международной конференции молодых ученых по математическому моделированию и информационным технологиям (Новосибирск, Россия, 2002г), на Конференции молодых ученых СО РАН, посвященной М М Лаврентьеву (Новосибирск, Россия, 2001 г), на Международной конференции 'Современные проблемы прикладной математики и механики теория, эксперимент и практика', посвященной 80-летию академика Н Н Яненко (Новосибирск, Россия, 2001 г), на Конференции молодых ученых, посвященной 10-летию ИВТ СО РАН (Новосибирск, Россия, 2000г), на XV конференции по интервальной математике в рамках конференции 'Вычислительные технологии 2000' (ИВТ СО РАН, Новосибирск, Россия, 2000г)

Публикации. По материалам диссертации опубликовано 10 работ статья в рецензируемом журнале, рекомендованном ВАК для представления докторских диссертаций [1], статья в зарубежном рецензируемом журнале [2], статья в рецензируемом журнале СО РАН [3], публикации в трудах международных и российских конференций [4, 5, 6, 7, 8, 9, 10]

Личный вклад автора. В работе [1] автор участвовал в постановке задачи Автором лично разработаны и реализованы схема неполного спектрального исчерпывания, методы Ланцоша, бисекций и обратной итерации с оцениваемой точностью, проведен вычислительный эксперимент и написан первый вариант статьи В работах [2, 3] автором лично разработан и реализован метод Годунова-обратной итерации, проведен вычислительный эксперимент и написана статья В публикациях [4, 5, 6, 8, 10] автором лично разработаны и реализова-

ны представленные вычислительные схемы и алгоритмы, проведен вычислительный эксперимент, описаны результаты В публикациях [7, 9] автор участвовал в постановке задачи Автором лично разработаны и реализованы представленные алгоритмы, проведен вычислительный эксперимент, описаны результаты

Структура и объем диссертации. Диссертация состоит из введения, четырех глав, заключения, приложения и списка литературы Работа изложена на 154-х страницах, содержит 12 рисунков и 19 таблиц Список литературы содержит 94 наименования

СОДЕРЖАНИЕ РАБОТЫ

Во введении обоснована актуальность темы, указана цель и научная новизна исследования Вводятся следующие обозначения К™ - п-мерное вещественное евклидово пространство с фиксированным ортонормированным базисом (ео О, вектор х = (а:0, а;1, ,хп~1)т € К", скалярное произведение (х, у) = х%Уг1 х е У £ 11ЖИ = V(ж> х) - норма вектора х, ||Л|| = (шах|Л(Лт А)!)1/2 - спектральная норма матрицы А € Япхп, где \{АТ А) - это собственное значение матрицы Ат А, етась - относительная погрешность машинного представления единицы

В первой главе вводятся основные понятия и определения, проводится анализ прямых и итерационных методов определения спектральных и сингулярных характеристик, рассматриваются особенности их реализации Описывается проблема Уилкинсона (1965), возникающая при расчете собственных значений в конечной арифметике Проблема Уилкинсона была впервые решена в 1983 году С К Годуновым и др , разработавшими теорию двусторонних последовательностей Штурма, элементы которых.определяются для собственных интервалов Хг € [А~, А*], симметричной трехдиагональной мат-

рицы Т е Rnxn где (Л+ — Аг ) < emaCh ||Т||, следующим образом Ро(Аг), ,РП_2(Л.) d¿f Р+(Аг"), Р+(Лг-),Р(;1(Л+), Р"^), где Р+(А~) - это левосторонняя последовательность Штурма, Р~(А*) -правосторонняя последовательность Штурма, а целое число / выбирается так, чтобы (Р+(АГ) - Pf(A+))(l/P;-г(А+) - 1/Р++1(К)) < О

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

Т Xi — Х-1

определения собственных значений Аг € R и собственных векторов хг ф 0, хг € Кп, 0 < г < п симметричной трехдиагональной матрицы Т € Rnxn с элементами do, 1 на главной диагонали и с элементами &о, , bn-2 на сопутствующих диагоналях Для расчета собственных интервалов [А~, А^] э Хг матрицы Т применяется новая реализация метода бисекций (Гивенс (1954), Уилкинсон (1965), Годунов и др (1988)), называемая методом бисекций с оцениваемой точностью Собственные интервалы [А~, А*] рассчитываются с погрешностью ег =f (А+ - Ä~) < £rnach 9Л(Т), где Ш(Т) = max{|<¿0| + Ь0|, max {К| + |Ь,|+|Ь1+1|}, |dn-i| + |b„-i|}, причем Аг =f (А++Аг~)/2

1<г<п-1

При расчете менее [гг/10] собственных интервалов метод бисекций с оцениваемой точностью совпадает с методом бисекций в реализации С К Годунова и др (1988) с гарантированной точностью При расчете [гг/10] и более интервалов [А~, А*] выполняется несколько шагов метода бисекций с гарантированной точностью для уточнения интервалов, [Хг — 2 emach ЯЯ(Т'), Хг + 2£тась9Я(Г)], где Аг вычисляются QR методом в реализации, исключающей извлечение квадратных корней, что

в несколько раз уменьшает вычислительную сложность метода бисек-ций

Собственные векторы хг матрицы Т рассчитываются в соответствии с алгоритмом 1, реализующим метод обратной итерации с оцениваемой точностью, который представляет собой новую реализацию метода обратной итерации (Уийландт (1944), Уилкинсон (1965), Ипсен (1997))

Алгоритм 1 Алгоритм вычисления собственных векторов хг, г = I, , то, 1 < I < т < п симметричной трехдиагональной матрицы Т € М"хп с элементами do, , cZ„-i на главной диагонали и с элементами bo, , 6„_2 на сопутствующих диагоналях

for (г = 0, , п - 2)

if (\Ьг | < О Ьг = ^(M^Lch

Выбрать векторы начального приближения € К™, г — I, ,тп пользуясь формулами (1) и (2) Заменить компоненты х[°\ непредставимые в выбранной модели арифметики, случайными числами

for (г = I, I + 1, , т)

if (к > 0 П |схг - cr^ij < emach ЩТ))

then <тг = шах {max {А~, <тг_1 + ег/2}, mm {crt_i + ег 2/3, Л^1"}} else <тг =

Факторизовать Т — аг I = L D LT while (|k(<:)|| < 1 Пк < 5)

Решить LDLT z[k+1) = х[к), к = к + 1

if (|<7"г—1 - | <

then for (j = 0,1, , г - 1) z[k) = z[k) - (x}, z<fc)) x3 if (K(fc)||>l) then

else Выдать сообщение 'сходимость не достигнута' В алгоритме 1 схема спектральной рекурсии

*? = 1, a:f = -x1fc-1sign(bh_1)/ft_1(A,), (1)

где к = 1,2, , п — 1, а (Рп(Аг)} - двусторонняя последовательность Штурма, применяется при определении начального вектора, соответствующего изолированному собственному значению Аг |Аг — Аг+11 > £mach ||Г|| Для построения системы неколлинеарных начальных векторов, соответствующих почти кратным собственным значениям Хэ, j > +1, , для которых |Аг — 1 < £mach ||Г||, j > г, предлагается новая вычислительная схема неполного спектрального исчерпывания

хп-1 = с°е„_ 1

хп_2 = С1 еп_2 + (О, ,0,<_2)т xn^_k = Cfcen_i_fc + (0, , 0, uj,_i-k)r (2)

хо = Сп~1 во + (0, и®, Х"3.«Г2)Т.

где числа uJk являются числами из случайного равномерного на (0,1) распределения, а Ск = Сп_2 С'п_з Со являются цепочками элементарных преобразований вращения С;, параметры которых определяются по элементам двусторонней последовательности Штурма Система (2) является грубой аппроксимацией ортогонального базиса матрицы Т, вычисляемого С К Годуновым и др (1988) методом спектрального исчерпывания

Теорема 1 Пусть сдвиг сгг в алгоритме 1 выбран таким образом, что о» € [Л~, Л+], причем аг ф Лг, Лг € [Аг~, А+], и ег = (Лг+ - А") < ¿шась ЩТ) Пусть = тгг{к)/\\г[к)\\, тг = ег и Ц^^Ц > 1 Тогда

||(Т - <тг 7)_1|| > т,/||(Т - <т, I) > 1/е.

Теорема 1 означает, что по достижении сходимости в алгоритме 1 рассчитанный собственный вектор х[к^ образует пару с близким к собственному числу Аг числом аг из -спектра || (Т—аг /)~1Ц > 1/£г матрицы Т Назовем <тг псевдособственным значением, х[к^ - псевдосбстве-ным вектором, а (сгг, з;^) - псевдособственной парой

Третья глава посвящена описанию метода Ланцоша с оцениваемой точностью, представляющего собой новую реализацию метода Коллум-Уилуби-Ланцоша (Коллум и Уилуби (1979, 1981, 2002)) Рассматриваются вопросы ускорения сходимости при помощи преобразования сдвиг-обращение (симметричные матрицы) и обращения (невырожденные несимметричные матрицы) для разреженных матриц Представлены алгоритмы, реализующие метод Ланцоша с оцениваемой точностью в двух вариантах для определения различных собственных значений и, при необходимости, собственных векторов симметричных вещественных разреженных матриц, а также для определения различных сингулярных чисел и, при необходимости, сингулярных векторов несимметричных вещественных разреженных матриц Проводится анализ вычислительной сложности разработанных алгоритмов

В основе метода Коллум-Уилуби-Ланцоша лежит следующий эмпирический факт ('феномен Ланцоша') при достаточно большом числе шагов в методе Ланцоша без реортогонализации каждое собственное значение симметричной матрицы с собственным вектором, имеющим нетривиальную проекцию на начальный вектор, является собственным значением генерируемой симметричной трехдиагональной

матрицы Ланцоша (Годунов и Прокопов (1970), ван Катц и ван дер Воет (1976), Коллум и Уилуби (1981, 2002)) Во избежание ошибочной идентификации почти кратных собственных значений, как ложных, в методе Ланцоша с оцениваемой точностью предлагается использовать собственные интервалы [А~, При проверке наличия у матрицы Тк ложных, близких к Аг собственных значений, предлагается рассматривать интервал [А~ —£г, +£г], а при идентификации кластеров из близких к Аг чисел Х3 предлагается использовать критерий I Aj — Аг| < maxtjJ {ег, е3 } При проверке критерия сходимости собственные векторы хг предлагается рассчитывать методом обратной итерации с оцениваемой точностью Определение наибольших сингулярных чисел в методе Ланцоша с оцениваемой точностью предлагается проводить непосредственным применением метода Голуба-Кахана-Ланцоша без реортогонализации для двухдиагонализации несимметричной матрицы А 6 RnXn, с последующим анализом симметричной трехдиагональной матрицы Голуба-Кахана

В четвертой главе дается описание вычислительных схем, реализованных в библиотеке Sturm, представлены таблицы и рисунки, иллюстрирующие результаты вычислений, проводится сравнительный анализ результатов тестирования программных модулей библиотеки Sturm с результатами аналогичных вычислений, выполненных в библиотеке LAPACK 3 0/3 1 и в пакете Matlab 7 2 Тестирование проводилось в вычислительной среде со следующими параметрами компьютер -четырех-процессорный Intel Xeon 3 20 GHz с 4 0 GB общей памяти, операционная система - Red Hat Linux 3 2 3, компилятор - gcc 3 2 3, уровень оптимизации - стандартный (-ОЗ), таймеры — clock (Sturm и LAPACK), cputime (Matlab) Логические связи между основными вычислительными схемами, реализованными в библиотеке Sturm, представлены на рисунке 1, где SVD - это сингулярное разложение несимметричной вещественной квадратной матрицы, SSD - спектральное

разложение симметричной вещественной матрицы, X- трехдиаго-нальная симметричная вещественная матрица, В к - двухдиагональная несимметричная вещественная матрица, а Gk - соответствующая Вк матрица Голуба-Кахана. Вычислительные схемы Sturm-LGAT, Sturm-LTCW и Sturm-LTICW реализуют метод Коллум-Уилуби-Ланцоша в симметричном случае, а Sturm-LBAT, Sturm-LBCW и Sturm-BTICW -в несимметричном. Собственные значения матрицы Тк (Gk) находятся методом бисекций с оцениваемой точностью по схемам Sturm-BSC и Sturm-QBSC. Собственные векторы матрицы Тк (Gk) рассчитываются методом обратной итерации с оцениваемой точностью в соответствии с одной из следующих схем, отличающихся выбором начальных векторов: Sturm-QIH, Sturm-QID, Sturm-QIR.

Рис. 1: Логические связи между основными вычислительными схемами, реализованными в библиотеке Sturm, и вспомогательными программными модулями, вызываемыми из библиотек LAPACK и SuperLU.

Задачи моделирования электромагнитного поля в нефтяной скважине. Рассмотрим аппроксимацию 3-х мерного уравнения Гельмгольца векторным методом конечных элементов на гексаэдраль-ной сетке в следующей вариационной постановке (Шурина и Нечаев (2005))

Найти Eh е Hj(rot, Q) такое, что дляУ vh € Hft(rot, П)

1 (3)

(— rotE*1, rotvh) + (k2Eh,vh) = -i(wJ0,v'1) Соответствующая система линейных уравнений имеет вид

( D + B -С \ ( Ег\^( fr\ \ С D + B )\Et) \ /, )'

л

где D + В = (D + В)Т G Rnx", С = Ст е Rnxn, d%3 = (¿rotN.\Nj), Ъг} = -(еа;2Г^,Г^), е., = (trwNj, N*), N^ - это векторные локальные базисные функции, ассоциированные с ребрами гексаэдра, а Ег и Ег - веса в разложении по базису действительной и мнимой компонент поля Е Рассмотрим задачу (3) в непроводящей области П = [—1,1]3 в которой находится генератор тока, при этом ш = 14 MHz, I = 1А, е = 8 85е - 12Ф/м В этом случае а = 0, а А = Ат

= D + В В таблице 1 представлены результаты расчета частичного спектрального разложения для пяти наименьших по модулю собственных значений матрицы А порядка п = 26460 Расчет проводился по схеме Sturm-LTICW в библиотеке Sturm и пакете Matlab при помощи процедуры EIGS по схеме Matlab-EIGS Время t\yX расчета по схеме Sturm-LTICW оказалось приблизительно в двадцать раз быстрее времени расчета по схеме Matlab-EIGS Более того, одно из собственных значений было рассчитано Matlab-EIGS неверно — характер векторной вариационной постановки (3) исключает наличие нулевых

собственных значений у матрицы А

Таблица 1 Частичное спектральное разложение для пяти различных наименьших по величине собственных значений матрицы А на сетке с шагом к = 0 05, (сЬт(Л) = 26460)

Лтт max \\Ахк - А к xfc|| max хк - 1|| «л,х

Sturm-LTICW -6 7404329919582245e - 07 -6 7404392292633789e - 07 -6 7404454665203500e - 07 -6 7404517037293810e - 07 -7 9996261463491110e - 07 1 46e — 11 1 83e — 07 343 66

Matlab-EIGS 0 0000000000000000e + 00° -6 7404329919552969e - 07 -6 7404392292524194e - 07 -6 7404454665091893e - 07 -6 7404517037259940e - 07 8 OOe — 07 1 81e — 13 6879 48

"неверное решение

В заключении представлено краткое описание полученных результатов разработаны методы обратной итерации, Ланцоша и би-секций с оцениваемой точностью, разработана ANSI С библиотека Sturm, программные модули которой реализуют методы с оцениваемой точностью, проведен вычислительный эксперимент, в ходе которого программные модули библиотеки Sturm применялись для расчета спектральных и сингулярных характеристик матриц, возникающих в конечно-разностных и конечно-элементных аппроксимациях задач моделирования физических процессов, сравнительный анализ полученных результатов с результатами решения этих же задач в библиотеке LAPACK и в пакете Matlab демонстрирует высокую эффективность и точность программных модулей библиотеки Sturm

В приложении А представлено описание основных программных модулей библиотеки Sturm

Основное содержание диссертации опубликовано в следующих работах-

[1] А М Мацех, Э П Шурина Нахождение спектрального разложения симметричных матриц и сингулярного разложения несимметричных матриц с оцениваемой точностью Автометрия, 43(2) 81-96, 2007

[2] А М Matsekh The Godunov-inverse iteration a fast and accurate solution to the symmetric tridiagonal eigenvalue problem Applied Numerical Mathematics, 54(2) 208-221, 2005

[3] A M Matsekh The Godunov-inverse iteration algorithm for symmetric tridiagonal matrices Bulletin of the Novosibirsk Computing Center Series Computer Science, (19) 39-50, 2003

[4] A M Matsekh Inverse iteration with guaranteed accuracy - a new method for computing eigenvectors of real symmetric tridiagonal matrices Program and Abstracts of the Minisymposium on Recent Advances in Dense Linear Algebra, Workshop on State-of-the-Art in Scientific and Parallel Computing, 72 Ume&, Sweden, 2006

[5] A M Matsekh Lanczos method with guaranteed accuracy - a new implementation of the Lanczos method with no reorthogonahzation Program and Abstracts of the 13-th International Linear Algebra Society Conference, 81-82 Amsterdam, Netherlands, 2006

[6] A M Matsekh Using spectral deflation to accelerate convergence of inverse iteration for symmetric tridiagonal eigenproblems Program and Abstracts of the Eighth Copper Mountain Conference on Iterative Methods, March 28-April 2, 61-62 Copper Mountain, CO, USA, 2004

[7] A M Matsekh, E P Shurma Using Godunov's two-sided Sturm sequences to accurately compute singular vectors of bidiagonal matrices Program and Abstracts of the The Eighth SIAM Conference on Applied Linear Algebra, 33 Williamsburg, VA, USA, 2003

[8] A M Мацех Анализ алгоритмов вычисления сингулярных чисел и сингулярных векторов несимметричных матриц по соответствующим матрицам в форме Голуба-Кахана Программа и тезисы докладов международной конференции молодых ученых по математическому модели-

рованию и информационным технологиям, 30 Академгородок, Новосибирск, Россия, 2002 [9] А М Matsekh, Е Р Shurma A new hybrid procedure for computing eigenvectors of tridiagonal symmetric unreduced matrices Вычислительные Технологии, 6, 2001 Special Issue, Part 2 Proceedings of the International conference 'Recent Developments m Applied Mathematics and Mechanics Theory, Experiment and Practice' [10] A M Мацех Вычислительная схема для решения полной проблемы собственных значений для плотных симметричных матриц Труды конференции молодых ученых, посвященной 10-ти летию ИВТ СО РАН, Т II математическое моделирование, 108 - 111 Новосибирск, Россия, 2000

Тираж 100 экз Заказ № 721 24 09 07г Отпечатано в ЗАО РИЦ «Прайс-курьер» Проспект ак Лаврентьева, 6 Тел 330-72-02

Оглавление

Введение

Глава 1 Современное состояние исследований

1.1 Алгебраическая проблема собственных значении.

1.2 Сингулярное разложение .1С

1.3 Проблема Уилкинсона.

1.4 Последовательности Штурма.

1.5 Методы определения спектральных и сингулярных характеристик и их реализация.

1.0 Выводы и постановка задачи.

Глава 2 Методы бисекций и обратной итерации с оцениваемой точностью

2.1 Метод бисекций с оцениваемой точностью.

2.2 Метод обратной итерации с оцениваемой точностью.

2.2.1 Определение собственного вектора

2.2.2 Выбор сдвига.

2.2.3 Определение начальных векторов

2.2.4 Расчет нескольких собственных векторов.

2.3 Выводы.

Глава 3 Метод Ланцоша с оцениваемой точностью

3.1 Особенности реализации.

3.2 Алгоритм определения собственных чисел

3.3 Особенности расчета сингулярных чисел

3.4 Алгоритм определения сингулярных чисел.

3.5 Выводы.

Глава 4 Библиотека Sturm

4.1 Общая характеристика

4.2 Вычислительные схемы.

4.3 Верификация и вычислительный эксперимент.

4.3.1 Полное спектральное разложение.

4.3.2 Частичное спектральное разложение.

4.3.3 Частичное сингулярное разложение.

4.4 Выводы.

Актуальность работы

Проблема разработки быстрых алгоритмов, гарантирующих точность при определении спектральных и сингулярных характеристик матриц конечно-разностных и конечно-элементных аппроксимаций дифференциальных операторов, является центральной проблемой современной вычислительной алгебры. Знание спектральных характеристик необходимо при решении задач моделирования электромагнитных полей |40, 3|, в задачах структурной механики [5, 21|, ядерной физики, квантовой химии, при решении задач, описывающих свойства диэлектрических волноводов, при анализе океанографических моделей и во многих других приложениях математической физики [47, 59]. Сингулярное разложение (БУБ) используется для распознавания, сжатия и восстановления изображений в компьютерной графике, а в генетике и биологии позволяет делать выводы о эволюции генов и связях между протеинами. Знание нулевых сингулярных чисел произвольной матрицы позволяет определить ее ранг, а также необходимо для построения обобщенного нормального решения систем уравнений, у которых число уравнений не совпадает с числом неизвестных.

Сложность определения спектральных и сингулярных характеристик матриц, возникающих при моделирования физических процессов как правило обусловлена необходимостью определения близких, почти кратных собственных значений и соответствующих им почти коллинеарных собственных векторов, при расчете которых стандартными вычислительными методами наблюдается медленная сходимость, обусловленная потерей ортогональности и точности. Не менее важной является проблема реали-л' зации методов определения спектральных и сингулярных характеристик в виде библиотек подпрограмм и прикладных пакетов для современных вычислительных систем. Первая широко известная библиотека подпрограмм для решения задач на собственные значения - Е1ЭРАСК [70] - была разработана в середине 70-х годов. На смену ей пришла библиотека подпрограмм ЬАРАСК [1, 15], представляющая собой набор процедур, реализующих алгоритмы решения задач линейной алгебры с плотными, трех- и двухдиагональпыми матрицами. Одной из наиболее распространенных современных библиотек предназначенных для определения спектральных и сингулярных характеристик матриц большой размерности является библиотека А11РАСК [45]. Среди современных пакетов прикладных программ, позволяющих решать задачи на собственные значения и задачу определения спектрального разложения, отметим пакет МаЫаЬ [75]. При определении спектрального и сингулярного разложений пакет Ма11аЬ вызывает откомпилированные модули из библиотек ЬАРАСК и А11РАСК, при этом все вычисления производятся в стандартной модели арифметики 1ЕЕЕ-754 [2, 62, 48]. Точность в перечисленных комплексах программ не гарантируется.

В основе современных методов определения спектральных и сингулярных характеристик с гарантированной точностью лежат работы С.К. Годунова, А.Г. Антонова, О.П. Киршпока, В.И. Костина и А.Д. Митчен-ко [84, 83, 28]. В методах с гарантированной точностью используется элементы интервального анализа и специально разработанная модель арифметики с направленным округлением. Помимо приближенных собственных чисел, в методах с гарантированной точностью используются наименьшие машинно-нредставимые интервалы, гарантированно содержащие собственные значения рассматриваемой матрицы. Будем называть такие интервалы собственными интервалами. С.К. Годунов и др. строят численно ортогональную систему приближенных собственных векторов симметричной трехдиагоналыюй вещественной матрицы Т методом спектрального исчерпывания, применяя преобразования вращения, параметры которых строятся при помощи двусторонних последовательностей Штурма [84, 83, 82]. Двусторонняя последовательность Штурма иитериретпруется [60] как решение проблемы Уилкинсона [77] о корректном исключении уравнения из переопределенной системы уравнений Т х — Хх = 0. Индекс исключаемого уравнения в проблеме Уилкинсопа совпадает с индексом исключаемого элемента двусторонней последовательности Штурма. В силу высокой вычислительной сложности [GG] метод спектрального исчерпывания не приемлем для определения собственных и сингулярных векторов разреженных матриц большой размерности, возникающих в аппроксимациях дифференциальных операторов многих задач математической физики. Такие задачи решаются итерационными методами.

Метод Ланцоша [44, 58J является одним из наиболее эффективных итерационных методов решения симметричной проблемы собственных значений. В основе современных формулировок метода Ланцоша лежат работы Пэйджа [63, 04]. Существует несколько формулировок метода Ланцоша: с полной реортогонализациеп [31, 83], с выборочной и частичной реорто-гоиализацией [G], без реортогонализации [11, 37, 12, 80], с квазиминимизацией [27], с неявным перезапуском [71, G]. Метод Ланцоша с неявным перезапуском реализован в библиотеке ARPACK и применяется для определения спектрального и сингулярного разложений больших разреженных матриц. Метод Ланцоша с неявным перезапуском является вариантом метода. Арнольди с неявным перезапуском [71], в котором может происходить сходимость к ложным собственным значениям [23]. Вычислительный эксперимент показывает (см. раздел 4.3), что ложная сходимость возможна и в методе Ланцоша с неявным перезапуском. Реализация метода Ланцоша без реортогонализации в соответствии со схемой Коллум и Уилуби, называемая также методом Коллум-Уилуби-Ланцоша, позволяет идентифицировать ложные собственные значения и наиболее приемлема при работе с большими разреженными матрицами. Однако в методе Коллум-Уилуби-Ланцоша возможна ошибочная идентификация истинных почти кратных собственных значений как ложных, а скорость сходимости определяется спектральными свойствами задачи, в результате чего число внешних итераций может оказаться непредсказуемо большим.

Собственные векторы симметричных трехдиагональных матриц часто определяются методом обратной итерации, который имеет низкую вычислительную сложность порядка п2 при расчете хорошо обусловленного базиса из собственных векторов. При наличии почти кратных собственных значений метод обратной итерации требует дополнительной реортогона-лизации. М11ГШ, метод, или метод относительных робастпых представлений, является повой реализацией метода обратной итерации, в которой при выборе вектора начального приближения используется 'составная факторизация', представляющая собой альтернативное решение проблемы Уил-кинсона. МШ111 метод был впервые реализован без реортогонализации в библиотеке ЬАРАСК 3.0 в виде процедуры хЭТЕСЛ, имеющей сложность 0(п2). Нами установлено [52, 53], что в случаях, когда матрица имеет почти кратные собственные значения, соответствующие им собственные векторы вычисляются процедурой хЭТЕСИ, в ЬАРАСК 3.0 с недопустимо большими погрешностями. В версии ЬАРАСК 3.1 данная проблема решается включением выборочной реортогонализации в новую реализацию МШЖ метода.

В методе обратной итерации используются собственные значения, рассчитываемые методом бисекций, вычислительная сложность которого при расчете всех собств енпых значений составляет 5* ~ 2/с п арифметических операций [18], где к представляет собой число 7-ичных цифр мантиссы. В 1ЕЕЕ-754 модели арифметики с двойной точностью вычислительная сложность метода бисекций составляет Я « 106 п2, что почти в 10 раз медленнее С}11 метода в реализации, исключающей извлечение квадратных корней, вычислительная сложность которого составляет 5 ~ 12 п2 [18| арифметических операций.

Заметим, что проблемы, связанные с оптимальностью и точностью при определении спектральных и сингулярных характеристик, систематически исследовались в работах Хоффмаиа и Уайлаидта [35], Голуба и Кахана [29], Парлетта [05], Стюарта [72, 74], Кахана [13] и Деммела [14], Трефетена и Эмбри [70], Хайэма [33|, Князева [42, 43, 41], Малышева [40, 19]. Из недавних работ, посвященных определению собственных значений матриц большой размерности, отмстим работы Коллум [10, 8], Саада [08] и Сорспсе-на [45, 71].

Цель работы

Целыо диссертации является разработка методов определения спектральных характеристик симметричных вещественных матриц и сингулярных характеристик несимметричных вещественных матрице оцениваемой точностью и реализация этих методов в виде комплекса программ, предназначенного для решения задач моделирования физических процессов современными вычислительными средствами. Под реализацией методов решения задач на собственные значения с оцениваемой точностью подразумевается реализация этих методов в стандартной IEEE-754 модели арифметики, основанная на использовании собственных интервалов и, при необходимости, двусторонних последовательностей Штурма.

Защищаемые положения

На защиту выносятся:

• метод обратной итерации с оцениваемой точностью, предназначенный для расчета заданного числа собственных векторов симметричных вещественных трехдиагональпых и трехдиагонализованных матриц и сингулярных векторов несимметричных вещественных двухдиагопаль-пых и двухдиагопализованных матриц;

• метод Лаицоша с оцениваемой точностью, предназначенный для расчета заданного числа различных собственных значений и собственных векторов симметричных вещественных разреженных матриц большого размера, а также для расчета заданного числа различных сингулярных чисел и сингулярных векторов несимметричных вещественных разреженных матриц большого размера;

• библиотека С-программ Sturm, программные модули которой реализуют методы обратной итерации, Лаицоша и бисекций с оцениваемой точностью, предназначенные для определения спектральных и сингулярных характеристик больших разреженных матриц, возникающих в конечно-элементных и конечно-разностных аппроксимациях задач моделирования физических процессов.

Научная новизна

D методе обратной итерации с оцениваемой точностью

• собственные интервалы используются при выборе сдвига; расчет собственных интервалов производится методом бисекций с оцениваемой точностью, представляющем собой новую реализацию метода бисекций с гарантированной точностью в сочетании с QR-методом;

• формулируется и доказывается теорема о выборе сдвига в методе обратной итерации с оцениваемой точностью, гарантирующее что сдвиг и приближенно рассчитанный собственный вектор образуют псевдособственную пару;

• новая вычислительная схема неполного спектрального исчерпывания используется при определении неколлинеарной системы векторов начального приближения, соответствующих почти кратным и численно кратным собственным значениям;

• вычислительная схема спектральной рекурсии используется для расчета вектора начального приближения, соответствующего изолированному собственному значению;

• расчет сингулярных векторов двухдиагональной вещественной матрицы сводится к расчету собственных векторов соответствующем! симметричной трехдиагопальной матрицы Голуба-Кахаиа.

В методе Ланцоша с оцениваемой точностью

• собственные интервалы используются для корректной идентификации истинных собственных значений; расчет собственных интервалов производится методом бисекций с оцениваемой точностью;

• метод обратной итерации с оцениваемой точностью применяется для оценки точности собственных и сингулярных чисел при расчете критерия останова итерационного процесса;

• метод обратной итерации с оцениваемой точностью применяется при расчете собственных и сингулярных векторов;

• частичная двухдиагонализация Голуба-Кахаиа-Ланцоша без реортогона-лизации применяется при расчете сингулярных характеристик несимметричных матриц, после чего решается эквивалентная частичная проблема собственных значений с матрицей Голуба-Кахаиа.

Практическая значимость работы

Алгоритмы, описывающие методы обратной итерации, Ланцоша и бисекций с оцениваемой точностью, были реализованы в виде программных модулей ANSI С библиотеки Sturm, реализованной в стандартной модели арифметики IEEE-754 с двойной точностью для операционной системы Linux. Библиотека Sturm является отчуждаемым программным комплексом, который может использоваться либо самостоятельно, либо может интегрироваться в другие библиотеки и комплексы программ. Библиотека Sturm предназначена для определения с оцениваемой точностью спектральных характеристик симметричных вещественных матриц и сингулярных характеристик несимметричных вещественных матриц дискретных аппроксимаций задач моделирования физических процессов. Вычислительный эксперимент демонстрирует высокую эффективность и точность программных модулей библиотеки Sturm.

Достоверность полученных результатов

Достоверность полученных результатов проверялась в ходе сравнительного анализа результатов решения задач моделирования физических процессов в библиотеке Sturm с результатами решения этих же задач в библиотеке LAPACK и пакете Matlab. Тестирование проводилось на матрицах, возникающих в конечно-разностных и конечно-элементных аппроксимациях задач моделирования конвективно-диффузионных процессов, на матрицах из наборов данных Matrix Market [59|, Harwell-Boeing [20|, NEP [7|, SPARSKIT [G9] и ТОКАМАК [59], возникающих в дискретных аппроксимациях задач структурной механики, океанографии, ядерной физики и электромагнетизма, а также на серии матриц большой размерности, генерируемых при решении задачи моделирования электрического ноля в нефтяной скважине векторным методом конечных элементов [94].

Апробация работы

Результаты проведенных исследований докладывались автором на семинаре С.К. Годунова 'Математика в приложениях' в ИМ СО РАН, на семинарах ИВМиМГ СО РАН и ИВТ СО РАН, на Всероссийской конференции по вычислительной математике 'КВМ-07' (Академгородок, Новосибирск, Россия, 2007г.) [93], на мииисимпозиуме 'Последние достижения в плотной линейной алгебре' в рамках конференции 'Современный уровень развития научных и параллельных вычислений' (Умеа, Швеция, 2006 г.) |54] по приглашению организаторов минисимпозиума; на 13-ой Конференции международного общества линейной алгебры (Амстердам, Нидерланды,

2000 г.) (55]; в виде технического доклада [53] по запросу группы разработчиков библиотеки ЬАРАСК; на Восьмой конференции по итерационным методам (Коппер Маунтейн, Колорадо, США, 2004 г.) [51]; на Восьмой конференции общества индустриальной и прикладной математики (Уильямс-бург, Вирджиния, США, 2003 г.) [57]; на Шестом международном симпозиуме но итерационным методам в научных вычислениях (третье призовое место в конкурсе статей аспирантов, Денвер, Колорадо, США, 2003 г.) [49]; па Международной конференции молодых ученых по математическому моделированию и информационным технологиям (Академгородок, Новосибирск, Россия, 2002 г.) [90]; па Конференции молодых ученых СО РАН, посвященной М.М. Лаврентьеву (Академгородок, Новосибирск, Россия,

2001 г.) [89]; на Международной конференции 'Современные проблемы прикладной математики и механики: теория, эксперимент и практика', посвященной 80-летию академика Н. Н.Яиепко (Академгородок, Новосибирск, Россия, 2001г.) [50]; на Конференции молодых ученых, посвященной 10-летию ИВТ СО РАН (Академгородок, Новосибирск, Россия, 2000г.) [88]; па XV конференции по интервальной математике в рамках конференции 'Вычислительные Технологии 2000' (ИВТ СО РАН, Академгородок, Новосибирск, Россия, 2000 г.) [91].

Публикации

По материалам диссертации опубликовано 10 работ, в том числе статья 15 рецензируемом журнале, рекомендованном ВАК для представления докторских диссертаций |92], статья в зарубежном рецензируемом журнале [52], статья в рецензируемом журнале СО РАН [50], а также; публикации [54, 55, 51, 57, 90, 50, 88] в трудах международных и российских конференций.

Личный вклад автора

В работе [92] автор участвовал в постановке задачи, автором лично разработаны и реализованы схема неполного спектрального исчерпывания, методы Ланцоша, бисекций и обратной итерации с оцениваемой точностью, проведен вычислительный эксперимент и написан первый вариант статьи. В работах [52, 50] автором лично разработан и реализованы метод Годунова-обратной итерации, проведен вычислительный эксперимент и написана статья. В публикациях [54, 55, 51, 90, 88] автором лично разработаны и реализованы представленные вычислительные схемы и алгоритмы, проведен вычислительный эксперимент, описаны результаты. В публикациях [57, 56] автор участвовал в постановке задачи, автором лично разработаны и реализованы представленные алгоритмы, проведен вычислительный эксперимент, описаны результаты.

Структура работы

Во введении обосновывается актуальность темы, указывается цель и научная новизна исследования, вводятся обозначения. В первой главе проводится обзор современных исследований, дается формальное определение основных понятий, проводится анализ современных методов определения спектральных характеристик симметричных трехдиагопальных матриц и методов подпространства Крылова. Вторая глава посвящена описанию повой реализации методов бисекций и обратной итерации с оцениваемой точностью, представлено детальное описание алгоритмов, формулируется и доказывается теорема о выборе сдвига в методе обратной итерации с оцениваемой точностью, описываются схемы спектральной рекурсии и неполного спектрального исчерпывания. В третьей главе; представлена новая реализация метода Ланцоша с оцениваемой точностью в двух вариантах - для определения произвольного числа различных собственных значений симметричных вещественных матриц и для определения произвольного числа сингулярных чисел несимметричных вещественных матриц. В четвертой главе описываются вычислительные схемы, реализованные в библиотеке Sturm, представлены результаты вычислительного эксперимента. В заключении дано краткое описание полученных результатов. Описание основных программных модулей библиотеки Sturm представлено в приложении А.

Обозначения

• С" п-мерное комплексное евклидово пространство;

• Ж" п-мерное вещественное евклидово пространство;

• (со, С],., еп-\) - фиксированный ортопормированный базис в1", где еЛ = (0,.,0,^0,.,0)г, А; = 0,1,., п — 1; к+1

• х - (ж0, х\ ., хп~1)т - вектор х € (х Е Сп);

• (ж, у) = жт у скалярное произведение ж £ С" и у Е С";

• (ж, у) = хт у - скалярное произведение я Е К" и у Е К";

• ||а;|| = у/(х,у) - норма вектора х Е М™ (х Е С");

• С"хп - п-мерное комплексное векторное пространство;

• Мпхп - п-мерное вещественное векторное пространство;

• А = (а^), г,^ — 0,1,. ,п - 1 - матрица А Е М,гхп (А Е С'1ХИ);

• Ат = (ау?:), г, ,7 = 0,1,., п — 1 - транспонированная матрица А Е Шпхп ■>

• А* — (а^г), г, ,7 = 0,1,., п — 1 - комплексно-сопряженная транспонированная матрица А Е Спх?г;

• с^^ - диагональная матрица;

• I - единичная матрица в ШпХп (С/гхп)

• \iiii = Ао < А1 < • • • < Л,г1 = Ашах - собственные значения матрицы А Е Епхп {А Е Сих");

• 0 < (тт\п = сто < о"1 < - • • < <7„1 = <тшах - сингулярные числа матрицы

ЛеМ"*" (Ле€'да1);

• ||Л|| = о"пшх(Л) - спектральная норма матрицы А Е К7гх" (А Е Спхп);

• ^шасЬ ~ относительная погрешность машинного представления единицы в выбранной модели конечной арифметики;

• £шш ~ минимальное положительное машинное число в выбранной модели конечной арифметики;

• £шах ~ максимальное машинное число в выбранной модели конечной арифметики.

4.4 Выводы

Методы бисекций, обратной итерации и Ланцоша оцениваемой точностью реализуются в библиотеке Sturm соответствии с разработанными обобщенными вычислительными схемами, предназначенными для решения задач, требующих определения полного и частичного спектрального разложения симметричных и сингулярного разложения несимметричных вещественных матриц, хранящихся как в плотных, так и в разреженных форматах. Сравнительный анализ результатов тестирования иллюстрирует высокую точность и эффективность методов с оцениваемой точностью.

Заключение

• Разработан метод обратной итерации с оцениваемой точностью, представляющий собой новый вариант метода обратной итерации, предназначенный для расчета заданного числа собственных векторов симметричных вещественных трехдиагональпых матриц и сингулярных векторов вещественных двухдиагональных матриц.

• Сформулирована и доказана теорема о выборе сдвига в методе обратной итерации с оцениваемой точностью, гарантирующее что сдвиг и приближенно рассчитанный собственный вектор образуют псевдособственную пару.

• Разработай метод Ланцоша с оцениваемой точностью, представляющий собой новый вариант метода Ланцоша без реортогоиализации, предназначенный для расчета заданного числа различных собственных значений симметричных вещественных матриц и сингулярных чисел несимметричных разреженных вещественных матриц большой размерности.

• Разработан метод бисекций с оцениваемой точностью, представляющий собой новый вариант метода бисекций с гарантированной точностью в сочетании с QR-методом.

• Разработана отчуждаемая ANSI С библиотека Sturm, программные модули которой реализуют методы бисекций, обратной итерации и Ланцоша с оцениваемой точностью в модели арифметики IEEE-754.

• Проведено тестирование основных программных модулей библиотек Sturm, LAPACK и пакета Matlab па матрицах возникающих в конечно-разностных и конечно-элементных аппроксимациях задач моделирования тепловых и электромагнитных полей, задач структурной механики, океанографии, ядерной физики. Вычислительный эксперимент демонстрирует высокую эффективность и точность программных модулей библиотеки Sturm.

1. E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney, D. Sorensen. LAPACK Users' Guide. S1.M, Philadelphia, third edition, 1999.

2. ANSI/IEEE. IEEE Standard for Binary Floating Point Arithmetic. New York, 7541985 edition, 1985.

3. P. Arbenz, R. Geus, S. Adam. Solving Maxwell eigenvalue problems for accelerating cavities. Physical Review Special Topics Accelerators and Beams, 4(022001), 2001.

4. Susanne Balk;, Jane Culluin. A parallel algorithm for c, omputing eigenvalues of very large real symmetric matrices on message passing architectures. Applied Numeriail Mathematics, 30:341 3G5, 1999. Theory.

5. F.L. Bauer, C.T. Fike. Norms and exclusion theorems. Numerische Math., 2:137 141, 1960.

6. Jane K. Cullum, Ralph A. Willoughby. Lanczos algorithms for large symmetric eigenvalue computations, volume 1. Society for Industrial and Applied Mathematics, Philadelphia, 2002.

7. J. Demmel, W. Kalian. Accurate singular values of bidiagonal matrices. SIAM J. Sci. Stat. Cornput., 11:873-912, 1990.

8. James Demmel. Accurate singular values decompositions of structured matrices. SIAM J. Matrix Anal. Appl, 21:502 580, 1999.

9. I.S. Dhillon, A.N. Malyshev. Inner deflation for symmetric tridiagonal matrices. Linear Algebra and its Applications, 358:139-141, 2003.

10. Juan C. Egana, Nelson M. Kuhl, Luis C. Santos. An inverse eigenvalue method for frequency isolation in spring-mass systems. Numerical Linear Algebra with Applications, 9(l):65-79, January/February 2002.

11. Siever Ellen, Inc. the staff of O'Reilly & Associates. Linux in a nutshell. O'Reilly k Associates, Inc., second edition, 1999.

12. Mark Embree. Misconvergence of Arnoldi Eigenvalue Iterations. 13-th International Linear Algebra Society Conference. Amsterdam, Netherlands, 2006. URL http:// staff.science.uva.nl/~brandts/ILAS06/.

13. J.G.F. Francis. The qr transformation: a unitary analogue to the lr transformation, parts i and ii. Com,put. ,/., 4:205-271, 332 345, 1961.

14. G. Golub, J. Wilkinson. Ill-conditioned eigensystems and the computation of the Jordan canonical form. SIAM Review, 18:578 619, 1976.

15. Gene H. Golub, Charles F. Van Loan. Matrix Computations. The Johns Hopkins University Press, third edition, 1996.

16. G.H. Golub, H.A. van der Vorst. Eigenvalue computation in the 20th century. Journal of Computational and Applied Mathematics, 123(1—2):35 65, 2000. Nicholas J. Higham. Accuracy and stability of numerical algorithms. SIAM, second edition, 2002.

17. R. Hiptmair. Finite elements in computational electromagnetism. Acta Numerica, 11:237- 339, 2002.

18. A.J. Hoffman, H.W. Wielandt. The variation of the spectrum of a normal matrix. Duke Math. J., 20:37 40, 1953.

19. Use C.F. Ipsen. Computing an eigenvector with inverse iteration. SIAM Review, 39(2):254-291, 1997.

20. Cullurn Jane, Willoughby Ralph A. Computing eigenvalues of very large symmetric matrices an implementation of a Lanczos algorithm with no reorthogonalization. Journal of Computational Physics, 44:329-358, 1981.

21. E.R. Jessup, I.C.F. Ipsen. Improving the accuracy of inverse iteration. SIAM Journal on Scientific and Statistical Computing, 13(2):550—72, 1992.

22. Cullum Jane K., Willoughby Ralph A. Lanczos algorithms for large symmetric eigenvalue computations. URL www.netXib.org, source code.

23. C. Lanczos. An iterative method for the solution of the eigenvalue problem of linear differential and integral operators. J. Res. Nat. Bur. Standards, Sect. D, 45:255-282, 1950.

24. Osni Marques. Eigensolvers and applications in finite element analyses. Advanced Solution Procedures on Innovative Computer Architectures (M. Papadrakakis and G. Dugeda, Eds.), 6G 79, 1996.

25. A.M. Matsekh. The Godunov-inverse iteration algorithm for symmetric tridiagonal matrices. Bulletin of the Novosibirsk Computing Center. Series: Computer Science, (19) :39—50, 2003.

26. A.M. Matsekh. The Godunov-inverse iteration: a fast and accurate solution to the symmetric tridiagonal eigenvalue problem. Applied Numerical Mathematics, 54(2):208-221, 2005.

27. A.M. Matsekh. Numerical comparison of symmetric eigenvalue and SVD solvers basedon the two-sided Sturm sequences (Sturm library) and the corresponding LAPACK functionality. Technical report, Los Alamos National Laboratory, '2005. LA-UR-05-8-101.

28. Gerard Meurant, Zdenek Strakos. The lanczos and conjugate gradient algorithms in finite precision arithmetic. Acta Numenca, 15:471-542, 2000.

29. National Institute of Standards and Technology. Matrix Market. URL http://math. nist.gov/MatrixMarket/, a visual repository of test data for numerical linear algebra.

30. J.C. Nedclec. Mixed finite elements in IR„ Numer. Math., 35:315-3-11, 1980.

31. J.M. Ortega, H.F. Kaiser. The Z/' and qr methods for symmetric tridiagonal matrices. Numer. Math., 5:211-225, 1963.

32. M. Overton. Numerical Computing with IEEE Floating Point, Arithmetic. SIAM, Philadelphia, PA, USA, 2001.

33. G3. C.C. Paige. Computational variants of the Lanczos method for the eigenproblem. Journal of the Institute of Mathematics and Its Applications, 10:373-381, 1972.

34. G4. C.C. Paige. Error analysis of the Lanczos algorithm for tridiagonalizing a symmetric matrix. Journal of the Institute of Mathematics and Its Applications, 18:3-11-349, 1976.

35. G5. Bcresford N. Parlett. The Symmetric Eigenvalue Problem. Prentice-Hall Inc., Englewood Cliffs, N.J, 1980.

36. G6. B.N. Parlett, I.S. Dhillon. Fernando's solution to Wilkinson's problem: An application of double factorization. Lin. Alg. Appl., 267:247-270, December 1997.

37. G7| B.N. Parlett, O.A. Marques. An implementation of the dqds algorithm (positive case).1.near Algebra and its Applications, 309:217-259, 2000.

38. Y. Saad. Numerical Methods for Large Eigenvalue Problems. Manchester University Press, Manchester, UK, 1992. URL http://www-users.cs.umn.edu/~saad/books. html.

39. Yousef Saad. SPAR.SKIT. URL http://www-users.cs.umn.edu/~saad/software/ SPARSKIT/sparskit.html, a basic tool-kit for sparse matrix computations.

40. B.T. Smith, J.M. Boyle, J.J. Dongarra, B.S. Garbow, Y. Ikebe, V.C. Klema, C.I3. Moler. Matrix Eigensystem Routines EISPACK Guide, volume 6 of Lecture Notes in Computer Science. Springer-Verlag, Berlin, 1976.

41. Danny C. Sorensen. Numerical methods for large eigenvalue problems. Acta Numerica, 11(00) :519 584, 2002.

42. GAV. Stewart. Introduction to m,atrix computalions. Academic Press, New York, 1973. G.W. Stewart. On the early history of the singular value decomposition. SI AM Review, 35(4):551-56G, 1993.

43. G.W. Stewart. Matrix Algorithms, volume II, Eigensystems. G.W. Stewart, 1999. URL http://www.cs.umd.edu/~stewart.

44. J.II. Wilkinson. The Algebraic Eigenvalue Problem. Oxford University Press, 1965. J.H. Wilkinson. Rounding Errors in Algebraic Processes. Dover Publications, reprint edition, 1994.

45. Д.К. Фаддеев, В.II. Фаддеева. Вычислительные методы липейной алгебры,, 3-е стереотипное издание. Лань, Санкт-Петербург, 2002.

46. C.К. Годунов, Прокопов Г.П. Применение метода минимальных итераций для вычисления собственных значений эллиптических операторов. Журнал вычислительной математики и математической физики, 10(о): 1180 -1190, 1970.

47. С.К. Годунов. Совремеиные аспекты линейной алгебры. Новосибирск, Научная книга, 1997.

48. С.К. Годунов. Лекции по современным аспектам линейной алгебры. Новосибирск, Научная книга, 2002. Университетская серия, том 12.

49. С.К. Годунов, А.Г. Антонов, О.П. Киршпок, В.И. Костин. Гарантированная точность решения систем линейных уравнений в евклидовых пространствах. Наука, Новосибирск, 1988.

50. С.К. Годунов, В.И. Костин, А.Д. Митченко. Вычисление собственного вектора симметричной трехдиагональной матрицы. Институт Математики СОАН СССР, Новосибирск, 1983. Препринт 44.

51. B.II. Кублановская. О некоторых алгоритмах для решения полной проблемы собственных значений. Журнал вычислительной математики и математической физики, 1(4):555 570, 1961.

52. В.II. Ильин. Методы неполной факторизации для решения алгебраических систем. Москва, Иаука-Физматлит, 1995.

53. В.П. Ильин. Методы, конечных разностей и конечных объем,ов для эллиптических уравнений. Новосибирск, Издательство Института Математики, 2000.

54. A.M. Мацех, Э.П. Шурина. Нахождение спектрального разложения симметричных матриц и сингулярного разложения несимметричных матриц с оцениваемой томностью. Автометрия, 43(2):81-96, 2007.

55. O.B. Нечаев, Э.П. Шурина. Многоссточный алгоритм решения векторным методом конечных элементов трехмерного уравнения Гельмгольца. Математическое моделирование, 17:92-102, 2005.