High Performance Conjugate Gradient (HPCG) benchmark

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

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

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

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

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

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

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

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

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

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

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

• После этого вычисляются вектор $p_i$ и вспомогательные скаляры по формулам:

если $i = 1$, то

$\vec{p}_1 = \vec{z}_1$

$\alpha_1 = (\vec{r}_0, \vec{z}_1)$

иначе

$\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$

• Затем вычисляются новые приближения к вектору решения и невязки:

$\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$

• В заключение цикла вычисляется $|| \vec{r}_i ||_2$ и сравнивается с установленным пределом ошибки. Если норма ниже предела, процесс прерывается.

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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