Уровень алгоритма

High Performance Conjugate Gradient (HPCG) benchmark: различия между версиями

Материал из Алговики
Перейти к навигации Перейти к поиску
[непроверенная версия][досмотренная версия]
(Перенесено из HPCG-0.1.9.docx.)
 
 
(не показано 19 промежуточных версий 5 участников)
Строка 1: Строка 1:
== Описание свойств и структуры алгоритма ==
+
{{level-a}}
  
=== Словесное описание алгоритма ===
+
Основные авторы описания: [[Участник:Frolov|А.В.Фролов]].
 +
 
 +
== Свойства и структура алгоритма ==
 +
 
 +
=== Общее описание алгоритма ===
  
 
HPCG является тестом, разработанным для тестирования производительности компьютерных систем. По сути это генерирование некоторой системы линейных алгебраических уравнений (СЛАУ) <math>A \vec{x} = \vec{b}</math> с разрежённой квадратной положительно определённой симметричной матрицей <math>A</math> и вектором <math>\vec{b}</math>, с последующим решением этой СЛАУ. В качестве метода решения выбран один из методов сопряжённых направлений, основанных на ортогонализации последовательностей Крылова: метод сопряжённых градиентов с предобуславливателем Гаусса-Зейделя. В качестве тестируемой СЛАУ используется дискретизация уравнения диффузии в 3-мерной области с такими условиями на границе и в правой части, чтобы правильное решение СЛАУ состояло из одних единиц. При этом на одну строку (столбец) приходится по 27 ненулевых элементов («внутри» матрицы) и от 7 до 18 на «приграничных» строках и столбцах. Всего в тесте 48 384 000 линейных уравнений, причём матрица содержит 1 298 936 872 ненулевых элемента.
 
HPCG является тестом, разработанным для тестирования производительности компьютерных систем. По сути это генерирование некоторой системы линейных алгебраических уравнений (СЛАУ) <math>A \vec{x} = \vec{b}</math> с разрежённой квадратной положительно определённой симметричной матрицей <math>A</math> и вектором <math>\vec{b}</math>, с последующим решением этой СЛАУ. В качестве метода решения выбран один из методов сопряжённых направлений, основанных на ортогонализации последовательностей Крылова: метод сопряжённых градиентов с предобуславливателем Гаусса-Зейделя. В качестве тестируемой СЛАУ используется дискретизация уравнения диффузии в 3-мерной области с такими условиями на границе и в правой части, чтобы правильное решение СЛАУ состояло из одних единиц. При этом на одну строку (столбец) приходится по 27 ненулевых элементов («внутри» матрицы) и от 7 до 18 на «приграничных» строках и столбцах. Всего в тесте 48 384 000 линейных уравнений, причём матрица содержит 1 298 936 872 ненулевых элемента.
  
=== Математическое описание ===
+
=== Математическое описание алгоритма ===
  
 
Исходные данные: разрежённая квадратная положительно определённая симметричная матрица <math>A</math> (элементы <math>a_{ij}</math>), вектор правой части <math>\vec{b}</math> (элементы <math>b_i</math>). В реальности они не вполне входные (элементы матрицы генерирует сама программа), но для разбора общей схемы алгоритма удобно считать их входными данными.
 
Исходные данные: разрежённая квадратная положительно определённая симметричная матрица <math>A</math> (элементы <math>a_{ij}</math>), вектор правой части <math>\vec{b}</math> (элементы <math>b_i</math>). В реальности они не вполне входные (элементы матрицы генерирует сама программа), но для разбора общей схемы алгоритма удобно считать их входными данными.
Строка 25: Строка 29:
 
* После этого вычисляются вектор <math>p_i</math> и вспомогательные скаляры по формулам:
 
* После этого вычисляются вектор <math>p_i</math> и вспомогательные скаляры по формулам:
 
:если <math>i = 1</math>, то
 
:если <math>i = 1</math>, то
:<math>\vec{p}_1 = \vec{z}_1</math>
+
:<math>
:<math>\alpha_1 = (\vec{r}_0, \vec{z}_1)</math>
+
\begin{align}
 +
