Обратная подстановка (вещественный вариант)

Обратная подстановка для неособенной верхней треугольной матрицы
Последовательный алгоритм
Последовательная сложность	[math]O(n^2)[/math]
Объём входных данных	[math]O(n^2)[/math]
Объём выходных данных	[math]n[/math]
Параллельный алгоритм
Высота ярусно-параллельной формы	[math]O(n)[/math]
Ширина ярусно-параллельной формы	[math]O(n)[/math]

Основные авторы описания: А.В.Фролов, Вад.В.Воеводин (раздел 2.2), А.М.Теплов (раздел 2.4).

1 Свойства и структура алгоритма

1.1 Общее описание алгоритма

Обратная подстановка - решение системы линейных алгебраических уравнений (СЛАУ) [math]Ux = y[/math] с верхней треугольной матрицей [math]U[/math]. Матрица [math]U[/math] может быть одной из составляющих матрицы [math]A[/math] в каких-либо разложениях и получается либо из [math]LU[/math]-разложения последней каким-либо из многочисленных способов (например, простое разложение Гаусса, разложение Гаусса с выбором ведущего элемента, компактная схема Гаусса, разложение Холецкого и др.), либо из других (например из QR-разложения). В силу треугольности [math]U[/math] решение СЛАУ является одной из модификаций общего метода подстановки и записывается простыми формулами.

В методом обратной подстановки назван также и метод решения СЛАУ с нижней треугольной матрицей. Там же отмечено, что в литературе иногда под обратной подстановкой имеют в виду, как и здесь, только решения СЛАУ с верхней треугольной матрицей, а решение нижних треугольных систем называют прямой подстановкой. Такой же системы названий будем придерживаться и здесь, во избежание одноимённого названия различных алгоритмов. Кроме того, обратная подстановка, представленная здесь, одновременно может быть частью метода Гаусса для решения СЛАУ, а именно - его обратным ходом, чего нельзя сказать про прямую подстановку.

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

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

Исходные данные: верхняя треугольная матрица [math]U[/math] (элементы [math]u_{ij}[/math]), вектор правой части [math]y[/math] (элементы [math]y_{i}[/math]).

Вычисляемые данные: вектор решения [math]x[/math] (элементы [math]x_{i}[/math]).

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

[math] \begin{align} x_{n} & = y_{n}/u_{nn} \\ x_{i} & = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}, \quad i \in [1, n - 1]. \end{align} [/math]

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

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

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

[math] \sum_{j = i+1}^{n} u_{ij} x_{j} [/math]

в режиме накопления или, в зависимости от требований задачи, без его использования, с их последующим вычитанием из компоненты вектора [math]y[/math] и деления на диагональный элемент матрицы [math]U[/math]. В отечественных реализациях, даже в последовательных версиях, упомянутый способ представления не используется. Дело в том, что даже в этих реализациях метода вычисление сумм вида

[math] y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} [/math]

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

[math] y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} [/math]

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

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

Как уже упоминалось в описании ядра алгоритма, основную часть метода обратной подстановки составляют множественные (всего их [math]n-1[/math]) вычисления сумм

[math]y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} [/math]

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

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

Последовательность действий описываемого алгоритма:

1. [math]x_{n} = y_{n}/u_{nn}[/math]

Далее для всех [math]i[/math] от [math]n-1[/math] до [math]1[/math] по убыванию выполняются

2. [math]x_{i} = \left (y_{i} - \sum_{j = i+1}^{n} u_{ij} x_{j} \right ) / u_{ii}[/math].

Особо отметим, что вычисления сумм вида [math]y_{i} - \sum{}_{j = i+1}^{n} u_{ij} x_{j}[/math] производят в режиме накопления вычитанием из [math]y_{i}[/math] произведений [math]u_{ij} x_{j}[/math] для [math]j[/math] от [math]n[/math] до [math]i + 1[/math], c убыванием [math]j[/math]. Использование другого порядка выполнения суммирования приводит к резкому ухудшению параллельных свойств алгоритма, хотя, к сожалению, это кое-где встречается в литературе и пакетах программ. В качестве примера использования другого порядка вычислений можно привести фрагмент программы из, где обратная подстановка является обратным ходом в методе Гаусса, а возрастание индекса суммирования связано, в основном, с ограничениями используемого авторами книги старого диалекта Фортрана.

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

