Решатели СЛАУ с блочно-ленточными матрицами

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

СЛАУ с блочно-ленточной матрицей при численном решении задач математического моделирования возникает непосредственно при использовании сеточных методов или метода конечных элементов в случае прямоугольных областей и ортогональной сетки. Но во многих случаях сложных областей можно построить криволинейную, но ортогональную сетку [1] и прийти от СЛАУ с разреженной матрицей общего вида к СЛАУ с блочно-ленточной матрицей.

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

Для получения быстрого кода решателя предлагается разрабатывать и использовать специальный предкомпилятор. Предкомпилятор – это конвертор, преобразующий исходный высокоуровневый код в другой высокоуровневый код, который более эффективно отобразится на архитектуру вычислительной системы. Разработанный на основе ОРС (Оптимизирующей распараллеливающей системы [3]) предкомпиля-тор ускоряет код решателя более чем на 30% для тестового набора данныx.

1.    Итерационные алгоритмы для решения целевых СЛАУ

Для полноты изложения приведем необходимую информацию об итерационных алгоритмах [4]. Поскольку предполагается рассматривать решатели для СЛАУ большой размерности, чтобы обеспечить приемлемую точность и не выходить за пределы имеющейся памяти, необходимо использовать итерационные алгоритмы. Как обычно, при решении задач большой размерности необходимо сбалансировать соблюдение точности (не должна накапливаться погрешность), ограничения по памяти и время работы вычислительной системы. Время работы итерационного процесса зависит от количества итераций и времени выполнения одной итерации.

Систему уравнений Ax = b o будем решать итерационным алгоритмом

x ^(k+1) = Bx ( k + b,

B = I — tC ^-1 A, b = tC ^-1 b o .

Здесь 0,1, 2,..., k - номер итерации, A - квадратная матрица размера n х n, b - вектор размера n, x ⁽⁰⁾ - начальное приближение, которое задает пользователь, x ^(k) -результат выполнения k -ой итерации, C – произвольная невырожденная матрица, t – положительный числовой параметр. Для сходимости итерационного процесса спектральный радиус матрицы B должен быть меньше единицы. Основная отличительная особенность процессоров последних лет состоит в том, что время выполнения арифметических операций более чем на порядок дольше времени считывания аргументов таких операций из оперативной памяти [5].

Можно сформулировать условия для выбора параметров итерационного процесса C и t >  0:1) вычисление C ^-1 x должно требовать мало времени; 2) объем памяти для матриц B и C должен быть равен или ненамного больше объема памяти для хранения матрицы A ; 3) спектральная норма матрицы B должна быть меньше 1.

Матрица B обратима для всех малых чисел t таких, что число t — 1 больше наибольшего собственного значения матрицы C ^-1 A. В качестве t можно взять число ¹

t = ^-1||—|pjjj — ePS (здесь eps - некоторое заданное положительное число), поскольку спектральная норма матрицы C-1A не больше любой нормы этой матрицы, которая не больше ||C-1|| • ||A||. Используется тот факт, что спектр полинома от линейного оператора равен полиному от его спектра (см. [3]) применительно к оператору I — tC-1A, т.е. Sp{I — tC-1A} = (1 — t Sp{C-1A}). Выбор параметра t позволит соблюсти требование 3). От этого параметра зависит количество шагов итерационного процесса.

В качестве матрицы C будем рассматривать блочную ленточную матрицу с малым количеством ненулевых диагоналей 2 • к + 1 близких к главной. Эта матрица получается из матрицы A обнулением всех элементов A ij , для которых | i — j | > к. Матрицу A представим в виде A = C + O, где O - матрица, состоящая из оставшихся элементов матрицы A после обнуления (вычитания) элементов по главной диагонали и k соседних (с двух сторон) диагоналей. Сделав LU -разложение матрицы C (или C = LDL* для симметричной матрицы), можно быстро вычислять C ^{- 1} x на каждой итерации алгоритма. Во многих приложениях матрица A является симметрической или положительно определенной. В этих случаях матрица C тоже может быть симметрической или положительно определенной и, как следствие, матрица O – симметрическая или положительно определенная соответственно. Может быть, что A – не