\vec{p}_1 & = \vec{z}_1 \\
 +
\alpha_1 & = (\vec{r}_0, \vec{z}_1)
 +
\end{align}
 +
</math>
 
:иначе
 
:иначе
:<math>\alpha_i = (\vec{r}_{i-1}, \vec{z}_i)</math>
+
:<math>
:<math>\beta_i = \frac{\alpha_i}{\alpha_{i - 1}}</math>
+
\begin{align}
:<math>\vec{p}_i = \beta_i \vec{p}_{i - 1} + \vec{z}_i</math>
+
\alpha_i & = (\vec{r}_{i-1}, \vec{z}_i) \\
 +
\beta_i   & = \frac{\alpha_i}{\alpha_{i - 1}} \\
 +
\vec{p}_i & = \beta_i \vec{p}_{i - 1} + \vec{z}_i
 +
\end{align}
 +
</math>
  
 
* Затем вычисляются новые приближения к вектору решения и невязки:
 
* Затем вычисляются новые приближения к вектору решения и невязки:
:<math>\gamma_i = \frac{\alpha_i}{(\vec{p}_i,A \vec{p}_i)}</math>
+
:<math>
:<math>\vec{x}_i = \vec{x}_{i - 1} + \gamma_i \vec{p}_i</math>
+
\begin{align}
:<math>\vec{r}_i = \vec{r}_{i - 1} - \gamma_i A \vec{p}_i</math>
+
\gamma_i & = \frac{\alpha_i}{(\vec{p}_i,A \vec{p}_i)} \\
 +
\vec{x}_i & = \vec{x}_{i - 1} + \gamma_i \vec{p}_i \\
 +
\vec{r}_i & = \vec{r}_{i - 1} - \gamma_i A \vec{p}_i
 +
\end{align}
 +
</math>
  
 
* В заключение цикла вычисляется <math>|| \vec{r}_i ||_2</math> и сравнивается с установленным пределом ошибки. Если норма ниже предела, процесс прерывается.
 
* В заключение цикла вычисляется <math>|| \vec{r}_i ||_2</math> и сравнивается с установленным пределом ошибки. Если норма ниже предела, процесс прерывается.
Строка 69: Строка 85:
 
** вычисление квадратичной нормы вектора.
 
** вычисление квадратичной нормы вектора.
  
=== Описание схемы реализации последовательного алгоритма ===
+
=== Схема реализации последовательного алгоритма ===
  
 
Математическое описание алгоритма описано выше. При этом в ходе вычислений, реализующих какое-либо вычисление произведений матрицы на вектор (включая вычисления невязок или решения треугольных СЛАУ) следует учитывать, что для вектора размерности <math>n</math> каждая его компонента умножается не более чем на 27 элементов матрицы.
 
Математическое описание алгоритма описано выше. При этом в ходе вычислений, реализующих какое-либо вычисление произведений матрицы на вектор (включая вычисления невязок или решения треугольных СЛАУ) следует учитывать, что для вектора размерности <math>n</math> каждая его компонента умножается не более чем на 27 элементов матрицы.
Строка 83: Строка 99:
 
Граф состоит из информационных графов своих основных частей, поэтому здесь нет смысла приводить его целиком: совместное представление только усложнило бы общую картину, которую, вне сомнения, следует изучать по графам отдельных частей - методов, составляющих этот алгоритм. При этом следует отметить, что одной из таких отдельных частей будет вычисление предобусловленной невязки, а не её внутренние части типа решения треугольных СЛАУ и невязки треугольной матрицы. Это связано с тем, что эти внутренние части вычисления предобусловленной невязки хорошо стыкуются друг с другом.
 
Граф состоит из информационных графов своих основных частей, поэтому здесь нет смысла приводить его целиком: совместное представление только усложнило бы общую картину, которую, вне сомнения, следует изучать по графам отдельных частей - методов, составляющих этот алгоритм. При этом следует отметить, что одной из таких отдельных частей будет вычисление предобусловленной невязки, а не её внутренние части типа решения треугольных СЛАУ и невязки треугольной матрицы. Это связано с тем, что эти внутренние части вычисления предобусловленной невязки хорошо стыкуются друг с другом.
  
