Linpack benchmark

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

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

Linpack benchmark является тестом, разработанным как приложение к библиотеке Linpack для тестирования производительности компьютерных систем. По сути это генерирование некоторой системы линейных алгебраических уравнений (СЛАУ) [math]A \vec{x} = \vec{b}[/math] со случайной плотной квадратной матрицей [math]A[/math] и вектором [math]\vec{b}[/math] таким, чтобы решением был вектор [math]\vec{x}[/math], у которого все элементы равны 1, с последующим решением этой СЛАУ. Решение в силу неспецифичности (с рядом оговорок) матрицы делается в два этапа. Сначала для матрицы системы выполняется [math]PLU[/math]-разложение ([math]P[/math] — матрица перестановок, [math]L[/math] — левая треугольная с единичными диагональными элементами, [math]U[/math] — правая треугольная матрица) методом Гаусса с выбором ведущего элемента по столбцу. На втором этапе с помощью полученного разложения последовательно решается СЛАУ вида [math]P L U \vec{x} = \vec{b}[/math]. На третьем этапе вычисляется невязка [math]\vec{r} = A \vec{x} - \vec{b}[/math] и её характеристики, после чего выводятся данные о полученных точности и производительности вычислений. Интересно, что производительность вычисляется только для основной части алгоритма, т. е. туда не включены ни вычисление невязки, ни вычисление норм.

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

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

Исходные данные: невырожденная (ГСЧ подобран так, чтобы это обеспечивалось, в версии от 1992 г.) квадратная матрица [math]A[/math] (элементы [math]a_{ij}[/math]), вектор правой части [math]\vec{b}[/math] (элементы [math]b_i[/math]). В реальности они не вполне входные (элементы матрицы генерирует псевдослучайно сама программа, все они находятся в диапазоне (-1/2, +1/2), [math]b_i = \sum\limits_{j = 1}^i a_{ij}[/math]), но для разбора общей схемы алгоритма удобно считать их входными данными.

Вычисляемые данные: левая треугольная матрица [math]L[/math] (элементы [math]l_{ij}[/math]), правая треугольная матрица [math]U[/math] (элементы [math]u_{ij}[/math]), матрица перестановок [math]P[/math] (вычисляется не в развёрнутом виде, а как набор номеров ведущих элементов), вектор решения [math]x[/math] (элементы [math]x_i[/math]).

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

Первая часть — разложение матрицы в произведение двух треугольных и матрицы перестановок: [math]A = P L U[/math] — выполняется методом Гаусса с выбором ведущего элемента по столбцу. При этом матрица перестановок хранится с помощью вектора, характеризующего выполненные перестановки. Матрица [math]L[/math] левая треугольная с единичными диагональными элементами, матрица [math]U[/math] — правая треугольная.

Вторая часть состоит из двух решений треугольных СЛАУ. Сначала c помощью «прямой-обратной» подстановки решается система [math]L \vec{y} = P T \vec{b}[/math], затем с помощью обратной подстановки — система [math]U \vec{x} = \vec{y}[/math]. При этом в тесте треугольные матрицы хранятся в соответствующих местах матрицы [math]A[/math]. Поэтому перед началом третьей части её генерируют заново (с теми же стартовыми параметрами для ГСЧ).

Третья часть (для которой не измеряется скорость работы) состоит в вычислении сначала невязки решения [math]\vec{r} = A \vec{x} - \vec{b}[/math] и затем её равномерной нормы, а также нормы решения. По окончании этой части алгоритм выдаёт данные по точности полученного решения, а также вычисляет производительность системы.

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

Основная часть алгоритма — разложение матрицы методом Гаусса с выбором ведущих элементов по столбцу — определяет основную часть операций. Следующая по значимости часть — вычисление невязки и её нормы — реализовано подпрограммой вычисления суммы вектора и произведения матрицы на другой вектор и функцией вычисления равномерной нормы. Третьи по значимости части — обратный ход метода Гаусса и «прямой-обратный» ход метода Гаусса.

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

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

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

Математическое описание алгоритма описано выше. При этом в стандартном тесте из-за того, что разложение Гаусса вычисляется на месте исходной матрицы, а также из-за того, что при прямом-обратном и обратном ходе решение заанимает место правой части, после вычисления решения приходится заново перевычислять как матрицу системы, так и правую часть — для последующего вычисления невязки. После вычисления невязки её норма в принципе может быть вычислена в любом порядке с нормой решения, это не влияет на алгоритм и его результаты.

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

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

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

