Алгоритмы решения систем линейных уравнений с блочно-ленточными матрицами

Ростов-на-Дону — Таганрог
Издательство Южного федерального университета
2022

МИНИСТЕРСТВО НАУКИ И ВЫСШЕГО ОБРАЗОВАНИЯ РОССИЙСКОЙ ФЕДЕРАЦИИ
ФЕДЕРАЛЬНОЕ ГОСУДАРСТВЕННОЕ АВТОНОМНОЕ  
ОБРАЗОВАТЕЛЬНОЕ УЧРЕЖДЕНИЕ ВЫСШЕГО ОБРАЗОВАНИЯ 
«ЮЖНЫЙ ФЕДЕРАЛЬНЫЙ УНИВЕРСИТЕТ»

Б. Я. Штейнберг, 
О. Б. Штейнберг

АЛГОРИТМЫ РЕШЕНИЯ  
СИСТЕМ ЛИНЕЙНЫХ УРАВНЕНИЙ  
С БЛОЧНО-ЛЕНТОЧНЫМИ 
МАТРИЦАМИ

Монография

УДК 517.983:[512.644:519.612] (035.3)
ББК 22.143+22.162.3 я44
Ш 88
Печатается по решению Комитета при Ученом совете  
Южного федерального университета по естественнонаучному 
и математическому направлению науки  и образования  
(протокол № 10 от 9 июня 2021 г.)
Рецензенты:
профессор кафедры информатики и вычислительного эксперимента  
ИММ и КН им. И. И. Воровича Южного федерального университета,  
доктор физико-математических наук Г. В. Муратова;
заведующий лабораторией функционально-градиентных и композитных 
материалов научно-образовательного центра «Материалы» ДГТУ  
доктор физико-математических наук С. М. Айзикович

 
Штейнберг, Б. Я.
Ш 88  
Алгоритмы решения систем линейных уравнений с блочно-ленточными 
матрицами : монография / Б. Я. Штейнберг, О. Б. Штейнберг  ; 
Южный федеральный университет. — Ростов-на-Дону ; Таганрог : 
 Издательство Южного федерального университета, 2022. — 138 с.

ISBN 978-5-9275-4061-7 
DOI 10.18522/801287963

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

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

ISBN 978-5-9275-4061-7 
УДК 517.983:[512.644:519.612] (035.3) 
 
ББК 22.143+22.162.3 я44

 
© Южный федеральный университет, 2022
 
© Штейнберг Б. Я., Штейнберг О. Б., 2022
 
© Оформление. Макет. Издательство 
Южного федерального университета, 2022

СОДЕРЖАНИЕ

Список сокращений .........................................................................................5

Введение .............................................................................................................6

1. Влияние вычислительных архитектур и компиляторов 
на быстродействие алгоритмов и их программных реализаций .......9
1.1. Разнообразие и особенности вычислительных систем .......................9
1.2. Модели параллельных вычислений ....................................................16
1.3. Анализ узких мест производительности программ ..........................19
1.4. Отставание развития оптимизирующих компиляторов  
от развития процессоров ......................................................................20
1.5. Информационные зависимости в программе ....................................22
1.6. Гнезда циклов и пространство итераций ............................................24
1.7. Блочно-ленточные матрицы .................................................................25

2. Преобразования программ и предкомпилятор для ускорения 
решения СЛАУ с блочно-ленточными матрицами .............................29
2.1. Некоторые эквивалентные преобразования программ .....................29
2.2. Предкомпилятор для ускорения алгоритмов ......................................36

