High Performance Conjugate Gradient (HPCG) benchmark

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

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

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

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

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

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

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

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

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

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

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

• Выполняется вычисление предобусловленной невязки [math]\displaystyle{ z_i := M^{-1} \vec{r}_{i - 1} }[/math]. Это происходит в три этапа, в которых используется разбиение матрицы A на две треугольные и диагональную ([math]\displaystyle{ L }[/math] - нижний треугольник [math]\displaystyle{ A }[/math], [math]\displaystyle{ U }[/math] - верхний треугольник [math]\displaystyle{ A }[/math], [math]\displaystyle{ D }[/math] — диагональ [math]\displaystyle{ A }[/math] (а также диагональ матриц [math]\displaystyle{ L }[/math] и [math]\displaystyle{ U }[/math]), так что [math]\displaystyle{ A = L + U - D }[/math], [math]\displaystyle{ U = L^T }[/math]):
  1. Вычисляется [math]\displaystyle{ \vec{t} = L^{-1} \vec{r}_{i - 1} }[/math]
  2. Вычисляется [math]\displaystyle{ \vec{s} = \vec{r}_{i-1} - (L - D) \vec{t} }[/math]
  3. Вычисляется [math]\displaystyle{ \vec{z}_i = U^{-1} \vec{s} }[/math].

• После этого вычисляются вектор [math]\displaystyle{ p_i }[/math] и вспомогательные скаляры по формулам:

если [math]\displaystyle{ i = 1 }[/math], то

[math]\displaystyle{ \begin{align} \vec{p}_1 & = \vec{z}_1 \\ \alpha_1 & = (\vec{r}_0, \vec{z}_1) \end{align} }[/math]

иначе

[math]\displaystyle{ \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]\displaystyle{ \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]\displaystyle{ || \vec{r}_i ||_2 }[/math] и сравнивается с установленным пределом ошибки. Если норма ниже предела, процесс прерывается.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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