Algorithm level

Difference between revisions of "Linpack benchmark"

From Algowiki
Jump to navigation Jump to search
[unchecked revision][checked revision]
 
(45 intermediate revisions by 2 users not shown)
Line 1: Line 1:
 
{{level-a}}
 
{{level-a}}
  
Primary authors of this description: [[:ru:Участник:Frolov|A.V.Frolov]], [[:ru:Участник:VadimVV|Vad.V.Voevodin]] ([[#Locality of data and computations|Section 2.2]]), [[:ru:Участник:Teplov|A.M.Teplov]] (раздел [[#Scalability of the algorithm and its implementations|Section 2.4]])
+
Primary authors of this description: [[:ru:Участник:Frolov|A.V.Frolov]].
  
 
== Properties and structure of the algorithm ==
 
== Properties and structure of the algorithm ==
Line 7: Line 7:
 
=== General description of the algorithm ===
 
=== General description of the algorithm ===
  
'''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''' is a test developed as a supplement to the Linpack library for testing the performance of computer systems. Actually, this test generates a system of linear algebraic equations (SLAE) <math>A \vec{x} = \vec{b}</math> with a random dense square matrix <math>A</math> and a vector <math>\vec{b}</math> chosen so that the solution of the system is the vector <math>\vec{x}</math> all of whose components are ones; then this SLAE is solved. In view of the general (with certain provisos) nature of the matrix, the solution of the system is divided into several stages. First, the <math>PLU</math>-decomposition of <math>A</math> is calculated by using Gaussian elimination with column pivoting (<math>P</math> is a permutation matrix, <math>L</math> is a lower triangular matrix with unit diagonal, and <math>U</math> is an upper triangular matrix). At the second stage, the calculated decomposition is used to solve the SLAE  <math>P L U \vec{x} = \vec{b}</math>. At the third stage, the residual <math>\vec{r} = A \vec{x} - \vec{b}</math> and the related characteristics are calculated. Then the data concerning the resulting accuracy and performance are produced. The performance is only estimated for the main part of the algorithm; that is, the calculations spent on the residual and norms are ignored.
  
Не следует путать Linpack benchmark, разработанный для тестирования компьютеров традиционной архитектуры, и High Performance Linpack benchmark, который разработан для тестирования компьютерных систем с параллельной архитектурой и в основе  которого лежит хоть и сходный, но другой алгоритм.
+
The Linpack benchmark, which was developed for testing computers of traditional architecture, should not be confused with High Performance Linpack benchmark. The latter was devised for testing parallel computer systems and is based on a different (albeit similar) algorithm.
  
=== Mathematical description of the algorithm ===
+
=== Mathematical description of the algorithm ===  
 +
Input data: square matrix <math>A</math> (with entries <math>a_{ij}</math>) (in 1992 version, the random number generator is chosen so that to ensure the non-singularity of <math>A</math>), the right-hand side vector <math>\vec{b}</math> (with components <math>b_i</math>). Actually, they are not exactly input data (because the entries of the matrix are pseudo-random numbers generated by the program itself; these numbers are contained in the interval (-1/2, +1/2), and the right-hand sides are calculated as <math>b_i = \sum\limits_{j = 1}^i a_{ij}</math>). However, in order to analyze the general scheme of the algorithm, it is convenient to regard <math>A</math> and <math>\vec{b}</math> as the input data.
  
Исходные данные: невырожденная (ГСЧ подобран так, чтобы это обеспечивалось, в версии от 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>), но для разбора общей схемы алгоритма удобно считать их входными данными.
+
Data to be calculated: lower triangular matrix <math>L</math> (with entries <math>l_{ij}</math>), upper triangular matrix <math>U</math> (with entries <math>u_{ij}</math>), permutation matrix <math>P</math> (which is calculated as a collection of pivot indices rather than a full square array), and the solution vector <math>x</math> (with components <math>x_i</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>).
+
Here, we do not give detailed formulas of individual parts of the algorithm. These parts are independent algorithms, and their detailed descriptions should be presented separately. Consequently, we only discuss the general scheme of the algorithm.  
  
Детальные формулы отдельных частей алгоритма здесь описывать не будем, поскольку они являются самостоятельными алгоритмами и их детальные описания должны быть выполнены раздельно. Поэтому здесь будет приведена только общая схема алгоритма.
+
The first part — decomposing the matrix into the product of two triangular matrices and a permutation matrix: <math>A = P L U</math> — is performed by using Gaussian elimination with column pivoting. The permutation matrix is stored as a vector containing information on the performed permutations. The matrix <math>L</math> is lower triangular with unit diagonal, while <math>U</math> is an upper triangular matrix.  
  
Первая часть — разложение матрицы в произведение двух треугольных и матрицы перестановок: <math>A = P L U</math> — выполняется методом Гаусса с выбором ведущего элемента по столбцу. При этом матрица перестановок хранится с помощью вектора, характеризующего выполненные перестановки. Матрица <math>L</math> левая треугольная с единичными диагональными элементами, матрица <math>U</math> — правая треугольная.  
+
The second part is solving two triangular SLAEs. First, the system <math>L \vec{y} = P T \vec{b}</math> is solved by forward substitution, then back substitution is used to solve the system <math>U \vec{x} = \vec{y}</math>. This test stores both triangular matrices in the corresponding portions of the matrix <math>A</math>. For this reason, <math>A</math> is generated again (with the same starting parameters for the random number generator) before passing to the third part.  
  
Вторая часть состоит из двух решений треугольных СЛАУ. Сначала c помощью «прямой-обратной» подстановки решается система <math>L \vec{y} = P T \vec{b}</math>, затем с помощью обратной подстановки — система <math>U \vec{x} = \vec{y}</math>. При этом в тесте треугольные матрицы хранятся в соответствующих местах матрицы <math>A</math>. Поэтому перед началом третьей части её генерируют заново (с теми же стартовыми параметрами для ГСЧ).
+
The third part (for which the performance is not measured) calculates first the residual of the solution <math>\vec{r} = A \vec{x} - \vec{b}</math> and then its uniform norm, as well as the norm of the solution. On the completion of this stage, the algorithm outputs data related to the accuracy of the calculated solution and evaluates the performance of the system.
 
 
Третья часть (для которой не измеряется скорость работы) состоит в вычислении сначала невязки решения <math>\vec{r} = A \vec{x} - \vec{b}</math> и затем её равномерной нормы, а также нормы решения. По окончании этой части алгоритм выдаёт данные по точности полученного решения, а также вычисляет производительность системы.
 
  
 
=== Computational kernel of the algorithm ===
 
=== Computational kernel of the algorithm ===
  
Основная часть алгоритма — разложение матрицы методом Гаусса с выбором ведущих элементов по столбцу — определяет основную часть операций. Следующая по значимости часть — вычисление невязки и её нормы — реализовано подпрограммой вычисления суммы вектора и произведения матрицы на другой вектор и функцией вычисления равномерной нормы. Третьи по значимости части — обратный ход метода Гаусса и «прямой-обратный» ход метода Гаусса.
+
The main part of the algorithm, which is the decomposition of a matrix via Gaussian elimination with column pivoting, is responsible for the bulk of operations. The second most important part is the calculation of the residual and its norm. It is implemented by the subroutine that calculates the sum of a vector and a matrix-vector product and by the function calculating the uniform norm. Less important parts are forward substitution in the Gauss method and back substitution in this method.
  
 
=== Macro structure of the algorithm ===
 
=== Macro structure of the algorithm ===
  
Как уже записано в описании ядра алгоритма, разложение матрицы методом Гаусса с выбором ведущих элементов по столбцу — определяет основную часть операций и при этом является первой макрооперацией. После неё идёт «прямой-обратный» ход метода Гаусса и сразу за ним — обратный ход метода Гаусса. После повторного заполнения матрицы следующей макрооперацией является вычисление невязки решения СЛАУ, а затем — вычисление её и вектора решения равномерных норм.
+
As already said in the description of the computational kernel of the algorithm, the decomposition of a matrix via Gaussian elimination with column pivoting is responsible for the bulk of operations. This decomposition is the first macro operation. It is followed by the forward substitution and, right away, by the back substitution. On the regeneration of the matrix, the next macro operation is calculating the residual for the solution of the SLAE, followed by the calculation the uniform norms of the residual and the solution.
  
 
=== Implementation scheme of the serial algorithm ===
 
=== Implementation scheme of the serial algorithm ===
  
Математическое описание алгоритма описано выше. При этом в стандартном тесте из-за того, что разложение Гаусса вычисляется на месте исходной матрицы, а также из-за того, что при прямом-обратном и обратном ходе решение заанимает место правой части, после вычисления решения приходится заново перевычислять как матрицу системы, так и правую часть — для последующего вычисления невязки. После вычисления невязки её норма в принципе может быть вычислена в любом порядке с нормой решения, это не влияет на алгоритм и его результаты.
+
The mathematical description of the algorithm was given above. The standard test calculates the triangular decomposition in-place, and the solutions to the triangular systems in the forward and back substitution are written to the locations of the right-hand sides. For this reason, on the calculation of the solution, the coefficient matrix and the right-hand side are regenerated for the subsequent calculation of the residual. As soon as the residual is calculated, its norm and the norm of the solution can be calculated in any order, which does not affect the algorithm and its results.
  
 
=== Serial complexity of the algorithm ===
 
=== Serial complexity of the algorithm ===
  
Поскольку основной частью алгоритма является разложение Гаусса с выбором ведущего элемента по столбцу, то он, как и эта часть, относится к алгоритмам с ''кубической'' сложностью. Остальные части квадратичны по сложности и потому не дают сравнимого вклада в сложность алгоритма.
+
Since the main part of the algorithm is the decomposition of a matrix via Gaussian elimination with column pivoting, the algorithm, as well as this decomposition, has ''cubical'' complexity. Other parts have quadratic complexity; consequently, they do not give a comparable contribution to the overall complexity of the algorithm.
  
 
=== Information graph ===
 
=== Information graph ===
  
Граф состоит из информационных графов своих основных частей. Первая часть (метод Гаусса с выбором ведущего элемента по столбцу) имеет самый сложный из графов. На рис.1 отражён упрощённый граф этой части, post factum по результатам выборов ведущих элементов. Зеленым отображены операции определения максимумов, с запоминанием позиций. Желтым цветом выделено обращение очередного ведущего элемента, коричневым – умножения в текущем столбце, фиолетовым – операции <math>a b + c</math> в других столбцах. На рис.2 приведён граф одного из слоёв (одного шага метода). Обозначения – те же самые, что в полном графе. Пунктир означает рассылку данных из ведущей строки.
+
The graph consists of the information graphs of its basic parts. The first part (Gaussian elimination with column pivoting) has the most complex graph. Figure 1 shows the simplified graph of this part, post factum the choice of pivots. The operations of finding the maxima in the columns and saving the corresponding positions are given in green. The inversion of the current pivot, multiplications in the current column, and operations of the type <math>a b + c</math> in the other columns are depicted in yellow, brown, and purple, respectively. Figure 2 shows the graph of a layer (one step of the method). The notation is the same as in the full graph. Dotted lines denote data communication from the pivot row.
  
 
{| align="center"
 
{| align="center"
 
     |- valign="top"
 
     |- valign="top"
     | [[File:Gaussian_elimination.png|700px|thumb|Рисунок 1. Упрощённый граф метода Гаусса с выбором ведущего элемента по столбцу]]
+
     | [[File:Gaussian_elimination.png|700px|thumb|Figure 1. Simplified graph of Gaussian elimination with column pivoting]]
     | [[File:Gaussian_elimination_step.png|500px|Рисунок 2. Граф одного шага метода Гаусса|thumb]]
+
     | [[File:Gaussian_elimination_step.png|500px|Figure 2. Graph of one step of Gaussian elimination|thumb]]
 
|}
 
|}
  
Вторая часть состоит из почти аналогичных фрагментов. Первый из них – «прямой ход», а именно решение СЛАУ с левой треугольной матрицей <math>L \vec{y} = \vec{b}</math> (рис.3). Диагональные элементы матрицы единичны, поэтому в этом куске нет операций деления. Все зелёные вершины отвечают операциям <math>a b + c</math>. Второй – «обратный ход», а именно решение СЛАУ с правой (верхней) треугольной матрицей <math>U \vec{x} = \vec{y}</math> (рис.4). Для её матрицы нет информации о специфичности, и потому в графе появляются также операции деления , они обозначены желтым цветом.
+
The second part consists of almost identical fragments. The first fragment describes the forward substitution, namely, the solution of the SLAE with a lower triangular matrix <math>L \vec{y} = \vec{b}</math> (Fig.3). The diagonal entries of this matrix are ones; consequently, this fragment contains no divisions. All the green nodes correspond to the operations of the type <math>a b + c</math>. The second fragment describes the back substitution, namely, the solution of the SLAE with an upper triangular matrix <math>U \vec{x} = \vec{y}</math> (Fig.4). This matrix has no specific features; hence, the associated graph contains division operations, which are shown in yellow.
  
 
{| align="center"
 
{| align="center"
 
     |- valign="top"
 
     |- valign="top"
     | [[File:DirectL.png|500px|thumb|Рисунок 3. Граф решения СЛАУ с левой треугольной матрицей]]
+
     | [[File:DirectL.png|500px|thumb|Figure 3. Graph of solving the SLAE with a lower triangular matrix]]
     | [[File:DirectU.png|450px|thumb|Рисунок 4. Граф решения СЛАУ с правой треугольной матрицей]]
+
     | [[File:DirectU.png|450px|thumb|Figure 4. Graph of solving the SLAE with an upper triangular matrix]]
 
|}
 
|}
  
 
=== Parallelization resource of the algorithm ===
 
=== Parallelization resource of the algorithm ===
  
При изучении ресурса параллелизма теста видно, что основная сложность (''квадратичная'') будет у разложения Гаусса с выбором ведущего элемента по столбцу (именно из-за наличия выбора элемента по столбцу). Остальные части алгоритма имеют не более линейной сложности. Поэтому для оптимизации вычислений возможному исследователю нужно заняться именно оптимизацией разложения Гаусса с выбором ведущего элемента по столбцу. Кроме этого, первый кусок второй части (решение первой из треугольных СЛАУ) может быть пристыкован к первой части исполнен почти полностью «на её фоне», ибо направление перевычисления данных у них совпадают. Однако второй кусок второй части, несмотря на его схожесть с первым, может быть вычислен не ранее, чем будет выполнен весь первый кусок, поскольку у них направления перевычисления данных противоположны.
+
The analysis of parallelization resource of the algorithm shows that Gaussian elimination with column pivoting has the highest, namely, ''quadratic'' complexity (exactly because of the pivoting). The other parts of the algorithm have at most linear complexity. Therefore, to optimize calculations, the virtual researcher must optimize Gaussian elimination with column pivoting. Besides, the first fragment of the second part (that is, the solution of the first triangular SLAE) can be attached to the first part and executed almost entirely in the background of the latter. Indeed, both have the same direction of recalculating the data. However, the rest of the second part, despite its similarity to the first fragment, can only be executed on completion the latter. The reason is that these two fragments have opposite directions of recalculating the data.
  
 
=== Input and output data of the algorithm ===
 
=== Input and output data of the algorithm ===
  
По сути теста у него всегда одни и те же данные, причём вычисляемые в самом же тесте. Поэтому в строгом смысле у алгоритма нет входных данных.  
+
Actually, this test has always the same data, which are calculated in the test itself. Consequently, in the strict sense, it has no input data.  
  
Выходными данными являются вычисляемые и выдаваемые тестом нормы невязки, решения, матрицы и их отношения, а также вычисляемая мощность вычислений целевой вычислительной системы.
+
The output data, produced by the test, are the norms of the residual, solution, and matrix, the ratios of these norms, and the computational power of the target computer system.
  
 
=== Properties of the algorithm ===
 
=== Properties of the algorithm ===
Line 75: Line 74:
 
=== Implementation peculiarities of the serial algorithm ===
 
=== Implementation peculiarities of the serial algorithm ===
  
Строго говоря, при реализации теста испытывающий компьютерную систему испытывает и версии Линпака или другого пакета, откуда он должен вставлять вызовы разложения Гаусса и решения СЛАУ с помощью полученного разложения. Поэтому находящиеся в сети версии теста не являются догмой, поскольку с их помощью можно тестировать и другие версии реализации метода Гаусса. Как уже отмечено выше, из-за того, что разложение Гаусса вычисляется на месте исходной матрицы, а также из-за того, что при прямом-обратном и обратном ходе решение занимает место правой части, после вычисления решения приходится заново перевычислять как матрицу системы, так и правую часть — для последующего вычисления невязки. Это решение вполне можно поменять с целью оптимизации теста.
+
Strictly speaking, the investigator implementing the test for a computer system tests at the same time LINPAC or any other package from which procedures for triangular decomposition and solving a SLAE with the help of this decomposition are called. The web versions of the test can also be used for testing other implementations of the Gauss method. It was already noted above that the triangular decomposition of the original matrix is calculated in-place, and, in the course of forward and back substitution, the solutions are written to the locations of right-hand sides. This causes the necessity to regenerate both the coefficient matrix and the right-hand side for the subsequent calculation of the residual. This organization of the algorithm may well be changed in order to optimize the test.
 
 
=== Locality of data and computations ===
 
==== Locality of implementation ====
 
===== Structure of memory access and a qualitative estimation of locality =====
 
 
 
[[file:linpack 1.PNG|thumb|center|700px|Рисунок 5. Тест Linpack. Общий профиль обращений в память]]
 
 
 
На рис.5 представлен общий профиль обращений в память для теста Linpack. В данном тесте задействовано 3 массива. Обращения к первым двум выделены на рис.5 зеленым цветом (фрагменты 1 и 2), остальные обращения выполняются к элементам третьего массива – матрицы, содержащей коэффициенты СЛАУ. Данный тест состоит из двух этапов – факторизации матрицы и собственно решения СЛАУ. Разделение между этапами отмечено на рисунке оранжевой линией. Видно, что первый этап обладает итерационной структурой, при этом на каждой следующей итерации отбрасываются из рассмотрения несколько элементов матрицы с наименьшими индексами.  Также можно заметить, что на втором этапе выполняется два прохода по элементам матрицы с возрастающим шагом, сначала последовательный, затем обратный. Однако в целом профиль устроен достаточно сложно, и для понимания свойств локальности необходимо детальное рассмотрение.
 
 
 
Фрагмент 1 показан отдельно на рис.6. При подобном приближении хорошо видно, что данный фрагмент состоит из двух последовательных переборов всех элементов массива, при этом размер фрагмент очень мал по сравнению с общим профилем. Такой последовательный перебор характеризуется высокой пространственной и достаточно низкой временной локальностью, поскольку к каждому элементу обращение выполняется всего дважды.
 
 
 
[[file:linpack 2.PNG|thumb|center|400px|Рисунок 6. Фрагмент 1 (профиль обращений к первому массиву)]]
 
 
 
Далее перейдем к рассмотрению фрагмента 2, представленного на рис.7. Данный профиль состоит из двух похожих этапов, состоящих из набора итераций. На каждой итерации выполняется последовательный перебор элементов массива, при этом на следующей итерации отбрасывается их рассмотрения один из элементов: на первом этапе – элемент с минимальным индексом, на втором – с максимальным индексом.
 
 
 
Подобный профиль обладает очень высокой как пространственной, так и временной локальностью, поскольку почти всегда обращения выполняются к соседним элементам, причем элементы достаточно часто используются повторно. Однако заметим, что данный фрагмент состоит всего из около 2000 элементов, что значительно меньше общего числа обращений в программе. Значит, большая часть (и, скорее всего,  основное влияние) приходится на последний массив.
 
 
 
[[file:linpack 3.PNG|thumb|center|400px|Рисунок 7. Фрагмент 2 (профиль обращений к второму массиву)]]
 
 
 
На рис.8 представлен фрагмент 3, выделенный на рис.5, соответствующий одной итерации цикла, в котором выполняются обращения к данному массиву. Из данного рисунка хорошо видно, что итерация состоит из двух частей, соответствующих двум внутренним циклам. В обоих случаях число обращений на любой итерации внутреннего цикла примерно равно, однако в первом внутреннем цикле на следующей итерации перебирается следующая часть элементов массива, в то время как во втором цикле перебираются одни и те же элементы.
 
 
 
[[file:linpack 4.PNG|thumb|center|700px|Рисунок 8. Фрагмент 3 (профиль обращений к третьему массиву, одна итерация)]]
 
 
 
Остается выяснить структуру обращений в память на итерациях данных двух циклов. Рассмотрим более подробно часть фрагмента, выделенную на рис.8 зеленым (рис.9).
 
 
 
[[file:linpack 5.PNG|thumb|center|700px|Рисунок 9. Одна итерация фрагмент 3, небольшая часть]]
 
 
 
Видно, что первые два этапа отличаются от остальных – на них выполняются некоторые начальные действия. Затем начинаются подобные итерации внутренних циклов, причем в обоих случаях (выделенные зеленым части 1 и 2) выполняется последовательный перебор элементов. Основная разница заключается лишь в том, что в первом цикле обращения к каждому элементу выполняются дважды.
 
 
 
Обращения к данному массиву в рамках внутренних циклов характеризуются высокой пространственной локальностью. При этом во втором цикле используются одни и те же элементы, что приводит также и к достаточно высокой временной локальности. Таким образом, в рамках одной итерации (рис. 4) локальность обращений высока. А при условии, что на разных итерациях происходят обращения к одним и тем же элементам (и число обращений на каждой итерации не так уж велико), это приводит к дальнейшему повышению локальности.
 
 
 
===== Quantitative estimation of locality =====
 
 
 
Основной фрагмент реализации, на основе которого были получены количественные оценки, приведен [http://git.algowiki-project.org/Voevodin/locality/blob/master/benchmarks/linpack/linpack_np.cpp здесь] (оригинал исходного кода взят [http://people.sc.fsu.edu/~%20jburkardt/cpp_src/linpack_bench/linpack_bench_s.cpp отсюда]). Условия запуска описаны [http://git.algowiki-project.org/Voevodin/locality/blob/master/README.md здесь].
 
 
 
Первая оценка выполняется на основе характеристики daps, которая оценивает число выполненных обращений (чтений и записей) в память в секунду. Данная характеристика является аналогом оценки flops применительно к работе с памятью и является в большей степени оценкой производительности взаимодействия с памятью, чем оценкой локальности. Однако она служит хорошим источником информации, в том числе для сравнения с результатами по следующей характеристике cvg.
 
 
 
На рис.10 приведены значения daps для реализаций распространенных алгоритмов, отсортированные по возрастанию (чем больше daps, тем в общем случае выше производительность). Можно увидеть, что данный показатель у теста Linpack один из самых высоких, что говорит об эффективном взаимодействии с памятью.
 
 
 
[[file:linpack daps ru.PNG|thumb|center|700px|Рисунок 10. Сравнение значений оценки daps]]
 
 
 
Вторая характеристика – cvg – предназначена для получения более машинно-независимой оценки локальности. Она определяет, насколько часто в программе необходимо подтягивать данные в кэш-память. Соответственно, чем меньше значение cvg, тем реже это нужно делать, тем лучше локальность.
 
 
 
На рис.11 приведены значения cvg для того же набора реализаций, отсортированные по убыванию (чем меньше cvg, тем в общем случае выше локальность). Можно увидеть, что, согласно данной оценке, локальность обращений в память теста Linpack очень высока, что соответствует полученным ранее оценкам и выводам.
 
 
 
[[file:linpack cvg.PNG|thumb|center|700px|Рисунок 11. Сравнение значений оценки cvg]]
 
  
 
=== Possible methods and considerations for parallel implementation of the algorithm ===
 
=== Possible methods and considerations for parallel implementation of the algorithm ===
=== Scalability of the algorithm and its implementations ===
+
=== Run results ===
==== Scalability of the algorithm ====
 
==== Scalability of of the algorithm implementation ====
 
Проведём исследование масштабируемости реализации HPL теста Linpack согласно [[Scalability methodology|методике]]. Исследование проводилось на суперкомпьютере "Ломоносов" [http://parallel.ru/cluster Суперкомпьютерного комплекса Московского университета].
 
 
 
Набор и границы значений изменяемых [[Глоссарий#Параметры запуска|параметров запуска]] реализации алгоритма:
 
 
 
* число процессоров [8 : 128] с шагом 8;
 
* размер матрицы [1000 : 100000] c шагом 1000.
 
 
 
В результате проведённых экспериментов был получен следующий диапазон [[Глоссарий#Эффективность реализации|эффективности реализации]] алгоритма:
 
 
 
* минимальная эффективность реализации 0.039%;
 
* максимальная эффективность реализации 86,7%.
 
 
 
На следующих рисунках приведены графики [[Глоссарий#Производительность|производительности]] и эффективности теста HPL в зависимости от изменяемых параметров запуска.
 
 
 
[[file:Масштабируемость Linpack performance.png|thumb|center|700px|Рисунок 12. Тест HPL. Изменение производительности в зависимости от числа процессоров и размера матрицы.]]
 
[[file:Linpack масштабируемость эффективность.png|thumb|center|700px|Рисунок 13. Тест HPL. Изменение производительности в зависимости от числа процессоров и размера матрицы.]]
 
 
 
Построим оценки масштабируемости теста Linpack:
 
* По числу процессов: -0.061256. При увеличении числа процессов эффективность на рассмотренной области изменений параметров запуска уменьшается, однако в целом уменьшение не очень быстрое. Малая интенсивность изменения объясняется высокой общей эффективностью работы при больших размерах задачи, однако при росте числа процессов эффективность равномерно уменьшается. Это объясняется также тем, что с ростом [[Глоссарий#Вычислительная сложность|вычислительной сложности]] падение эффективности становится не таким быстрым. Уменьшение эффективности на рассмотренной области работы параллельной программы объясняется быстрым ростом накладных расходов на организацию параллельного выполнения. С ростом вычислительной сложности задачи эффективность снижается так же быстро, но при больших значениях числа процессов. Это подтверждает предположение о том, что накладные расходы начинают сильно превалировать над вычислениями.
 
* По размеру задачи: 0.010134. При увеличении размера задачи эффективность возрастает. Эффективность возрастает тем быстрее, чем большее число процессов используется для выполнения. Это подтверждает предположение о том, что размер задачи сильно влияет на эффективность выполнения приложения. Оценка показывает, что с ростом размера задачи эффективность на рассмотренной области значений параметров запуска сильно увеличивается. Также, учитывая разницу максимальной и минимальной эффективности почти в 80%, можно сделать вывод, что рост эффективности при увеличении размера задачи наблюдается на большей части рассмотренной области значений.
 
* По двум направлениям: 0.0000284 При рассмотрении увеличения как вычислительной сложности, так и числа процессов на всей рассмотренной области значений эффективность увеличивается, однако скорость увеличения эффективности не очень большая.
 
 
 
[http://git.algowiki-project.org/Teplov/Scalability/tree/master/HPL Исследованная реализация теста]
 
 
 
=== Dynamic characteristics and efficiency of the algorithm implementation ===
 
 
=== Conclusions for different classes of computer architecture ===
 
=== Conclusions for different classes of computer architecture ===
=== Existing implementations of the algorithm ===
 
  
 
== References ==
 
== References ==
 +
 
<references />
 
<references />
  
[[Category:Started articles]]
+
[[Category:Articles in progress]]
  
 
[[Ru:Linpack benchmark]]
 
[[Ru:Linpack benchmark]]

Latest revision as of 16:15, 12 July 2022


Primary authors of this description: A.V.Frolov.

1 Properties and structure of the algorithm

1.1 General description of the algorithm

Linpack benchmark is a test developed as a supplement to the Linpack library for testing the performance of computer systems. Actually, this test generates a system of linear algebraic equations (SLAE) [math]A \vec{x} = \vec{b}[/math] with a random dense square matrix [math]A[/math] and a vector [math]\vec{b}[/math] chosen so that the solution of the system is the vector [math]\vec{x}[/math] all of whose components are ones; then this SLAE is solved. In view of the general (with certain provisos) nature of the matrix, the solution of the system is divided into several stages. First, the [math]PLU[/math]-decomposition of [math]A[/math] is calculated by using Gaussian elimination with column pivoting ([math]P[/math] is a permutation matrix, [math]L[/math] is a lower triangular matrix with unit diagonal, and [math]U[/math] is an upper triangular matrix). At the second stage, the calculated decomposition is used to solve the SLAE [math]P L U \vec{x} = \vec{b}[/math]. At the third stage, the residual [math]\vec{r} = A \vec{x} - \vec{b}[/math] and the related characteristics are calculated. Then the data concerning the resulting accuracy and performance are produced. The performance is only estimated for the main part of the algorithm; that is, the calculations spent on the residual and norms are ignored.

The Linpack benchmark, which was developed for testing computers of traditional architecture, should not be confused with High Performance Linpack benchmark. The latter was devised for testing parallel computer systems and is based on a different (albeit similar) algorithm.

1.2 Mathematical description of the algorithm

Input data: square matrix [math]A[/math] (with entries [math]a_{ij}[/math]) (in 1992 version, the random number generator is chosen so that to ensure the non-singularity of [math]A[/math]), the right-hand side vector [math]\vec{b}[/math] (with components [math]b_i[/math]). Actually, they are not exactly input data (because the entries of the matrix are pseudo-random numbers generated by the program itself; these numbers are contained in the interval (-1/2, +1/2), and the right-hand sides are calculated as [math]b_i = \sum\limits_{j = 1}^i a_{ij}[/math]). However, in order to analyze the general scheme of the algorithm, it is convenient to regard [math]A[/math] and [math]\vec{b}[/math] as the input data.

Data to be calculated: lower triangular matrix [math]L[/math] (with entries [math]l_{ij}[/math]), upper triangular matrix [math]U[/math] (with entries [math]u_{ij}[/math]), permutation matrix [math]P[/math] (which is calculated as a collection of pivot indices rather than a full square array), and the solution vector [math]x[/math] (with components [math]x_i[/math]).

Here, we do not give detailed formulas of individual parts of the algorithm. These parts are independent algorithms, and their detailed descriptions should be presented separately. Consequently, we only discuss the general scheme of the algorithm.

The first part — decomposing the matrix into the product of two triangular matrices and a permutation matrix: [math]A = P L U[/math] — is performed by using Gaussian elimination with column pivoting. The permutation matrix is stored as a vector containing information on the performed permutations. The matrix [math]L[/math] is lower triangular with unit diagonal, while [math]U[/math] is an upper triangular matrix.

The second part is solving two triangular SLAEs. First, the system [math]L \vec{y} = P T \vec{b}[/math] is solved by forward substitution, then back substitution is used to solve the system [math]U \vec{x} = \vec{y}[/math]. This test stores both triangular matrices in the corresponding portions of the matrix [math]A[/math]. For this reason, [math]A[/math] is generated again (with the same starting parameters for the random number generator) before passing to the third part.

The third part (for which the performance is not measured) calculates first the residual of the solution [math]\vec{r} = A \vec{x} - \vec{b}[/math] and then its uniform norm, as well as the norm of the solution. On the completion of this stage, the algorithm outputs data related to the accuracy of the calculated solution and evaluates the performance of the system.

1.3 Computational kernel of the algorithm

The main part of the algorithm, which is the decomposition of a matrix via Gaussian elimination with column pivoting, is responsible for the bulk of operations. The second most important part is the calculation of the residual and its norm. It is implemented by the subroutine that calculates the sum of a vector and a matrix-vector product and by the function calculating the uniform norm. Less important parts are forward substitution in the Gauss method and back substitution in this method.

1.4 Macro structure of the algorithm

As already said in the description of the computational kernel of the algorithm, the decomposition of a matrix via Gaussian elimination with column pivoting is responsible for the bulk of operations. This decomposition is the first macro operation. It is followed by the forward substitution and, right away, by the back substitution. On the regeneration of the matrix, the next macro operation is calculating the residual for the solution of the SLAE, followed by the calculation the uniform norms of the residual and the solution.

1.5 Implementation scheme of the serial algorithm

The mathematical description of the algorithm was given above. The standard test calculates the triangular decomposition in-place, and the solutions to the triangular systems in the forward and back substitution are written to the locations of the right-hand sides. For this reason, on the calculation of the solution, the coefficient matrix and the right-hand side are regenerated for the subsequent calculation of the residual. As soon as the residual is calculated, its norm and the norm of the solution can be calculated in any order, which does not affect the algorithm and its results.

1.6 Serial complexity of the algorithm

Since the main part of the algorithm is the decomposition of a matrix via Gaussian elimination with column pivoting, the algorithm, as well as this decomposition, has cubical complexity. Other parts have quadratic complexity; consequently, they do not give a comparable contribution to the overall complexity of the algorithm.

1.7 Information graph

The graph consists of the information graphs of its basic parts. The first part (Gaussian elimination with column pivoting) has the most complex graph. Figure 1 shows the simplified graph of this part, post factum the choice of pivots. The operations of finding the maxima in the columns and saving the corresponding positions are given in green. The inversion of the current pivot, multiplications in the current column, and operations of the type [math]a b + c[/math] in the other columns are depicted in yellow, brown, and purple, respectively. Figure 2 shows the graph of a layer (one step of the method). The notation is the same as in the full graph. Dotted lines denote data communication from the pivot row.

Figure 1. Simplified graph of Gaussian elimination with column pivoting
Figure 2. Graph of one step of Gaussian elimination

The second part consists of almost identical fragments. The first fragment describes the forward substitution, namely, the solution of the SLAE with a lower triangular matrix [math]L \vec{y} = \vec{b}[/math] (Fig.3). The diagonal entries of this matrix are ones; consequently, this fragment contains no divisions. All the green nodes correspond to the operations of the type [math]a b + c[/math]. The second fragment describes the back substitution, namely, the solution of the SLAE with an upper triangular matrix [math]U \vec{x} = \vec{y}[/math] (Fig.4). This matrix has no specific features; hence, the associated graph contains division operations, which are shown in yellow.

Figure 3. Graph of solving the SLAE with a lower triangular matrix
Figure 4. Graph of solving the SLAE with an upper triangular matrix

1.8 Parallelization resource of the algorithm

The analysis of parallelization resource of the algorithm shows that Gaussian elimination with column pivoting has the highest, namely, quadratic complexity (exactly because of the pivoting). The other parts of the algorithm have at most linear complexity. Therefore, to optimize calculations, the virtual researcher must optimize Gaussian elimination with column pivoting. Besides, the first fragment of the second part (that is, the solution of the first triangular SLAE) can be attached to the first part and executed almost entirely in the background of the latter. Indeed, both have the same direction of recalculating the data. However, the rest of the second part, despite its similarity to the first fragment, can only be executed on completion the latter. The reason is that these two fragments have opposite directions of recalculating the data.

1.9 Input and output data of the algorithm

Actually, this test has always the same data, which are calculated in the test itself. Consequently, in the strict sense, it has no input data.

The output data, produced by the test, are the norms of the residual, solution, and matrix, the ratios of these norms, and the computational power of the target computer system.

1.10 Properties of the algorithm

2 Software implementation of the algorithm

2.1 Implementation peculiarities of the serial algorithm

Strictly speaking, the investigator implementing the test for a computer system tests at the same time LINPAC or any other package from which procedures for triangular decomposition and solving a SLAE with the help of this decomposition are called. The web versions of the test can also be used for testing other implementations of the Gauss method. It was already noted above that the triangular decomposition of the original matrix is calculated in-place, and, in the course of forward and back substitution, the solutions are written to the locations of right-hand sides. This causes the necessity to regenerate both the coefficient matrix and the right-hand side for the subsequent calculation of the residual. This organization of the algorithm may well be changed in order to optimize the test.

2.2 Possible methods and considerations for parallel implementation of the algorithm

2.3 Run results

2.4 Conclusions for different classes of computer architecture

3 References