3. О параллельных прямых алгоритмах решения СЛАУ  
с ленточными и блочно-ленточными матрицами ...............................37
3.1. Параллельное решение СЛАУ с двухдиагональными и блочными 
двухдиагональными матрицами ..........................................................38
3.2. О последовательном алгоритме, соответствующем параллельному 
алгоритму решения СЛАУ с двухдиагональной матрицей ...............40
3.3. Применения алгоритма параллельного решения системы  
уравнения с ленточной двухдиагональной матрицей .......................43
3.4. СЛАУ с ленточными треугольными матрицами ................................45
3.5. СЛАУ с блочной трехдиагональной не треугольной матрицей .......47
3.6. СЛАУ с ленточной матрицей с преобладанием по одной 
из диагоналей ........................................................................................51
3.7. Параллельное вычисление дробнолинейных рекуррентных 
последовательностей ............................................................................55
3.8. Разложения трехдиагональных матриц в произведение 
двухдиагональных ................................................................................57
3.9. Оценки погрешностей решения СЛАУ с двухдиагональной 
матрицей ................................................................................................59

4. Распараллеливание рекуррентных циклов ..........................................64
4.1. Рекуррентно вычисляемые переменные и рекуррентные циклы .....64
4.2. Рекуррентно замкнутое множество отображений .............................66
4.3. Рекуррентные циклы с несколькими использованиями 
рекуррентно вычисляемой переменной..............................................69
4.4. О циклах с нелинейной рекуррентной зависимостью ......................72
4.5. Распараллеливание циклов с рекуррентной зависимостью  
на общей памяти ...................................................................................75
4.6. Векторизация рекуррентных циклов ..................................................82
4.7. Совместное применение распараллеливания и векторизации  
циклов с линейной рекуррентной зависимостью ..............................85

5. Алгоритмы итерационного типа и их ускорение .................................92
5.1. Итерационные алгоритмы решения СЛАУ  
с  блочно-ленточными матрицами .......................................................92
5.2. Умножение и итерационное умножение блочно-ленточной  
матрицы на вектор ................................................................................97
5.3. О сводимости программ, аффинно вычисляющих данные, 
к задачам линейной алгебры ................................................................99
5.4. Операторы сдвига, операторы умножения на функцию 
и соответствующие им матрицы .......................................................107
5.5. Операторы сдвига ...............................................................................109
5.6. Действия операторов сдвига на массивы и примеры 
соответствующих этим действиям матриц .......................................114
5.7. О восстановлении алгоритма типа Гаусса-Зейделя по СЛАУ  
с блочно-ленточной матрицей ...........................................................119

6. Параллельные алгоритмы для вычислительных систем 
с распределенной памятью ....................................................................122
6.1. Блочно-аффинные размещения массивов  
в распределенной памяти ...................................................................122
6.2. Межпроцессорные пересылки данных .............................................126
6.3. Оптимальное количество процессорных элементов 
при параллельном вычислении массивов данных рекуррентным 
программным циклом .........................................................................126
6.4. Изменение ускорения многопроцессорной вычислительной 
системы при повышении быстродействия процессоров ................128

Заключение ....................................................................................................130

Литература ....................................................................................................131

СПИСОК СОКРАЩЕНИЙ

 АЛУ (ALU) — арифметико-логическое устройство.
 
ВС — вычислительная система.
 
ВУ — вычислительное устройство.
 
ВЭ — вычислительный элемент.
 
МВС — многопроцессорная вычислительная система.
 
МП — модуль памяти.
 
ОРС — открытая распараллеливающая система.
 
ПЛИС — программируемая логическая интегральная схема.
 
ПО — программное обеспечение.
 
ПЭ — процессорный элемент.
 
СЛАУ — система линейных алгебраических уравнений.
 
ARM — advanced RISC machines.
 
CPU — central processing unit.
 
DVM — distributed virtual memory system — распараллеливающая 
система ИПМ им. Келдыша.
 
GPU — graphics processing unit.
 
LL — представление нижнетреугольной трехдиагональной матрицы 
в виде произведения двух нижнетреугольных двух-
диагональных матриц.
 
LU — представление квадратной матрицы в виде произведения 
нижнетреугольной и верхнетреугольной матриц.
 