Б.Я. Штейнберг, А.А. Василенко, В.В. Веселовский, Н.А. Живых симметрическая или не положительно определенная (например, с седловой особенностью), но C – симметрическая или положительно определенная. Если у A было преобладание по главной диагонали, то и C обладает этим же свойством.

Рассмотрим сложность вычисления одной итерации алгоритма:

Bx + b = ( I - tC - 1 A ) x + b = (1 - t)x - tC - 1 Ox + b.              (3)

Основное время вычисления одной итерации занимает вычисление вектора C ^-1 Ox .

2.    Хранение блочно-ленточных матриц

Разреженные матрицы хранятся в памяти [6] только ненулевыми элементами и номерами их строк и столбцов (итого, 3 числа). В некоторых случаях заранее известно, что у разреженной матрицы все ненулевые элементы сосредоточены на нескольких (м.б. блочных) диагоналях. Если хранить только эти диагонали, то расход памяти для хранения такой матрицы будет меньше, чем в предыдущем случае, поскольку для таких элементов не следует запоминать номера строк и столбцов.

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

Если построить сетку на (трехмерном) прямоугольном параллелепипеде, узлы которой занумеровать соответствующим образом, матрица конечно-разностной схемы также получится ленточной. Можно рассматривать отображения прямоугольного параллелепипеда с координатной сеткой на другие области, перенося на них сетку.

/□ □        □\

• □    □ □□

• □    □ □□

....

....

■•. ■•. ■•.

□

□

.

.

.

• □       □ □ □

• □       □ □)

3.    Компилятор для ускорения решателей и результатычисленных экспериментов

Пятидиагональная блочно-ленточная матрица

Если распределенную матрицу хранить тройками чисел ( i, j, x [ i,j ]), где целые (4 байта) i, j – номера строки и столбца, в которых хранится ненулевой элемент матрицы x [ i, j ] (тип double), то объем памяти для одного ненулевого элемента равен 16 байт (4 + 4 + 8). Ленточную матрицу можно хранить диагоналями. Каждая диагональ -это одномерный массив double. Объем памяти получается в 2 раза меньше.

Если рассматривается не одно дифференциальное уравнение, а система дифференциальных уравнений на регулярной области, то матрица СЛАУ получается блочно-ленточной. Элементы матрицы СЛАУ – блоки (матрицы), размерность которых равна количеству уравнений в исходной системе дифференциальных уравнений. Такие же матрицы возникают и при использовании метода конечных элементов.

Каждая блочная диагональ блочно-ленточной матрицы хранится в памяти одномерным массивом. Элементы каждого блока размещаются в оперативной памяти рядом построчно. Такое размещение приводит к минимизации кэш-промахов. Диагональ матрицы A состоит из всех элементов A ij , для которых (i — j ) = к. Это позволяет легко по номеру строки элемента находить номер столбца.

Программы решения СЛАУ с блочно-ленточной матрицей имеют особенности, учитывая которые можно повысить производительность алгоритма. После LU-разложения по главной диагонали у матриц L и U часть программных циклов имеет треугольный вид. Блоки по главной диагонали имеют треугольный вид. Разложение симметричной матрицы можно получить методом Холецкого [7]. Блоки по главной диагонали нижне-треугольной матрицы L обозначим D[i],i = 1, 2, ...,n. Эти блоки имеют нижне-треугольный вид. Блоки, расположенные на поддиагонали, обозначим B [ i ] ,i = 2, 3, ...,n. Обозначая x[i] и h [ i ] ,i = 1, 2, ...,n блочные векторы неизвестных и правой части соответственно, получим блочную систему уравнений:

D [1] • x[1] = h[1]; D[i] • x[i] + B[i] • x [ i — 1] = h[i], i = 2, 3,..., n. (4)

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

Блочные ленты (диагонали) хранятся в памяти в виде одномерных массивов. Это улучшает локализацию данных, но усложняет индексы элементов такого массива. Индексные выражения элементов массива аффинно зависят от счетчиков циклов, в которых находятся вхождения этих элементов. Оптимизация вычисления таких индексных выражений может достигаться преобразованиями ≪ вычисление общих под-выражений ≫ , ≪ вынос инвариантов из цикла ≫ [8–10].