Для алгоритма обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка [math]n[/math] в последовательном (наиболее быстром) варианте требуется:

• [math]n[/math] делений,
• [math]\frac{n^2-n}{2}[/math] сложений (вычитаний),
• [math]\frac{n^2-n}{2}[/math] умножений.

Выполнение умножений и сложений (вычитаний) составляет основную часть алгоритма.

При этом использование режима накопления требует выполнения умножений и вычитаний в режиме двойной точности (или использования функции аналогичной DPROD на языке Фортран).

Таким образом, при классификации по последовательной сложности, метод обратной подстановки относится к алгоритмам со сложностью [math]O(n^2)[/math].

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

Опишем граф алгоритма как аналитически, так и в виде рисунка.

Граф алгоритма обратной подстановки состоит из двух групп вершин, расположенных в целочисленных узлах двух областей разной размерности.

Рис. 1. Обратная подстановка

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

Делимое в этой операции:

• при [math]i = n[/math] — элемент входных данных, а именно [math]y_{n}[/math];
• при [math]i \lt n[/math] — результат срабатывания операции, соответствующей вершине из второй группы, с координатами [math]i[/math], [math]i+1[/math].

Делитель для этой операции - элемент входных данных, а именно [math]u_{nn}[/math].

Результат выполнения операции является выходным данным [math]x_{i}[/math].

Вторая группа вершин расположена в двумерной области, соответствующая ей операция [math]a-bc[/math]. Естественно введённые координаты области:

• [math]i[/math] — меняется в диапазоне от [math]n-1[/math] до [math]1[/math], принимая все целочисленные значения;
• [math]j[/math] — меняется в диапазоне от [math]n[/math] до [math]i+1[/math], принимая все целочисленные значения.

Аргументы операции:

• [math]a[/math]:
  • при [math]j = n[/math] элемент входных данных [math]y_{i}[/math];
  • при [math]j \lt n[/math] — результат выполнения операции, соответствующей вершине из второй группы, с координатами [math]i, j+1[/math];
• [math]b[/math] — элемент входных данных, а именно [math]u_{ij}[/math];
• [math]c[/math] — результат срабатывания операции, соответствующей вершине из первой группы, с координатой [math]j[/math].

Результат срабатывания операции является промежуточным данным алгоритма.

Описанный граф представлен на рис. 1, выполненном для случая [math]n = 5[/math]. Здесь вершины первой группы обозначены жёлтым цветом и знаком деления, вершины второй — зелёным цветом и буквой f. Изображена только подача входных данных из вектора [math]y[/math], подача элементов матрицы [math]U[/math], идущая во все вершины, на рисунке не представлена.

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

Для обратной подстановки в случае решения линейной системы с верхней треугольной матрицей порядка [math]n[/math] в параллельном варианте требуется последовательно выполнить следующие ярусы:

• [math]n[/math] ярусов делений (в каждом из ярусов одно деление),
• по [math]n - 1[/math] ярусов умножений и сложений/вычитаний (в каждом из ярусов — линейное количество операций, от [math]1[/math] до [math]n-1[/math].

Таким образом, в параллельном варианте, в отличие от последовательного, вычисления делений будут определять довольно значительную долю требуемого времени. При реализации на конкретных архитектурах наличие в отдельных ярусах ЯПФ отдельных делений может породить и другие проблемы. Например, при реализации метода обратной подстановки на ПЛИСах остальные вычисления (умножения и сложения/вычитания) могут быть конвейеризованы, что даёт экономию и по ресурсам на программируемых платах; деления из-за их изолированности приведут к занятию ресурсов на платах, которые будут простаивать большую часть времени.

При этом использование режима накопления требует совершения умножений и вычитаний в режиме двойной точности, а в параллельном варианте это означает, что практически все промежуточные вычисления для выполнения алгоритма в режиме накопления должны быть двойной точности. В отличие от последовательного варианта это означает некоторое увеличение требуемой памяти.

Таким образом, при классификации по высоте ЯПФ, алгоритм обратной подстановки относится к алгоритмам с линейной сложностью. При классификации по ширине ЯПФ его сложность также будет линейной.

1.9 Входные и выходные данные алгоритма

Входные данные: верхняя треугольная матрица [math]U[/math] (элементы [math]u_{ij}[/math]), вектор правой части [math]y[/math] (элементы [math]y_{i}[/math]).

Объём входных данных: [math]\frac{n (n + 3)}{2}[/math] (в силу треугольности достаточно хранить только ненулевые элементы матрицы [math]U[/math]).

Выходные данные: вектор решения [math]x[/math] (элементы [math]x_{i}[/math]).

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

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

Соотношение последовательной и параллельной сложности в случае неограниченных ресурсов, как хорошо видно, является линейным (отношение квадратической к линейной).

При этом вычислительная мощность алгоритма обратной подстановки, как отношение числа операций к суммарному объему входных и выходных данных – всего лишь константа.

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

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

У алгоритма обратной подстановки существует несколько блочных вариантов. Граф некоторых из них совпадает с графом точечного варианта, различия связаны в основном с порядком прохождения основных циклов алгоритма, а именно, с их развёртыванием и перестановкой. Эти приёмы могут помочь в оптимизации обменов на конкретных вычислительных системах.

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

2.1 Особенности реализации последовательного алгоритма

В простейшем варианте метод обратной подстановки на Фортране можно записать так:

        X(N) = Y(N) / U(N,N)
        DO  I = N-1, 1, -1
            S = Y(I)
            DO  J = N, I+1, -1
                S = S - DPROD(U(I,J), X(J))
            END DO
            X(I) = S / U(I,I)
	END DO

При этом для реализации режима накопления переменная [math]S[/math] должна быть двойной точности.

2.2 Локальность данных и вычислений

2.2.1 Локальность реализации алгоритма

2.2.1.1 Структура обращений в память и качественная оценка локальности

Рисунок 2. Обратный ход метода Гаусса решения СЛАУ. Общий профиль обращений в память

На рис.2 представлен профиль обращений в память для обратного хода метода Гаусса решения СЛАУ. Хорошо видно, что профиль состоит из двух этапов, идущих друг за другом (граница между ними выделена оранжевой линией). Это соответствует двум циклам, из которых собственно состоит исходный код обратного хода метода Гаусса. Из самого рисунка это заметить достаточно сложно, но анализ исходного кода показывает, что профиль образуется из обращений к 4 массивам. Три из них выделены на рис.2 зеленым; к четвертому относятся все остальные обращения. Сразу отметим, что в данном профиле общее число обращений всего в несколько раз больше числа задействованных элементов массивов (4500 против 1000), а в таких условиях сложно добиться высокой локальности.

Чтобы выяснить это, перейдем к подробному рассмотрению каждого из массивов в отдельности.

На рис.3 представлен фрагмент 1. Оранжевой линей также проведено разделение 2 этапов. Видно, что второй этап устроен очень просто и представляет собой последовательный перебор в обратном порядке. Первый этап имеет итерационный характер, причем на каждой итерации отбрасывается элемент массива с наименьшим номером. Подобный профиль характеризуется высокой пространственной локальностью, поскольку элементы перебираются подряд; достаточно высокой временной локальностью, поскольку на каждой итерации происходит повторное обращение к тем же элементам.

Рисунок 3. Фрагмент 1 (профиль обращений к первому массиву)

На рис.4 представлен фрагмент 2, показывающий обращения ко второму массиву. Сразу заметим, что данный фрагмент еще меньше предыдущего – здесь задействовано всего 600 обращений в память. По сути, данный фрагмент является аналогичным фрагменту 1 (рис.3), с единственной разницей – здесь итерации расположены в обратном порядке. Однако это изменение не оказывает особого влияния ни на пространственную, ни на временную локальность, поэтому данный фрагмент обладает теми же свойствами.

Рисунок 4. Фрагмент 2 (профиль обращений ко второму массиву)

Далее рассмотрим фрагмент 3 (рис.5). Здесь также оранжевой линией проведено разделение между двумя этапами, и также на второй этап приходится совсем немного обращений. Обращения на первом этапе образуют аналог случайного доступа, достаточно часто встречающийся, например, в случае косвенной адресации. При этом в некоторый момент к определенному элементу происходит достаточно много обращений подряд, после чего этот элемент более не используется. Такой профиль обычно характеризуется низкой пространственной и временной локальностью, что, однако, в данном случае в достаточной степени нивелируется малым числом задействованных элементов.

Рисунок 5. Фрагмент 3 (профиль обращений к третьему массиву)

Далее перейдем к фрагменту, занимающему основную часть рис.2. Данный профиль отображен на рис.6. В первую очередь стоит отметить особенность данного профиля – здесь задействовано более 1000 элементов, в то время как в остальных профилях – порядка 30. При этом число обращений даже меньше, чем обращений к первому или третьему массиву.

Отметим также, что здесь большая часть обращений сосредоточена на втором этапе. Первый же этап является подобием первого этапа предыдущего массива (рис.5), за тем лишь исключением, что в данном профиле отсутствуют множественные обращения подряд к одному элементу перед тем, как элемент перестанет использоваться. С учетом того, что первый этап состоит всего из порядка 500 обращений, достаточно равномерно распределенных на отрезке в 1000 элементов, это говорит об очень низкой как пространственной, так и временной локальности.

Другой характер обращений можно наблюдать на втором этапе. Видно, что он обладает более высокой пространственной локальностью, так как обращения явно сгруппированы в кластеры, при этом кластеры обладают подобным строением. Однако структура кластера видно плохо, поэтому требуется дальнейшее приближение.

Рисунок 6. Фрагмент 4 (профиль обращений к четвертому массиву)

На рис.7 показаны два кластера, выделенные на рис.6 зеленым цветом. Такой масштаб позволяет сразу увидеть, что каждый кластер представляет собой последовательный перебор с небольшим шагом определенного набора элементов массива.

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

Рисунок 7. Два кластера из фрагмента 4

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

2.2.1.2 Количественная оценка локальности

Основной фрагмент реализации, на основе которого были получены количественные оценки, приведен здесь (функция Kernel2). Условия запуска описаны здесь.

Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.

На рис.8 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Можно увидеть, что, в соответствии с высказанным в предыдущем разделе предположением о негативном влиянии одного из массивов, данная программа показывает достаточно низкую производительность, заметно ниже, чем у прямого хода метода Гаусса.

Рисунок 8. Сравнение значений оценки daps

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

На рис.9 приведены значения cvg для того же набора реализаций, отсортированные по убыванию (чем меньше cvg, тем в общем случае выше локальность). Можно увидеть, что, согласно данной оценке, профиль обращений обладает низкой локальностью, лишь немногим лучше профиля программы со случайным доступом в память. Это повторяет выводы, сделанные на основе оценки daps.

Рисунок 9. Сравнение значений оценки cvg

2.3 Возможные способы и особенности параллельной реализации алгоритма

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

        DO PARALLEL I = 1, N
           X(I) = Y(I)
        END DO
        DO J = N, 1, -1
           X(J) = X(J) / U(J,J)
           DO PARALLEL I = 1, J-1
              X(I) = X(I) - U(I,J)*X(J)
           END DO
        END DO

так и с использованием "скошенного параллелизма" в главном гнезде циклов.

2.4 Масштабируемость алгоритма и его реализации

2.4.2 Масштабируемость реализации алгоритма

Проведём исследование масштабируемости параллельной реализации алгоритма согласно методике. Исследование проводилось на суперкомпьютере "Ломоносов" Суперкомпьютерного комплекса Московского университета. Набор и границы значений изменяемых параметров запуска реализации алгоритма:

• число процессоров [8, 16. 24, 32, 48, 64, 96, 128, 192, 256];
• размер матрицы [1024 : 4096] с шагом 1024.

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

• минимальная эффективность реализации 1.152e-15%;
• максимальная эффективность реализации 5.39e-10%.

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

Рисунок 10. Параллельная реализация алгоритма. Изменение производительности в зависимости от числа процессоров и размера матрицы.

Рисунок 11. Параллельная реализация алгоритма. Изменение эффективности в зависимости от числа процессоров и размера матрицы.

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

Исследованная реализация алгоритма

2.5 Динамические характеристики и эффективность реализации алгоритма

Для проведения экспериментов использовалась реализация алгоритма обратной подстановки. Все результаты получены на суперкомпьютере "Ломоносов". Использовались процессоры Intel Xeon X5570 с пиковой производительностью в 94 Гфлопс, а также компилятор intel 13.1.0. На рисунках показана эффективность реализации алгоритма обратной подстановки на 32 процессах.

Рисунок 9. График загрузки CPU при выполнении алгоритма обратной подстановки

На графике загрузки процессора видно, что почти все время работы программы уровень загрузки составляет около 50%. Это указывает на равномерную загруженность вычислениями процессоров, при использовании 8 процессов на вычислительный узел и без использования Hyper Threading.

Рисунок 10. График операций с плавающей точкой в секунду при выполнении алгоритма обратной подстановки

На Рисунке 10 показан график количества операций с плавающей точкой в секунду. На графике видна общая низкая производительность вычислений около 150 Мфлопс в пике и около 75 Мфлопс в среднем по всем узлам.

Рисунок 11. График кэш-промахов L1 в секунду при работе алгоритма обратной подстановки

На графике кэш-промахов первого уровня видно, что число промахов не очень большое для нескольких ядер и находится на уровне 3,5 млн/сек в среднем (в пике до 8,5 млн/сек), что указывает на не очень интенсивные вычисления в большей части процессов. Это указывает на неравномерное распределение вычислений и не очень высокую производительность приложения.

Рисунок 12. График кэш-промахов L3 в секунду при работе алгоритма обратной подстановки

На графике кэш-промахов третьего уровня видно, что число промахов очень немного и находится на уровне 1,2 млн/сек в среднем, однако в пике по всем узлам значения доходят до 3 млн/сек. Соотношение кэш промахов L1|L3 для процессов с доходит до 3, однако в среднем около 2,5. Это указывает на очень плохую локальность вычислений как у части процессов, так и для всех в среднем, и это является признаком низкой производительности.

Рисунок 13. График количества чтений из оперативной памяти при работе алгоритма обратной подстановки

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

Рисунок 14. График количества записей в оперативную память при работе алгоритма обратной подстановки

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

Рисунок 15. График скорости передачи по сети Infiniband в байт/сек при работе алгоритма обратной подстановки

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

Рисунок 16. График скорости передачи по сети Infiniband в пакетах/сек при работе алгоритма обратной подстановки

На графике скорости передачи данных в пакетах в секунду наблюдается очень высокая интенсивность передачи данных выраженная в пакетах в секунду. Это говорит о том, что, вероятно, процессы обмениваются не очень существенными объемами данных, но очень интенсивно. Так же наблюдается почти меньший дизбаланс между процессами чем наблюдаемый в графиках использования памяти и вычислений и передачи данных в байтах/сек. Это указывает на то, что процессы обмениваются по алгоритму одинаковым числом пакетов, однако получают разные объемы данных и ведут неравномерные вычисления. Скорее всего активность Infiniband является одной из основных причин снижения производительности

Рисунок 17. График числа процессов, ожидающих вхождения в стадию счета (Loadavg), при работе алгоритма обратной подстановки

На графике числа процессов, ожидающих вхождения в стадию счета (Loadavg), видно, что на протяжении всей работы программы значение этого параметра постоянно и приблизительно равняется 8. Это свидетельствует о стабильной работе программы с загруженными вычислениями всеми узлами. Это указывает на очень рациональную и статичную загрузку аппаратных ресурсов процессами. И показывает достаточно хорошую эффективность выполняемой реализации. В целом, по данным системного мониторинга работы программы можно сделать вывод о том, что программа работала неэффективно, но стабильно. Использование памяти не очень интенсивное, а использование коммуникационной среды крайне интенсивное, при этом объемы передаваемых данных очень низкие. Это указывает на требовательность к латентности коммуникационной среды алгоритмической части программы. Низкая эффективность связана судя по всему с достаточно высоким объемом пересылок на каждом процессе, интенсивными обменами сообщениями.

2.6 Выводы для классов архитектур

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

2.7 Существующие реализации алгоритма

Вещественный вариант обратной подстановки реализован как в основных библиотеках программ отечественных организаций, так и в западных пакетах LINPACK, LAPACK, ScaLAPACK и др. При этом в отечественных реализациях, как правило, выполнены стандартные требования к методу с точки зрения ошибок округления, то есть, реализован режим накопления, и обычно нет лишних операций. Реализация точечного варианта алгоритма в современных западных пакетах обычно происходит из одной и той же реализации метода в LINPACK, использующей библиотеку BLAS.

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