MIMD — multiple instruction, multiple data.
 
MPI — massage passing interface.
 
RISC — reduced instruction set computing.
 
SIMD — single instruction, multiple data.
 
SoC — system-on-chip — система на кристалле.
 
SUIF — Stanford university intermediate format.
 
UU — представление верхнетреугольной трехдиагональной матрицы 
в виде произведения двух верхнетреугольных двух-
диагональных матриц.

ВВЕДЕНИЕ

В данной работе представлены параллельные алгоритмы решения систем 
линейных алгебраических уравнений (СЛАУ) с блочно-ленточными 
матрицами.
Блочно-ленточные матрицы возникают во многих задачах математического 
моделирования при конечноразностных или конечноэлементных численных 
методах. Например, при решении краевых задач в прямоугольной 
области или в трехмерном случае в прямоугольном параллелепипеде с декартовой 
сеткой [1; 2]. Но для более сложных областей при использовании 
ортогональной криволинейной сетки тоже можно прийти к решению СЛАУ 
с блочно-ленточной матрицей [3; 4].
Системы уравнений могут быть большой размерности, иметь много 
блочных диагоналей и решаться итерационными методами [5; 6]. Но итерационные 
алгоритмы решения больших СЛАУ могут в качестве вспомогательных 
содержать задачи решения СЛАУ с меньшим числом блочных 
диагоналей, например трехдиагональные. Быстрым алгоритмам решения 
СЛАУ с трехдиагональными матрицами посвящено много работ [7; 8; 9; 
10; 11; 12; 13]. Метод параллельно-циклической редукции для скалярного 
случая трехдиагональных матриц (к матричному случаю также применим) 
описан в [14]. Условия устойчивой сходимости матричной прогонки описаны 
в [15] и относятся в равной степени и к последовательным алгоритмам, 
и к параллельным.
Оценки сложности давних алгоритмов в большинстве публикаций приводятся 
через количество арифметических операций. Такая оценка сложности 
плохо отражает время работы алгоритма на современных вычислительных 
системах (в первую очередь, процессорах), время работы которых 
более зависит от количества обращений к оперативной памяти, а, еще точнее, 
от движения данных по иерархии памяти.
Алгоритмы, предлагаемые в данной работе, в значительной степени 
ориентированы на оптимизацию использования памяти.
Решение систем уравнений с трехдиагональными матрицами может 
сводиться к вычислению программных циклов с линейной или дробно- 

Алгоритмы решения СЛАУ с блочно-ленточными матрицами. Введение    
7

линейной рекуррентной зависимостью, что тоже допускает параллельное 
выполнение [16; 17; 18].
О параллельных численных методах решения СЛАУ написано много 
книг и обзорных статей, среди которых, в дополнение к предыдущим, можно 
упомянуть [19; 20; 21; 22; 23; 24; 25].
Параллельные алгоритмы решения СЛАУ с ленточными матрицами 
рассматривались еще в те времена, когда слово «распараллелить» было 
почти синонимом слова «ускорить». При этом такие алгоритмы учитывали 
конкретные параллельные архитектуры вычислительных машин, хотя и в 
те времена параллельные архитектуры отличались многообразием.
Иногда использовалось такое понятие, как параллельная сложность алгоритмов, 
т. е. количество одновременно выполняемых операций как функция 
от количества входных данных. При этом не учитывались не только 
операции пересылок и чтения/записи данных, но и количество вычислительных 
устройств [21]. Например, определитель квадратной матрицы n*n 
можно параллельно считать за n*log(n) шагов, используя определение, т. е. 
как сумму n! слагаемых, каждое из которых является произведением n элементов 
матрицы. Зачастую, если виделись в алгоритме операции, которые 
можно выполнить одновременно, то считалось, что их можно выполнить 
за один такт, не беспокоясь о пересылках данных или о времени доступа к 
памяти.
В настоящее время параллельные вычислительные архитектуры стали 
еще более разнообразными и, конечно, для конкретной задачи не может 
быть одного параллельного алгоритма, который был бы одинаково хорош 
для всех. В данной работе рассматривается несколько параллельных алгоритмов 
для решения целевых СЛАУ на нескольких типах параллельных 
компьютеров. Для многоядерных процессоров проведены численные эксперименты 
для нескольких базовых алгоритмов.
По причине многих применений в практически важных задачах СЛАУ 
с блочно-ленточными матрицами исследовались давно, но в последние 
десятилетия появились эффективные на современных (2015–2020 годов) 
процессорах параллельные алгоритмы вычисления циклов с линейной рекуррентной 
зависимостью [16; 17; 18] и эффективные итерационные алгоритмы, 
сочетающие распараллеливание со скошенным (или ромбическим) 
тайлингом [26; 27; 28; 29]. Появляются процессоры (системы на кристалле) 
с более чем 1000 ядер, для которых нужны алгоритмы, использующие соответствующий 
уровень параллелизма.