Граф состоит из информационных графов своих основных частей. Первая часть (метод Гаусса с выбором ведущего элемента по столбцу) имеет самый сложный из графов. На первом рисунке отражён упрощённый граф этой части, post factum по результатам выборов ведущих элементов. Красным отображены операции определения максимумов, с запоминанием позиций. Голубым цветом выделено обращение очередного ведущего элемента, зелёным – умножения в текущем столбце, синим – операции [math]a b + c[/math] в других столбцах. Приведём также граф одного из слоёв (одного шага метода). Обозначения – те же самые, что в полном графе. Пунктир означает рассылку данных из ведущей строки.

Вторая часть состоит из почти аналогичных фрагментов. Первый из них – «прямой ход», а именно решение СЛАУ с левой треугольной матрицей [math]L \vec{y} = \vec{b}[/math]. Диагональные элементы матрицы единичны, поэтому в этом куске нет операций деления. Все зелёные вершины отвечают операциям [math]a b + c[/math]. Второй – «обратный ход», а именно решение СЛАУ с правой (верхней) треугольной матрицей [math]U \vec{x} = \vec{y}[/math]. Для её матрицы нет информации о специфичности, и потому в графе появляются также операции деления , они обозначены желтым цветом.

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

При изучении ресурса параллелизма теста видно, что основная сложность (квадратичная) будет у разложения Гаусса с выбором ведущего элемента по столбцу (именно из-за наличия выбора элемента по столбцу). Остальные части алгоритма имеют не более линейной сложности. Поэтому для оптимизации вычислений возможному исследователю нужно заняться именно оптимизацией разложения Гаусса с выбором ведущего элемента по столбцу. Кроме этого, первый кусок второй части (решение первой из треугольных СЛАУ) может быть пристыкован к первой части исполнен почти полностью «на её фоне», ибо направление перевычисления данных у них совпадают. Однако второй кусок второй части, несмотря на его схожесть с первым, может быть вычислен не ранее, чем будет выполнен весь первый кусок, поскольку у них направления перевычисления данных противоположны.

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

По сути теста у него всегда одни и те же данные, причём вычисляемые в самом же тесте. Поэтому в строгом смысле у алгоритма нет входных данных.

Выходными данными являются вычисляемые и выдаваемые тестом нормы невязки, решения, матрицы и их отношения, а также вычисляемая мощность вычислений целевой вычислительной системы.

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

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

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

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

2.2.2 Описание локальности реализации алгоритма

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

Рисунок 12.1. Тест Linpack. Общий профиль обращений в память

На рисунке 12.1 представлен общий профиль обращений в память для теста Linpack. В данном тесте задействовано 3 массива. Обращения к первым двум выделены на рис. 12.1 зеленым цветом (фрагменты 1 и 2), остальные обращения выполняются к элементам третьего массива – матрицы, содержащей коэффициенты СЛАУ. Данный тест состоит из двух этапов – факторизации матрицы и собственно решения СЛАУ. Разделение между этапами отмечено на рисунке оранжевой линией. Видно, что первый этап обладает итерационной структурой, при этом на каждой следующей итерации отбрасываются из рассмотрения несколько элементов матрицы с наименьшими индексами. Также можно заметить, что на втором этапе выполняется два прохода по элементам матрицы с возрастающим шагом, сначала последовательный, затем обратный. Однако в целом профиль устроен достаточно сложно, и для понимания свойств локальности необходимо детальное рассмотрение.

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

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

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

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

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

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

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

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

Рисунок 12.5. Одна итерация фрагмент 3, небольшая часть

Видно, что первые два этапа отличаются от остальных – на них выполняются некоторые начальные действия. Затем начинаются подобные итерации внутренних циклов, причем в обоих случаях (выделенные зеленым части 1 и 2) выполняется последовательный перебор элементов. Основная разница заключается лишь в том, что в первом цикле обращения к каждому элементу выполняются дважды.

Обращения к данному массиву в рамках внутренних циклов характеризуются высокой пространственной локальностью. При этом во втором цикле используются одни и те же элементы, что приводит также и к достаточно высокой временной локальности. Таким образом, в рамках одной итерации (рис. 4) локальность обращений высока. А при условии, что на разных итерациях происходят обращения к одним и тем же элементам (и число обращений на каждой итерации не так уж велико), это приводит к дальнейшему повышению локальности.

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

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

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

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

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

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

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