High Performance Conjugate Gradient (HPCG) benchmark

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]\vec{p}_1 = \vec{z}_1[/math]

[math]\alpha_1 = (\vec{r}_0, \vec{z}_1)[/math]

иначе

[math]\alpha_i = (\vec{r}_{i-1}, \vec{z}_i)[/math]

[math]\beta_i = \frac{\alpha_i}{\alpha_{i - 1}}[/math]

[math]\vec{p}_i = \beta_i \vec{p}_{i - 1} + \vec{z}_i[/math]

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

[math]\gamma_i = \frac{\alpha_i}{(\vec{p}_i,A \vec{p}_i)}[/math]

[math]\vec{x}_i = \vec{x}_{i - 1} + \gamma_i \vec{p}_i[/math]

[math]\vec{r}_i = \vec{r}_{i - 1} - \gamma_i A \vec{p}_i[/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 Описание входных и выходных данных

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

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

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

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

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