Б. Я. Штейнберг, О. Б. Штейнберг  

Данная монография отличается направленностью на решение СЛАУ с 
блочно-ленточными матрицами. Анализ сходимости численных методов 
решения СЛАУ описывается, например, в [5; 15; 30].
Для создания быстрых алгоритмов и программ с блочно-ленточными 
матрицами необходимо учитывать не только численные методы линейной 
алгебры, но и особенности вычислительных архитектур и компиляторные 
оптимизирующие преобразования. Необходимая информация по этим вопросам 
приводится в первых главах монографии и носит ознакомительный 
характер. Один из основных результатов монографии — это параллельные 
алгоритмы вычисления программных циклов с рекуррентными 
зависимостями, к которым сводятся некоторые прямые методы решения 
СЛАУ с блочно-ленточными матрицами. Эффективность представленных 
алгоритмов подтверждается результатами численных экспериментов, демонстрирующих 
ускорение в несколько раз. Другой основной результат — 
анализ и сводимость алгоритмов итерационного типа для решения СЛАУ 
с блочно-ленточными матрицами к возможности применения сочетания 
скошенного (ромбического) тайлинга с распараллеливанием. Для метода 
Гаусса-Зейделя решения двумерной задачи Дирихле и некоторых задач математической 
физики этот подход дает ускорение в десятки раз [18; 26; 27; 
28; 29]. В этих же источниках приводится описание метода сочетания скошенного 
тайлинга и распараллеливания, который не приводится в данной 
монографии. Во всех этих случаях матрица СЛАУ блочно-ленточная. Описать 
класс задач, к которым можно эффективно применить сочетание скошенного 
тайлинга и распараллеливания — интересная и полезная задача, 
которая сформулирована и существенно продвинута в данной монографии.
Большинство представленных в данной работе алгоритмов описаны на 
псевдокоде, близком к языку С.

1. ВЛИЯНИЕ ВЫЧИСЛИТЕЛЬНЫХ АРХИТЕКТУР 
И КОМПИЛЯТОРОВ НА БЫСТРОДЕЙСТВИЕ 
АЛГОРИТМОВ И ИХ ПРОГРАММНЫХ 
РЕАЛИЗАЦИЙ

1.1. Разнообразие и особенности 
вычислительных систем