=== Описание ресурса параллелизма алгоритма ===
+
=== Ресурс параллелизма алгоритма ===
  
 
При изучении ресурса параллелизма видно, что основная ''сложность (линейная)'' будет у вычислений скалярных произведений. Остальные части алгоритма в силу разреженности матриц имеют не более постоянной сложности с небольшим коэффициентом. Поэтому для оптимизации вычислений нужно заняться именно оптимизацией вычисления скалярных произведений.
 
При изучении ресурса параллелизма видно, что основная ''сложность (линейная)'' будет у вычислений скалярных произведений. Остальные части алгоритма в силу разреженности матриц имеют не более постоянной сложности с небольшим коэффициентом. Поэтому для оптимизации вычислений нужно заняться именно оптимизацией вычисления скалярных произведений.
  
=== Описание входных и выходных данных ===
+
=== Входные и выходные данные алгоритма ===
  
 
По сути теста у него всегда одни и те же данные, причём вычисляемые в самом же тесте. Поэтому в строгом смысле у алгоритма нет входных данных.  
 
По сути теста у него всегда одни и те же данные, причём вычисляемые в самом же тесте. Поэтому в строгом смысле у алгоритма нет входных данных.  
Строка 98: Строка 114:
  
 
=== Особенности реализации последовательного алгоритма ===
 
=== Особенности реализации последовательного алгоритма ===
=== Описание локальности данных и вычислений ===
+
 
==== Описание локальности алгоритма ====
+
=== Возможные способы и особенности параллельной реализации алгоритма ===
==== Описание локальности реализации алгоритма ====
+
 
===== Описание структуры обращений в память и качественная оценка локальности =====
+
Существующая реализация теста одна. При этом следует отметить важный момент: для параллельных систем версия программы представляет собой алгоритм, отличающийся от последовательного. Разработчиками заявлено, что это сделано для оптимизации распараллеливания программы, и при этом количество требующихся для сходимости к нужной точности итераций увеличивается в сравнении с последовательной версией. Таким образом, под одним названием существуют как минимум 2 алгоритма.
===== Количественная оценка локальности =====
+
 
===== Анализ на основе теста Apex-Map =====
+
=== Результаты прогонов ===
=== Возможные способы и особенности реализации параллельного алгоритма ===
 
=== Масштабируемость алгоритма и его реализации ===
 
==== Описание масштабируемости алгоритма ====
 
==== Описание масштабируемости реализации алгоритма ====
 
=== Динамические характеристики и эффективность реализации алгоритма ===
 
 
=== Выводы для классов архитектур ===
 
=== Выводы для классов архитектур ===
=== Существующие реализации алгоритма ===
 
  
Существующая реализация теста одна. При этом следует отметить важный момент: для параллельных систем версия программы представляет собой алгоритм, отличающийся от последовательного. Разработчиками заявлено, что это сделано для оптимизации распараллеливания программы, и при этом количество требующихся для сходимости к нужной точности итераций увеличивается в сравнении с последовательной версией. Таким образом, под одним названием существуют как минимум 2 алгоритма.
+
== Литература ==
 +
 
 +
<references />
 +
 
 +
[[Категория:Статьи в работе]]
 +
 
 +
[[en:High Performance Conjugate Gradient (HPCG) benchmark]]

Текущая версия на 11:33, 12 июля 2022


Основные авторы описания: А.В.Фролов.

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

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

HPCG является тестом, разработанным для тестирования производительности компьютерных систем. По сути это генерирование некоторой системы линейных алгебраических уравнений (СЛАУ) [math]A \vec{x} = \vec{b}[/math] с разрежённой квадратной положительно определённой симметричной матрицей [math]A[/math] и вектором [math]\vec{b}[/math], с последующим решением этой СЛАУ. В качестве метода решения выбран один из методов сопряжённых направлений, основанных на ортогонализации последовательностей Крылова: метод сопряжённых градиентов с предобуславливателем Гаусса-Зейделя. В качестве тестируемой СЛАУ используется дискретизация уравнения диффузии в 3-мерной области с такими условиями на границе и в правой части, чтобы правильное решение СЛАУ состояло из одних единиц. При этом на одну строку (столбец) приходится по 27 ненулевых элементов («внутри» матрицы) и от 7 до 18 на «приграничных» строках и столбцах. Всего в тесте 48 384 000 линейных уравнений, причём матрица содержит 1 298 936 872 ненулевых элемента.

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

