Вычислительные методы алгебры и оценивания. И.В. Семушкин (2011) (1185350), страница 13
Текст из файла (страница 13)
Еще одной сложностью является внедрениепроцедуры выбора главного элемента в эти алгоритмы.754 Алгоритмы окаймления в LU-разложении4.3Окаймление неизвестной части разложенияОсновная работа в алгоритмах окаймления приходится на решение треугольных систем (4.3). Это матрично-векторные операции, которые можнореализовать в виде подпрограмм по типу рис.
4.1, добиваясь в них максимальной для данной машины эффективности. Еще один способ организации вычислений, который называют алгоритмом Донгарры–Айзенштата,имеет то преимущество, что его основной операцией является матричновекторное умножение. Математически алгоритм можно описать следующимобразом. Выпишем из равенства (4.1) три других соотношения, на этот раздля ajj , a3j и aTj3 .
Отсюда получимujj = ajj − lTj1 u1j ,uTj3 = aTj3 − lTj1U13,l3j = (a3j − L31u1j ) /ujj .(4.4)Характер доступа к данным при таком вычислении показан на рис. 4.3.U11U13L̄11⇐ uTj3L31A33⇑ l3jРис. 4.3. Доступ к данным в алгоритмах окаймления неизвестной части разложения.L̄11 , U11 — вычисление закончено, обращений больше нет. L31 , U13 — вычислениезакончено, но обращения происходят. A33 — обращений не было. Вычисляются: j-йстолбец l3j матрицы L̄ и j-я строка uTj3 матрицы UВидно, что блоки U13 и L31, необходимые для вычислений (4.4), на каждый такой момент времени уже известны, так же как и все другие величиныв правых частях равенств (4.4), поэтому здесь нет решения уравнений.Основной операцией в (4.4) является умножение вектора на прямоугольную матрицу.
Можно реализовать такие умножения посредством скалярныхпроизведений или линейных комбинаций, что приводит к двум различнымформам алгоритма, показанным на рис. 4.4.764.3 Окаймление неизвестной части разложенияДля j = 1 до nДля k = 1 до j − 1Для i = j до naji = aji − ljk akiДля k = 1 до j − 1Для i = j + 1 до naij = aij − lik akjДля s = j + 1 до nlsj = asj /ajjДля j = 1 до nДля i = j + 1 до nДля k = 1 до j − 1aij = aij − lik akjДля i = j до nДля k = 1 до j − 1aji = aji − ljk akiДля s = j + 1 до nlsj = asj /ajjРис. 4.4. Алгоритмы Донгарры–Айзенштата окаймления неизвестной части L̄Uразложения: алгоритм линейных комбинаций (слева) и алгоритм скалярных произведений (справа)Первый цикл по k,i на рис.
4.4 (слева) производит последовательные модификации j-й строки матрицы A, которая по окончании циклаk превращается в j-ю строку матрицы U . Эти модификации можно рассматривать как вычитание векторно-матричного произведения lTj1U13 изaTj3 во второй формуле uTj3 = aTj3 − lTj1 U13 в (4.4) c помощью линейных комбинаций строк U . В случае j = i результатом модификациибудет первая величина ujj в (4.4). Во второй паре циклов по k и iвыполняются модификации j-го столбца матрицы A по формуле l3j == (a3j − L31u1j ), т.
е. по второй формуле (4.4) с точностью до деленияна ujj c вычислением матрично-векторного произведения L31u1j в (4.4) посредством линейных операций столбцов матрицы L. Обратите внимание, —в отличие от алгоритма отложенных модификаций на рис. 3.6, теперь длинывекторов, участвующих в линейных комбинациях, одинаковы. Второй оператор цикла с индексом k можно удалить, причем программа снова будетверна; мы вставили этот оператор, чтобы подчеркнуть наличие линейныхкомбинаций.
Отметим еще, что в первом цикле по k,i происходит обращениек строкам матрицы A, а во втором цикле по k,i — к ее столбцам. Следовательно, эта форма неэффективна для векторных компьютеров, требующихразмещения элементов вектора в смежных ячейках памяти.Алгоритм на рис. 4.4 (справа) использует скалярные произведения. Наj-м шаге первый цикл по i, k вычисляет, с точностью до финального деления, j-й столбец матрицы L; с этой целью j-й столбец в A модифицируетсяпосредством скалярных произведений строк L и j-го столбца U . Во втором774 Алгоритмы окаймления в LU-разложениицикле по i, k вычисляется j-я строка матрицы U , для чего j-я строка в Aмодифицируется посредством скалярных произведений j-й строки L и столбцов U . Снова требуется доступ к строкам и столбцам матрицы A.
Потенциальное преимущество алгоритма Донгарры–Айзенштата заключается в том,что в некоторых векторных компьютерах матрично-векторные умножениявыполняются весьма эффективно. Вычисления по этим формам показаны втабл. 4.2.Таблица 4.2. Вычисления по алгоритмам на рис. 4.4 для примера 3.3. Позицииэлемента–делителя столбца L̄ показаны выделенным шрифтомj=12 4 −41/2423/2812/2506115j=224 −41/2243/22/212/21/206−215j=324 −41/2243/22/232/21/22/36−2−65j=424 −41/2243/22/232/21/22/36−2−64Как и в алгоритмах окаймления известной части разложения (подразд. 4.2), в данных алгоритмах дополнительную сложность представляетвнедрение процедуры выбора главного элемента.4.4Задание на лабораторный проект № 3Написать и отладить программу, реализующую ваш вариант методаисключения с выбором главного элемента, для численного решения системлинейных алгебраических уравнений Ax = f , вычисления det A и A−1.Предусмотреть сообщения, предупреждающие о невозможности решенияуказанных задач с заданной матрицей A.Отделить следующие основные части программы:а) подпрограмму факторизации матрицы A, отвечающую вашему варианту метода исключения;б) подпрограмму решения систем линейных алгебраических уравнений;в) подпрограмму вычисления определителя матриц;г) подпрограммы обращения матриц;д) сервисные подпрограммы.784.4 Задание на лабораторный проект № 3Уделить особое внимание эффективности программы (в смысле экономии оперативной памяти).
Предусмотреть пошаговое выполнение алгоритмаисключения с выводом результата на экран.Выполнить следующие пункты задания.1. Провести подсчет фактического числа выполняемых операций умножения и деления при решении системы линейных алгебраических уравнений,сравнить его с оценочным числом (n3/3).2. Определить скорость решения задач (решение систем линейных алгебраических уравнений, обращение матриц) с учетом времени, затрачиваемогона разложение матрицы.
Для этого спроектировать и провести эксперимент,который охватывает матрицы порядка от 5 до 100 (через 5 порядков). Представить результаты в виде таблицы и графика зависимости времени выполнения (в минутах и секундах) от порядка матриц. Таблицу и график вывестина экран.3. Оценить точность решения систем линейных алгебраических уравнений, имеющих тот же самый порядок, что и задачи из пункта 2. Для этогосгенерировать случайные матрицы A, задать точное решение x∗ и образоватьправые части f = Ax∗. Провести анализ точности вычисленного решения xот порядка матрицы. Результаты представить в виде таблицы и графика.Для заполнения матрицы A использовать случайные числа из диапазонаот −100 до 100.
В качестве точного решения взять вектор x∗ = (1, 2, . . . , n)T ,где n — порядок матрицы. Для оценки точности использовать норму вектораkxk∞ = max(|xi|).i(4.5)4. Повторить пункт 3 задания для плохо обусловленных матриц (см. подразд. 2.6), имеющих порядок от 4 до 40 с шагом 4.5. Вычислить матрицу A−1 следующими двумя способами.Способ 1 — через решение системы AX = I, где I — единичная матрица.Способ 2 — через разложение матрицы A в произведение элементарныхматриц, обращение которых осуществляется аналитически, а их произведение дает матрицу A−1.Сравнить затраты машинного времени и точность обращения матриц прииспользовании указанных способов 1 и 2.
Эксперименты провести для случайных матриц порядков от 5 до 100 через 5. Для оценки точности в обоихспособах использовать оценочную формулу−1−1−1kA−1т − Aпр k ≤ kI − AAпр k · kAk .(4.6)794 Алгоритмы окаймления в LU-разложенииНорму матрицы следут вычислять в соответствии со следующим определением:!nX|aij | ,(4.7)kAk∞ = maxij=1−1где A−1т — точное значение обратной матрицы, а Aпр — приближенное значение, полученное в результате обращения каждым из способов 1 и 2.6.
Провести подсчет фактического числа выполняемых операций умножения и деления при обращении матриц первым и вторым способами, сравнитьего с оценочным числом (n3).Замечание 4.1. По ходу проведения численных экспериментов наэкран дисплея должны выводиться следующие таблицы.Решение систем линейных алгебраических уравнений:ПорядокВремяТочностьТеоретическоеРеальноечисло операций число операцийАналогичная таблица должна быть построена для плохо обусловленныхматриц.Обращение матриц:ПорядокВремяТочностьЧисло операцийспос. 1 спос. 2 спос. 1 спос.