Разложение трёхдиагональной матрицы, вещественная версия, последовательный вариант без корней (первая стадия прогонки)

1 Описание свойств и структуры алгоритма

1.1 Словесное описание алгоритма

1.1.1 Решаемая задача

Дана трёхдиагональная симметричная матрица [math]A[/math]. Нужно найти её разложение в произведение двудиагональной [math]L[/math] (левая с единичной диагональю), диагональной [math]D[/math] и [math]L^T[/math]. Для задачи решения СЛАУ с матрицей [math]A[/math] это разложение является первой стадией решения (известного под именем прогонки).

1.1.2 Общая схема

Реализуется модификация компактной схемы Гаусса, с учётом структуры матриц.

1.2 Математическое описание

Исходные данные: одномерный массив [math]n - 1[/math] чисел [math]b_k[/math] (поддиагональ исходной матрицы, равная наддиагонали), одномерный массив [math]n[/math] чисел [math]a_k[/math] (диагональ исходной матрицы).

Вычисляемые данные: одномерный массив [math]n - 1[/math] чисел [math]l_k[/math] (поддиагональ матрицы [math]L[/math]), одномерный массив [math]n[/math] чисел [math]d_k[/math] (диагональ матрицы [math]D[/math]).

Формулы метода:

[math] \begin{align} d_1 & = a_1, l_k & = \frac{b_k}{d_k}, \quad k = 1, \dots, n - 1, d_k & = a_k - l_{k - 1} b_{k - 1}, \quad k = 2, \dots, n. \end{align} [/math]

1.3 Вычислительное ядро алгоритма

Вычислительное ядро последовательной схемы решения задачи можно представить в качестве последовательного набора [math]n - 1[/math] «тройных» операций (деление элементов одного из входных массивов на элементы одного из получаемых массивов, умножение элементов получаемого массива на элемент одного из входных массивов и вычитание результата из следующего элемента другого входного массива).

1.4 Макроструктура алгоритма

Как уже записано в описании ядра алгоритма, основную часть вычисления составляет массовый последовательный набор «тройных» операций (деление элементов одного из входных массивов на элементы одного из получаемых массивов, умножение элементов получаемого массива на элемент одного из входных массивов и вычитание результата из следующего элемента другого входного массива).

1.5 Описание схемы реализации последовательного алгоритма

Формулы метода описаны выше. Последовательность исполнения — по возрастанию индексов.

1.6 Последовательная сложность алгоритма

Для вычисления разложения матрицы порядка [math]n[/math], количество операций умножения равно количеству операций вычитания и количеству операций деления и равно [math]n - 1[/math]. Поэтому алгоритм должен быть отнесён к алгоритмам линейной сложности по количеству последовательных операций.

1.7 Информационный граф

Опишем граф алгоритма в виде рисунка. Изображён граф для [math]n = 10[/math]. Как видно, граф чисто последовательный. Вершине соответствует «тройная» операция (деление-умножение-вычитание).

1.8 Описание ресурса параллелизма алгоритма

Последовательный вариант вычисления не имеет ресурсов параллелизма. Его ЯПФ единственна и совпадает с последовательным алгоритмом. Таким образом, в получившемся алгоритме высота параллельной формы будет равна [math]n - 1[/math] операций умножения плюс [math]n - 1[/math] операций вычитания плюс [math]n - 1[/math] операций деления. В таком виде алгоритм должен быть отнесён к алгоритмам линейной сложности по высоте параллельной формы. Ширина параллельной формы равна [math]1[/math], что даёт нам постоянную сложность по ширине параллельной формы.

1.9 Описание входных и выходных данных

Входные данные: массив [math]a[/math] (элементы [math]a_i[/math] с номерами от [math]1[/math] до [math]n[/math]), массив [math]b[/math] (элементы [math]b_i[/math] с номерами от [math]1[/math] до [math]n - 1[/math]).

Дополнительные ограничения: отсутствуют.

Объём входных данных: [math]2n - 1[/math].

Выходные данные: массив [math]l[/math] (элементы [math]l_k[/math] с номерами от [math]1[/math] до [math]n - 1[/math]), массив [math]d[/math] (элементы [math]d_k[/math] с номерами от [math]1[/math] до [math]n[/math]).

Объём выходных данных: [math]2n - 1[/math].

1.10 Свойства алгоритма

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является константой (1). При этом вычислительная мощность алгоритма, как отношение числа операций к суммарному объему входных и выходных данных — всего-навсего 1 (входных и выходных данных почти столько же, сколько операций; если точнее — даже больше на 33 процента). При этом алгоритм полностью детерминирован. Дуги информационного графа локальны.

2 Программная реализация

На Фортране можно записать так:

	D(1) = A(1)
	DO I = 1, N-1
		L(I) = B(I)/D(I)
		D(I+1) = A(I+1)-L(I)*B(I)
	END DO

Если хранить выходные данные на местах входных, что часто практикуется, то программа будет выглядеть так:

	DO I = 1, N-1
		R=B(I)
		B(I) = R/A(I)
		A(I+1) = A(I+1)-R*B(I)
	END DO