Исходные данные: разрежённая квадратная положительно определённая симметричная матрица [math]A[/math] (элементы [math]a_{ij}[/math]), вектор правой части [math]\vec{b}[/math] (элементы [math]b_i[/math]). В реальности они не вполне входные (элементы матрицы генерирует сама программа), но для разбора общей схемы алгоритма удобно считать их входными данными.

Вычисляемые данные: вектор решения [math]\vec{x}[/math] (элементы [math]x_i[/math]), вектор [math]\vec{r}[/math] (элементы [math]r_i[/math]) невязки.

Детальные формулы отдельных частей алгоритма здесь описывать не будем, поскольку они являются самостоятельными алгоритмами и их детальные описания должны быть выполнены раздельно. Поэтому здесь будет приведена только общая схема алгоритма. Обращаем внимание на то, что в учебной литературе обычно приводится схема метода сопряжённых градиентов без предобуславливателя. Введение предобуславливателя объясняется тем, что без него матрица системы [math]A \vec{x} = \vec{b}[/math] может быть плохо обусловленной, и это увеличивает необходимое для окончания работы метода сопряжённых градиентов количество итераций сверх меры. Дело в том, что, хотя метод сопряжённых градиентов по своему построению и математическому обоснованию является прямым (то есть у него известно количество операций, после которого он наверняка сойдётся к решению), но из-за того, что он использует точные формулы, выполняющиеся в процессах ортогонализации только в точной арифметике, появляющаяся неустойчивость настолько критична, что становится необходимым уменьшить количество итераций по сравнению с верхней теоретической оценкой. Это возможно только существенным уменьшением обусловленности матрицы системы с помощью матрицы [math]M[/math], называемой пере(пред-)обуславливателем. Эта матрица должна быть «приближением» матрицы [math]A[/math] и в то же время легко обращаться. Тогда вместо плохо обусловленной СЛАУ [math]A \vec{x} = \vec{b}[/math] будет решаться хорошо обусловленная СЛАУ [math]M^{-1} A \vec{x} = M^{-1} \vec{b}[/math]. Само собой, фактической замены СЛАУ никто не производит, поскольку матрица [math]M^{-1} A[/math] не является разрежённой. Вместо этого «умножение на [math]M^{-1}[/math]» производят, когда в методе сопряжённых уравнений переходят от невязок [math]\vec{r}_i[/math] к векторам [math]M^{-1} \vec{r}_i[/math].

Стартовая часть: программа генерирует матрицу [math]A[/math] в некотором структурном виде и вектор правой части [math]b[/math].

Начало процесса: вычисляется стартовый вектор невязки [math]\vec{r}_0 = \vec{b} - A \vec{x}_0[/math], в качестве другого стартового вектора [math]p_0[/math] используется вектор [math]x_0[/math].

После этого стартует циклический процесс. Номер исполнения тела цикла обозначим [math]i[/math]. В цикле выполняются следующие операции:

  • Выполняется вычисление предобусловленной невязки [math]z_i := M^{-1} \vec{r}_{i - 1}[/math]. Это происходит в три этапа, в которых используется разбиение матрицы A на две треугольные и диагональную ([math]L[/math] - нижний треугольник [math]A[/math], [math]U[/math] - верхний треугольник [math]A[/math], [math]D[/math] — диагональ [math]A[/math] (а также диагональ матриц [math]L[/math] и [math]U[/math]), так что [math]A = L + U - D[/math], [math]U = L^T[/math]):
    1. Вычисляется [math]\vec{t} = L^{-1} \vec{r}_{i - 1}[/math]
    2. Вычисляется [math]\vec{s} = \vec{r}_{i-1} - (L - D) \vec{t}[/math]
    3. Вычисляется [math]\vec{z}_i = U^{-1} \vec{s}[/math].
  • После этого вычисляются вектор [math]p_i[/math] и вспомогательные скаляры по формулам:
если [math]i = 1[/math], то
[math] \begin{align} \vec{p}_1 & = \vec{z}_1 \\ \alpha_1 & = (\vec{r}_0, \vec{z}_1) \end{align} [/math]
иначе
[math] \begin{align} \alpha_i & = (\vec{r}_{i-1}, \vec{z}_i) \\ \beta_i & = \frac{\alpha_i}{\alpha_{i - 1}} \\ \vec{p}_i & = \beta_i \vec{p}_{i - 1} + \vec{z}_i \end{align} [/math]
  • Затем вычисляются новые приближения к вектору решения и невязки:
[math] \begin{align} \gamma_i & = \frac{\alpha_i}{(\vec{p}_i,A \vec{p}_i)} \\ \vec{x}_i & = \vec{x}_{i - 1} + \gamma_i \vec{p}_i \\ \vec{r}_i & = \vec{r}_{i - 1} - \gamma_i A \vec{p}_i \end{align} [/math]
  • В заключение цикла вычисляется [math]|| \vec{r}_i ||_2[/math] и сравнивается с установленным пределом ошибки. Если норма ниже предела, процесс прерывается.

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

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

  • Вычисление стартовой невязки
  • и затем циклически:
    • «прямая» обратная подстановка для решения СЛАУ с нижней треугольной матрицей,
    • новое вычисление «невязки», но не с полной, а с нижней треугольной матрицей,
    • обратная подстановка для решения СЛАУ с верхней треугольной матрицей,
    • вычисление скалярного произведения 2 векторов,
    • вычисление взвешенной суммы двух векторов,
    • вычисление произведения матрицы на вектор,
    • вычисление скалярного произведения двух векторов,
    • вычисление взвешенной суммы двух векторов,
    • вычисление квадратичной нормы вектора.

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

Как уже записано в описании ядра алгоритма, макроструктура алгоритма состоит в последовательном выполнении следующих макроопераций:

  • Вычисление стартовой невязки,
  • и затем циклически:
    • «прямая» обратная подстановка для решения СЛАУ с нижней треугольной матрицей,
    • новое вычисление «невязки», но не с полной, а с нижней треугольной матрицей,
    • обратная подстановка для решения СЛАУ с верхней треугольной матрицей,
    • вычисление скалярного произведения 2 векторов,
    • вычисление взвешенной суммы двух векторов,
    • вычисление произведения матрицы на вектор,
    • вычисление скалярного произведения двух векторов,
    • вычисление взвешенной суммы двух векторов,
    • вычисление квадратичной нормы вектора.

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

Математическое описание алгоритма описано выше. При этом в ходе вычислений, реализующих какое-либо вычисление произведений матрицы на вектор (включая вычисления невязок или решения треугольных СЛАУ) следует учитывать, что для вектора размерности [math]n[/math] каждая его компонента умножается не более чем на 27 элементов матрицы.

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

Поскольку для вектора размерности [math]n[/math] каждая его компонента умножается не более чем на 27 элементов матрицы, то все макрооперации, связанные с умножением матрицы на вектор, имеют линейную последовательную сложность с коэффициентом 27 (или менее). Вычисления скалярных произведений (в т. ч. квадратичной нормы невязки) также линейны по последовательной сложности. Таким образом, общая последовательная сложность алгоритма линейна. Коэффициент пропорциональности между [math]n[/math] и сложностью в общем случае зависит не только от разрежённости матриц, но и от количества итераций, требуемых для достижения требуемой точности решения.

В одном вычислении невязки для конкретных матриц теста сложность составит около [math]27 n[/math] операций типа «умножить и сложить». Вычисление скалярного произведения даёт порядка [math]n[/math] таких операций. Таким образом, для тела цикла последовательная сложность составит порядка [math]112 n[/math] операций умножения и сложения.

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

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

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

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

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

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

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

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

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

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

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

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

2.3 Результаты прогонов

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

3 Литература