4.    Численные эксперименты

Авторами статьи разработан экспериментальный решатель СЛАУ с блочноленточными матрицами. Решатель ориентирован на внедрение в конечно-элементный пакет ACELAN-COMPOS. Задача программы – нахождение решений СЛАУ для сетки 100 x 100 x 100. Для текущей задачи предполагаются блоки размера 4. Задача трехмерная, поэтому в простейшем случае количество блочных диагоналей равно 7. Вычисление нормы |x(k+1) — xk| занимает время, сопоставимое со временем одной итерации, поэтому проверка проводится через группу итераций (выбран размер группы 50). Итого, количество ненулевых элементов матрицы примерно 7 • 42 • 106 ^ 108. Матрица системы уравнений имеет седловую особенность, положительных собствен- ных чисел – 4 · 666666, отрицательных собственных чисел – 4 · 333333. Подматрица с положительными собственными числами и подматрица с отрицательными собственными числами имеют преобладание по главной диагонали. Учитывая симметрию, в памяти можно хранить половину элементов 5 · 107 восьмибайтных чисел (double). На каждой итерации итерационного алгоритма каждый элемент матрицы используется 1 раз. Итого, при выполнении алгоритма производится 50 · 108 = 5 · 109 операций над числами типа double.

Программа решения СЛАУ написана на языке С (С89). С помощью ОРС (оптимизирующей распараллеливающей системы) разработанная Программа (исходная) преобразована к Программе-2 (оптимизированная). Оптимизация ОРС выполняет последовательность преобразований из набора ≪ линеаризация ≫ , ≪ вынос инвариантных выражений из циклов ≫ , ≪ раскрутка циклов ≫ , ≪ канонизация ≫ .

Компьютер, на котором проводились численные эксперименты, имеет следующие характеристики: CPU - i7-9700; 3,00 GHz; L1 - 256 Kb; L2 - 2 Mb; L3 - 12 Mb; частота системной шины - 8 GT/s; максимальная пропускная способность - 41,6 GB/s; RAM - DDR4, 16Gb; мин. частота - 1600 МГц; макс. частота - 2666 МГц. Из таблицы численных экспериментов видно, что оптимизирующий предкомпилятор ускоряет решатель на 30% (при использовании компилятора ICC).

Таблица

Результаты численных экспериментов

Компилятор, ключ оптимизации	GCC, O2	GCC, O3	ICC, O2	ICC, O3	MSVS, O2
Программа , c	5,0357	5,3133	7,0142	6,7631	10,1573
Программа-2 , c	4,9197	4,6092	4,2593	4,2085	9,2703

Заключение

Распараллеливание пока не проводилось. При решении СЛАУ с треугольными матрицами возникают циклы с линейной рекуррентной зависимостью [11], распараллеливание которых предполагается реализовать в следующей работе. Высокое быстродействие представленного решателя достигнуто несколькими действиями. Некоторые оптимизации выполнены предкомпилятором. Существенное ускорение достигнуто за счет представления матрицы системы блочными диагоналями и за счет хранения блоков в одномерном массиве.

За многие годы накоплены сотни пакетов прикладных программ, как с открытым кодом, так и коммерческих. Большинство этих пакетов нуждаются в быстроте их решателей. Разумеется, решатели, написанные более 10 лет назад, не ориентированы на современные процессоры. Более того, решатель, оптимизированный под один процессор, может неэффективно отображаться на другой процессор, даже выпущенный тем же производителем. Такая ситуация приводит к проблеме эффективной переносимости кода. Создание новых процессоров с более сложными архитектурами обостряет проблему переносимости кода с сохранением высокой эффективности. Эту проблему отчасти могли бы решать оптимизирующие компиляторы нового поколения, а в первое время могут использоваться специальные предкомпиляторы для решателей пакетов.

Работа проводилась при финансовой поддержке гранта Правительства РФ № 075-15-2019-1928.