Чем сложнее архитектура, тем в большей степени эффективность высокоуровневой 
программы зависит от компилятора.
Будем рассматривать основные типы вычислительных систем (ВС) и их 
модулей памяти. Многие типы ВС допускают многообразие комбинаций 
(гибридные архитектуры). Выделяют основные типы параллельных вычислительных 
архитектур: SIMD, MIMD, кроме которых есть еще конвейерные, 
систолические, программируемые и др. Но высокая производительность 
существенно зависит еще и от структуры памяти ВС.
Основные типы памяти параллельных вычислительных систем: общая 
и распределенная. Но рассматривается также и гибридная память, и иерархия 
разных видов модулей памяти на микросхеме процессора.
Общая память. Общая память может иметь специфику доступа: многоканальная 
память, ортогональная память, ассоциативная память, иерархия 
модулей памяти, основная или вспомогательная (более быстрая, но с меньшим 
объемом). В микросхемы процессоров встраивается быстрая вспомогательная 
память: разных характеристик кэш-память или адресуемая локальная 
общая память.
Распределенная память и коммуникационные сети. При распределенной 
памяти обмен данными выполняется с помощью коммуникационной 
сети или коммутатора. Часто используются коммуникационные сети: 
шина, кольцо, решетка (mesh), многокольцевая сеть, тор (двумерный и 
трехмерный), гиперкуб, звезда, дерево (см., например: [31]), однородный 
граф малого диаметра (см., например: [32; 33; 34]).

Б. Я. Штейнберг, О. Б. Штейнберг  

При соединении компьютеров сетью «звезда» предполагается, что один 
из компьютеров головной, а каждый из остальных подчиненный, который 
соединяется только с головным.

Рисунок 1. Коммуникационная сеть «звезда».

Иерархическое соединение предполагает, что к головному компьютеру 
подсоединены подчиненные первого уровня, к которым подсоединены компьютеры 
второго уровня, и т. д.

Рисунок 2. Иерархическое (древовидное) соединение.

Звездные и иерархические сети не являются однородными и не обладают 
сильной симметрией. Количество связей звезды и иерархической сети 
равно количеству вершин минус 1.

Рисунок 3. Шина, соединяющая процессорные элементы.

Шина и звезда являются частными случаями дерева.

Алгоритмы решения СЛАУ с блочно-ленточными матрицами. Глава  1 
11

Рисунок 4. Кольцевая сеть.

Кольцевая сеть бывает однонаправленной и двунаправленной: однона-
правленное кольцо позволяет пересылать данные только в одну сторону, а 
двунаправленное — в обе стороны. И та, и другая кольцевая сеть обладает 
симметрией, а количество ребер равно количеству вершин. Двунаправлен-
ная кольцевая сеть использовалась, например, в суперкомпьютере ПС-2000.
Многомерная решетка (mesh) — это коммуникационная сеть, вершины 
которой расположены в узлах многомерной целочисленной декартовой решетки, 
а связи соединяют соседние узлы по координатным направлениям.

Рисунок 5. Двумерная решетка.

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

Б. Я. Штейнберг, О. Б. Штейнберг  

Рисунок 6. Двумерный тор.

Диаметр графа — это наибольшее расстояние между всеми парами вершин. 
От диаметра коммуникационной сети зависит время пересылки данных.
Диаметр звезды равен 2. Однонаправленная кольцевая сеть имеет диаметр, 
равный количеству вершин, а двунаправленная — равный половине 
от количества всех вершин (процессорных элементов). Диаметр 
m-мерной решетки размера N = n1*n2 *…* nm равен сумме размерностей 
(n1 + n2 + … + nm), а m-мерного тора — половине суммы (n1 + n2 +…+ nm) / 2. 
Если все стороны решетки равны, то диаметр m-мерной решетки равен 
m*N 1/m, а диаметр m-мерного тора — m*N 1/m / 2. Нетрудно заметить, что 
одномерная решетка оказывается кольцом.
Чтобы большой диаметр кольцевой шины сократить (и этим самым 
сократить время пересылок данных), можно соединить процессорные элементы 
еще одним кольцом (или несколькими кольцами). Если вершины 
первого кольца занумерованы в порядке их следования (например, по часовой 
стрелке), то, ввиду естественных требований регулярности (симметрии), 
второе кольцо должно соединять вершины, разность номеров между 
которыми постоянна (по модулю количества всех вершин N).

Рисунок 7. Двойное кольцевое соединение.