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

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

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

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

Дана трёхдиагональная матрица [math]A[/math]. Нужно найти её разложение в произведение двух двудиагональных [math]L[/math] (левая с единичной диагональю) и [math]U[/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]n - 1[/math] чисел [math]l_k[/math] (поддиагональ матрицы [math]L[/math]), одномерный массив [math]n[/math] чисел [math]u_k[/math] (диагональ матрицы [math]U[/math]), одномерный массив [math]n - 1[/math] чисел [math]r_k[/math] (наддиагональ матрицы [math]U[/math]).

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

[math] \begin{align} u_1 & = a_1, \\ r_k & = c_k, & k & = 1, \dots, n - 1, \\ l_k & = \frac{b_k}{u_k}, & k & = 1, \dots, n - 1, \\ u_k & = a_k - l_{k-1} r_{k-1}, & 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 = 7[/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]c[/math] (элементы [math]c_i[/math] с номерами от [math]1[/math] до [math]n - 1[/math]), массив [math]b[/math] (элементы [math]b_i[/math] с номерами от [math]1[/math] до [math]n - 1[/math]).

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

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

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

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

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

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

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

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

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

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

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