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

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]c_k[/math] (поддиагональ матрицы [math]L[/math]), одномерный массив [math]n[/math] чисел [math]l_k[/math] (диагональ матрицы [math]L[/math]).

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

[math] \begin{align} l_1 & = \sqrt{a_1}, \\ c_k & = \frac{b_k}{l_k}, \quad k = 1, \dots, n - 1, \\ l_k & = a_k - c^2_{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]. Количество извлечений квадратного корня равно [math]n[/math]. Поэтому алгоритм должен быть отнесён к алгоритмам линейной сложности по количеству последовательных операций.

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

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

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

Последовательный вариант вычисления не имеет ресурсов параллелизма. Его ЯПФ единственна и совпадает с последовательным алгоритмом. Таким образом, в получившемся алгоритме высота параллельной формы будет равна [math]n - 1[/math] операций умножения плюс [math]n - 1[/math] операций вычитания плюс [math]n - 1[/math] операций деления, а также [math]n[/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[/math]), массив [math]c[/math] (элементы [math]c_k[/math] с номерами от [math]1[/math] до [math]n - 1[/math]).

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

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

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

Обращает на себя внимание то, что для трёхдиагональных матриц метод Холецкого лучше компактной схемы Гаусса только в отношении требуемой памяти (за счёт симметричности), но проигрывает по скорости из-за «лишних» извлечений квадратного корня.

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

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

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

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

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