Физические задачи очень часто сводятся к СЛАУ с разреженными матрицами коэффициентов. Такие матрицы можно преобразовать в ленточные, работать с которыми – приятно и эффективно. Разберемся с различными форматами представления матриц и оценим производительность соответствующих решателей, попутно сделав для них одинаковый интерфейс.
Введение
СЛАУ с сильно разреженными матрицами, например, возникают при решении физических задач методом конечных элементов (МКЭ). При этом матрицы получается также симметрическими и положительно определенными.
Чтобы эффективно решать СЛАУ с сильно разреженными матрицами можно пойти двумя путями:
Использовать специальные библиотеки для разреженных матриц – такие как PARDISO, MUMPS, SPOOLES, TAUCS, SuperLU, UMFPACK, PSPASES, WSMP, PaStiX и другие.
Постараться все ненулевые элементы матриц сгруппировать около главной диагонали. Получится ленточная матрица. Если ширину ленты удастся сделать достаточно малой, то для нее, возможно, более эффективны будут ленточные решатели. Они имеются, например, в той же библиотеке LAPACK и по сравнению с первым вариантом – проще в использовании.
В обоих случаях в компьютерной памяти хранится только «содержательная» часть матрицы. Для варианта 1 – ненулевые компоненты и их «адреса», для варианта 2 – лента.
Получить ленточную матрицу можно тоже двумя путями:
В алгоритме решения базовой задачи позаботиться, чтобы матрица получалась ленточной «автоматом». Для МКЭ это означает, что нужно правильно организовать нумерацию узлов конечных элементов в системе. Это может быть чрезвычайно эффективно для всяких «вытянутых» и «длинных» систем. Скажем, мы на работе обсчитываем контактную подвеску железных дорог. Она – очень-преочень длинная, для нее (если позаботиться о правильной нумерации узлов), ленточная матрица получается естественным образом.
Второй путь – сначала сделать разреженную матрицу «как придется», а затем преобразовать ее в ленточную посредством перенумерации строк и столбцов. Для этого существуют эффективные алгоритмы, такие как GPS, GGPG, RCM, Sloan, JCL и другие. В случае с МКЭ можно автоматически преренумеровывать не столбцы и строки, а номера узлов конечных элементов т.е. номера блоков в блочном представлении матрицы. В одной из задач мы использовали для этого алгоритм JCL (сайт). Результаты превзошли всякие ожидания: программа из типа «Запустил – пойду пообедаю» стала интерактивной – расчет стал выполняться секунды.
Проиллюстрируем перенумерацию узлов алгоритмом JCL для простейшей системы и ее блочной матрицы. Слева – до, справа – после. Цифры в матрице – ненулевые элементы у блоков.
Различных подпрограмм для решения СЛАУ с ленточными матрицами – тоже множество. Для нас во-первых интересно оценить скорость их работы, а во-вторых – реализовать некоторый единый интерфейс для работы с ними, скрывающий различные форматы представления матриц. Это позволит пробовать разные решатели в наших «околонаучных» программах и менять их «на лету».
Почему это важно? Дело в том, что разные решатели отличаются не только скоростью. Они также различаются параметрами устойчивости алгоритмов и имеют разные «пределы терпимости» к плохо обусловленным матрицам, которые иногда получаются в жизни. В некоторых задачах мы сильно намучились с устойчивостью решений. Хочется опробовать разные решатели в реальных программах, чтобы выбрать самые «некапризные» в плане устойчивости и приемлемые по скорости. Сделать это без единого интерфейса – очень проблемно.
Как оценить устойчивость в синтетическом тесте – не очевидно. Устойчивость – скорее характеристика решения реальной задачи, поэтому «тест на устойчивость» в рамках данной заметки мы делать не будем. Ограничимся разбирательством с форматами ленточных матриц, формированием интерфейса решателей и «бытовым» тестом их скоростных качеств.
Все «интерфейсы» и тесты будут сделаны на Фортране с компилятором GFortran. В прошлый раз мы убедились, что в разных языках решатели СЛАУ вызываются одни и те же (в виде скомпилированных библиотек). Значит – переписывать одно и то же на разных языках нет смысла.
Будут рассмотрены решатели и соответствующие форматы матриц из следующих библиотек:
LAPACK (+BLAS) в трех реализациях: классической (версия 3.7.0), OpenBLAS (0.2.19) и Intel MKL (11.3 – он поставился у меня из дистрибутива Python Anaconda).
LAIPE2 от профессора Ляо.
Протестируем основные подпрограммы для решения СЛАУ из этих библиотек, доступные для различных форматов представления матриц и основанные на разложениях \(LU\), \(LL^T\) (\(U^TU\)) и \(LDL^T\) (\(U^TDU\)).
Дополнительно будет упомянут формат ленточных матриц, примененный в старейшей математической библиотеке IMSL. Когда-то, давным-давно :), я начинал изучать Фортран в одном университете, где использовали Microsoft FORTRAN PowerStation, в который входила IMSL. Это было очень крутое решение, мы делали такое, что все завидовали. :) Потом Microsoft FORTRAN приказал долго жить, IMSL, кажется, какое-то время входила «из коробки» в Intel FORTRAN, потом они заменили ее на MKL. Права на IMSL сейчас принадлежат фирме Rogue Wave Software. Это коммерческий продукт. Мне нужно вспомнить формат матриц оттуда, чтобы перевести некоторые старые программы на открытый LAPACK. Кстати, по IMSL есть прекрасная книга на русском языке О.В. Бартеньева (в трех томах).
Форматы представления симметрических матриц
Пусть у нас имеется квадратная симметрическая и (скорее всего) ленточная матрица A. Обозначим:
n – размер матрицы.
nc – число верхних (или нижних) кодиагоналей матрицы, если она ленточная. А если нет, то этот параметр примем по максимальному числу кодиагоналей: nc = n-1
Матрица симметрическая: A(i,j) = A(j,i).
Различные решатели «берут на на вход» матрицу в различном формате ее представления. Выгодно работать с «упакованными» форматами, где нулевые, а также одинаковые в силу симметрии элементы в памяти не хранятся. Формат представления матрицы у различных решателей свой. Основных вариантов (если не брать произвольно разреженные) здесь три: общий вид (будем обозначать его буквой «G»), формат в представлении упакованных треугольников («T») и упакованных лент («B»). С учетом некоторых особенностей к этим буквам будем добавлять и другие: «S» – симметричная и т.д.
Упакованную матрицу будем обозначать P (packed).
G – общий формат
Это самый простой – «обычный» формат. Ради общности рассуждений приведем его на примере условной матрицы 8х8:
A: Столбцы, j=1..n Размер матрицы A: (n, n)
1 2 3 4 5 6 7 8 Обращение к элементу: А(i, j)
┌─────────────────────────
1│ 11 12 13 14 15 16 17 18 n=8, nc=7.
2│ 21 22 23 24 25 26 27 28
3│ 31 32 33 34 35 36 37 38 <– внутри не значения,
Строки, 4│ 41 42 43 44 45 46 47 48 а цифры - номера
i=1..n 5│ 51 52 53 54 55 56 57 58 строк/столбцов (i, j)
6│ 61 62 63 64 65 66 67 68 там, где элементы ненулевые
7│ 71 72 73 74 75 76 77 78
8│ 81 82 83 84 85 86 87 88
На всякий случай, заметим, что в Фортране матрица в памяти хранится по столбцам, т.е. так: 11, 21, 31, 41... 12, 22, 32.... Для оптимизации скорости перебирать ее элементы в циклах надо так, чтобы более часто (во внутреннем цикле) менялся левый индекс (i), а более редко (во внешнем цикле) – правый (j).
ST – упакованный треугольник
Если матрица симметрическая, можно хранить в памяти только уникальные элементы, например, нижний треугольник матрицы. Он упаковывается в профиль, представляющий собой вектор, куда записываются нижние симметричные элементы – последовательно сверху вниз по столбцам:
A: Столбцы, j=1..n Размер матрицы A: (n, n)
1 2 3 4 5 6 7 8 Обращение к элементу: А(i, j)
┌─────────────────────────
1│ 11 n=8, nc=7.
2│ 21 22 Симметрично
3│ 31 32 33
Строки, 4│ 41 42 43 44
i=1..n 5│ 51 52 53 54 55
6│ 61 62 63 64 65 66
7│ 71 72 73 74 75 76 77
8│ 81 82 83 84 85 86 87 88
Профиль P: 1│ 11
2│ 21 Размер профиля P: ((n+1)*n)/2
3│ 31
Элементы, 4│ 41 Обращение к элементу, при j >= i:
ib = ((n+1)*n)/2 5│ 51
6│ 61 A(i, j) = P(i+(2*n-j)*(j-1)/2)
7│ 71
8│_81
9│ 22
10│ 32
11│ 42
12│ 52
13│ 62
14│ 72
15│_82
16│ 33
..│....
36│ 88
С таким форматом работают целый ряд функций библиотеки LAPACK.
STL – упакованный треугольник для LAIPE2
В принципе, это точно такой же формат, как и ST. Единственное отличие – библиотека LAIPE2 при этом дополнительно работает в вектором меток L, который упрощает адресацию к ячейкам матрицы:
A: Столбцы, j=1..n Размер матрицы A: (n, n)
1 2 3 4 5 6 7 8 Обращение к элементу: А(i, j)
┌─────────────────────────
1│ 11 n=8, nc=7.
2│ 21 22 Симметрично
3│ 31 32 33
Строки, 4│ 41 42 43 44
i=1..n 5│ 51 52 53 54 55
6│ 61 62 63 64 65 66
7│ 71 72 73 74 75 76 77
8│ 81 82 83 84 85 86 87 88
Профиль P: 1│ 11 1│ 1 <- Вектор меток L:
2│ 21 2│ 8 Размер: n
3│ 31 3│14 L(1) = 1
Элементы, 4│ 41 4│19 L(k) = L(k-1) + n-i+1, k=2..n
ib = ((n+1)*n)/2 5│ 51 5│23
6│ 61 6│26
7│ 71 7│28
8│_81 8│29
9│ 22
10│ 32 Размер профиля P: ((n+1)*n)/2
11│ 42
12│ 52 Обращение к элементу, при j >= i:
13│ 62
14│ 72 A(i, j) = P(L(j)+i-1)
15│_82
16│ 33
..│....
36│ 88
GB – упакованная ленточная матрица общего вида
Этот формат применяется в основном для несимметрических ленточных матриц в LAPACK. Применительно к симметрической матрице он будет выглядеть так:
A: Столбцы, j=1..n Размер матрицы A: (n, n)
1 2 3 4 5 6 7 8 Обращение к элементу: А(i, j)
┌─────────────────────────
1│ 11 12 13 14 n=8, nc=3.
2│ 21 22 23 24 25
3│ 31 32 33 34 35 36
Строки, 4│ 41 42 43 44 45 46 47
i=1..n 5│ 52 53 54 55 56 57 58
6│ 63 64 65 66 67 68
7│ 74 75 76 77 78
8│ 85 86 87 88
Столбцы, jp=1..n
P: 1 2 3 4 5 6 7 8
┌───────────────────────── Размер матрицы P: (3*nc+1, n)
1│ * * * * * * * *
2│ * * * * * * * *
3│ * * * * * * * * Обращение к элементу:
4│ * * * 14 25 36 47 58
5│ * * 13 24 35 46 57 68 A(i, j) = P(2*nc+1+i-j, j)
6│ * 12 23 34 45 56 67 78
Строки, 7│ 11 22 33 44 55 66 77 88
ip = 1..2*nc+1 8│ 21 32 43 54 65 76 87 *
9│ 31 42 53 64 75 86 * * * – не используется
10│ 41 52 63 74 85 * * *
LAPACK почему-то требует, чтобы у матрицы P было nc первых пустых (несодержательных) строк. Возможно, имея их для вспомогательных нужно можно ускорить какие-то алгоритмы.
SB – упакованная симметрическая ленточная матрица
Формат упаковки симметрических матриц предусматривает хранение только, например, нижних относительно диагонали элементов:
A: Столбцы, j=1..n Размер матрицы A: (n, n)
1 2 3 4 5 6 7 8 Обращение к элементу: А(i, j)
┌─────────────────────────
1│ 11 n=8, nc=3.
2│ 21 22 Симметрично
3│ 31 32 33
Строки, 4│ 41 42 43 44
i=1..n 5│ 52 53 54 55
6│ 63 64 65 66
7│ 74 75 76 77
8│ 85 86 87 88
Столбцы, jp=1..n Размер матрицы P: nc+1, n
P: 1 2 3 4 5 6 7 8
┌───────────────────────── Обращение к элементу, при i>=j:
1│ 11 22 33 44 55 66 77 88
Строки, 2│ 21 32 43 54 65 76 87 * A(i, j) = P(1+i-j, j)
ip = 1..nc+1 3│ 31 42 53 64 75 86 * *
4│ 41 52 63 74 85 * * * * – не используется
SBL – упакованная симметрическая ленточная матрица для LAIPE2
В LAIPE2 симметрическая ленточная матрица с постоянной шириной ленты упаковывается в профиль, представляющий собой вектор, куда записываются нижние симметричные элементы последовательно по столбцам:
Матрица A: Столбцы, j=1..n Размер матрицы A: (n, n)
1 2 3 4 5 6 7 8 Обращение к элементу: А(i, j)
┌─────────────────────────
1│ 11 n=8, nc=3.
2│ 21 22 Симметрично
3│ 31 32 33
Строки, 4│ 41 42 43 44
i=1..n 5│ 52 53 54 55
6│ 63 64 65 66
7│ 74 75 76 77
8│ 85 86 87 88
* * * – доп. элементы для записи профиля
*
Профиль P: 1│ 11
2│ 21 Размер профиля: (n-1)*nc+n
3│ 31
Элементы, 4│_41
ip = 1..(n-1)*nc+n 5│ 22
6│ 32
7│ 42 Обращение к элементу, при i >= j:
8│_52
9│ 33 A(i, j) = P((nc+1)*(j-1)+i-j+1)
10│ 43
11│ 53
12│_63
13│ 44
14│ 54
15│ 64
16│_74
17│ 55
18│ 65
19│ 75
20│_85
21│ 66
22│ 76
23│ 86
24│_ *
25│ 77
26│ 87
27│ *
28│_ *
29│ 88
SBLV – упакованная симметрическая ленточная матрица с переменной шириной ленты для LAIPE2
Ленточные матрицы часто получаются с переменной шириной ленты. Так зачем хранить в памяти ленту по ее максимальной ширине, если можно дополнительно откинуть также нулевые элементы в местах, где лента сужается?
Чтобы пояснить формат для таких матриц, «выгрызем» из исходной матрицы А еще несколько элементов. Такая матрица упаковывается в профиль, представляющий собой вектор, куда записываются верхние симметричные элементы последовательно по столбцам (сверху вниз), а информация о ширине ленты в каждой позиции хранится в дополнительном векторе меток L:
A: Столбцы, j=1..n Размер матрицы A: (n, n)
1 2 3 4 5 6 7 8 Обращение к элементу: А(i, j)
┌─────────────────────────
1│ 11 12 n=8, nc=3 (макс. ширина ленты)
2│ 22 23
3│ 33 36
Строки, 4│ 44 45 46
i=1..n 5│ 55 56 57
6│ 66 67
7│ Симметрично 77
8│ 88
Профиль P: 1│_11 1│ 1 <- Вектор меток L:
2│ 12 2│ 2 Размер: n
3│_22 3│ 3 L(1) = 1
Элементы, 4│ 23 4│ 3 L(k) = L(k-1) + nz(k), k=2..n,
ib = 1..sum(nz)+n 5│_33 5│ 4 где nz(k) - число ненулевых элементов
6│_44 6│ 7 над диагональю в k-том столбце
7│ 45 7│ 9
8│_55 8│ 9
9│ 36
10│ 46 Размер профиля: sum(nz)+n
11│ 56
12│_66 Обращение к элементу, при i >= j:
13│ 57 Если j < L(max(1,i-1)) - L(i) + i,
14│ 67 то A(i, j) = 0,
15│_77 иначе A(i, j) = P(L(i)+j-1).
16│_88
SBI – упакованная симметрическая ленточная матрица для IMSL
Ленточная матрица для самых быстрых процедур IMSL хранилась в оригинальном формате. В единую матрицу упаковывается как лента матрицы, так и вектор правых частей. Утверждается, что такая форма упаковки позволяет реализовать особенно эффективные алгоритмы:
A: Столбцы, j=1..n Размер матрицы A: (n, n)
1 2 3 4 5 6 7 8 Обращение к элементу: А(i, j)
┌─────────────────────────
1│ 11 n=8, nc=3.
2│ 21 22 Симметрично
3│ 31 32 33
Строки, 4│ 41 42 43 44
i=1..n 5│ 52 53 54 55
6│ 63 64 65 66
7│ 74 75 76 77
8│ 85 86 87 88
Столбцы, jp=1..nc+2 Размер матрицы P: (n+nc, nc+2)
1 2 3 4 5
P: ┌──────────────── Обращение к элементу, при i>=j:
1│ * * * * *
2│ * * * * * A(i, j) = P(nc+i, i-j+1)
3│ * * * * *
Строки, 4│ 11 * * * b1
ip=1..n+nc 5│ 22 12 * * b2 * – неиспользуемые и вспомог. ячейки
6│ 33 23 13 * b3
7│ 44 34 24 14 b4 b1:b8 – вектор правых частей
8│ 55 45 35 25 b5 и результатов решения
9│ 66 56 46 36 b6
10│ 77 67 57 47 b7
11│ 88 78 68 58 b8
Интерфейс решателей для симметрических матриц
Описание
Напишем маленький отдельный модуль SlaeSolSym для работы с симметрическими матрицами в различных форматах. Надо, чтобы процедуры модуля скрывали от «пользователя» особенности конкретной упаковки матриц, векторы меток LAIPE, а также различные тонкости применения разных решателей. «Пользователь» должен получить единый интерфейс для работы с упакованными матрицами всех видов, как будто он работает с обычными квадратными матрицами. Тогда решатели для разных задач можно будет переключать «на лету».
Общая последовательность работы с модулем будет такая:
Шаг 1. Вызываем процедуру, которая по размеру исходной матрицы, ширине ленты и заданному формату упаковки выделяет память под упакованную матрицу P, вектор правых частей B, а также делает некоторые другие начальные приготовления:
call SlaeBeginAlloc(P, B, n, nc, PackForm, nip, njp, Threads)
P | dimension(:,:) | Вход / выход | Матрица коэффициентов – возможно, что в упакованном формате, на входе – еще не размещенная в памяти, на выходе – размещенная |
B | dimension(:) | Вход / выход | Вектор правых частей, на входе еще не размещенный в памяти, на выходе – размещенный |
n | integer | Вход | Размер матрицы коэффициентов (как для квадратного варианта) |
nc | integer | Вход | Число верхних или нижних кодиагоналей ленты матрицы. Если она не ленточная – установить nc = n - 1 |
PackForm | character(*) | Вход | Формат упаковки матрицы: G, ST, STL, GB, SB, SBL, SBLV или SBI |
nip | integer | Выход | Число строк упакованной матрицы (для справки) |
njp | integer | Выход | Число столбцов упакованной матрицы |
Threads | integer | Вход | Число потоков для библиотек OpenBLAS, MKL, LAIPE2, которые будут им установлены. Для остальных – задавать необязательно. |
(Элементы матриц и векторов – real(8))
Шаг 2. Заполняем матрицу P, работая с ней как с обычной квадратной. Для этого имеется процедура, помещающая элемент в матрицу на место (i, j), соответствующее его адресу в квадратной матрице, вне зависимости от способа упаковки:
call PutToMatrix(P, n, nc, i, j, a, PackForm)
P | dimension(:,:) | Вход / выход | Матрица коэффициентов в упакованном формате |
n | integer | Вход | Размер матрицы коэффициентов при квадратном варианте |
nc | integer | Вход | Число верхних или нижних кодиагоналей ленты матрицы. Если она не ленточная – установить nc = n - 1 |
i | integer | Вход | Номер строки элемента в исходной квадратной матрице |
j | integer | Вход | Номер столбца элемента в исходной квадратной матрице |
a | real(8) | Вход | Помещаемый элемент |
PackForm | character(*) | Вход | Формат упаковки: G, ST, STL, GB, SB, SBL, SBLV или SBI |
Заполнять можно только одну сторону симметрической матрицы (любую).
Если нужно достать элемент из матрицы, предусмотрена аналогичная функция:
GetFromMatrix(P, n, nc, i, j, PackForm) – достает элемент (real(8)) из упакованной матрицы.
Шаг 3. Решаем СЛАУ:
call SlaeSolution(P, B, n, nc, Solver)
P | dimension(:,:) | Вход / выход | Матрица коэффициентов в упакованном формате. На выходе – не сохраняется. |
B | dimension(:) | Вход / выход | Вход: вектор правых частей. Выход – результат решения. |
n | integer | Вход | Размер матрицы коэффициентов при квадратном варианте |
nc | integer | Вход | Число верхних или нижних кодиагоналей ленты матрицы. Если она не ленточная – установить nc = n |
Solver | character(*) | Вход | Название процедуры решателя. Из библиотеки LAPACK: dgesv, dsysv, dposv, dspsv, dppsv, dgbsv или dpbsv. Из библиотеки LAIPE2: DAG_8, DSG_8, DSP_8, CSG_8, CSP_8, VSG_8 или VSP_8. |
Шаг 4. Освобождаем память и делаем завершающие установки:
call SlaeEndDealloc(P, B)
P | dimension(:,:) | Вход / выход | Матрица коэффициентов. На входе – размещенная в памяти, на выходе – освобожденная. |
B | dimension(:) | Вход / выход | Вектор правых частей или решений. На входе – размещенный в памяти, на выходе – освобожденный. |
Особенности:
Между вызовами SlaeBeginAlloc и SlaeEndDealloc можно работать только с одними экземплярами матриц (в модуле имеются внутренние данные, имеющиеся в единственном числе).
Формат SBI поддерживается только на уровне выделения памяти, упаковки и адресации к элементам.
Какие библиотеки применяются с модулем (классический LAPACK, OpenBLAS, MKL, LAIPE2) – управляется директивами компилятора.
В модулей введен тип данных, определяющий вариант решателя СЛАУ и имеется вспомогательная процедура его установки по условному номеру SetSlaeSolution.
При работе с форматом ленты переменной ширины SBLV процедура SlaeBeginAlloc выделяет память под ленту с максимальной шириной nc. После заполнения матрицы нужно вызвать дополнительную процедуру RepackSBLV, которая перепакует матрицу P на основе фактических значений ее элементов.
Код модуля
!> Модуль с интерфейсом для различных вариантов решения СЛАУ с симм. матрицей коэфф. real(8)
module SlaeSolSym
! Включить при любой реализации LAPACK (в том числе OpenBLAS/MKL):
#define LAPACK
! Включить если используются данные библиотеки:
#define OpenBLAS
!#define MKL
#define LAIPE2
implicit none
!> Тип, определяющий вариант решателя СЛАУ
type SlaeSolType
character(6) :: Library ! Библиотека: LAPACK или LAIPE2
character(4) :: PackForm ! Формат представления симметрической матрицы:
! G - общий, ST - упак.треуг., STL - упак.треуг.LAIPE
! GB - общ.лент, SB-упак.симм.лент., SBL-упак.симм.лент.LAIPE
! SBLV-упак.симм.лент.перем.шир.LAIPE, SBI-упак.симм.лент.IMSL
character(4) :: FactType ! Тип факторизации (справочно): LU, LLt, LDLt (или аналог.)
character(5) :: Solver ! Процедура решения -
! LAPACK: dgesv, dsysv, dposv, dspsv, dppsv, dgbsv, dpbsv.
! LAIPE2: DAG_8, DSG_8, DSP_8, CSG_8, CSP_8, VSG_8, VSP_8.
endtype
!> Внутренние вспомогательные данные модуля
integer, dimension(:), allocatable, private :: ipiv ! Вспомогательные матрицы
real(8), dimension(:), allocatable, private :: work ! для некоторых процедур LAPACK
integer, private :: lwork ! Вспомогат. параметр для LAPACK
integer, dimension(:), allocatable, private :: Lbls ! Вектор меток процедур LAIPE2
contains
!> Выделяет память для упакованной матрицы P и вектора правых частей B
!> Для форматов STL и SBLV, также выделяет память под вектор меток и заполняет его
!> Для формата SBLV считает ширину ленты пока постоянной
!> Делает нач.приготовления: для OpenBLAS, MKL, LAIPE2 устанавливает число потоков
subroutine SlaeBeginAlloc(P, B, n, nc, PackForm, nip, njp, Threads)
real(8), dimension(:,:), allocatable :: P !> Упакованная матрица (вх/вых)
real(8), dimension(:), allocatable :: B !> Вектор правых частей (вх/вых)
integer :: n !> Размер матрицы
integer :: nc!> Число нижн.кодиагоналей ленты (вх)
character(*) :: PackForm !> Тип формата матрицы (вх)
integer :: nip, njp !> Размеры упак.матрицы (вых)
integer, optional :: Threads !> Число потоков (опц.)
integer :: k
!---------------------------------------------------------------------------------
select case (trim(PackForm))
case ("G") ! --- Матрица общего вида ---
nip = n; njp = n
case ("ST") ! --- Упакованный треугольник ---
nip = ((n+1)*n)/2; njp = 1
case ("STL") ! --- Упакованный треугольник LAIPE2 ---
nip = ((n+1)*n)/2; njp = 1
allocate(Lbls(n))
Lbls(1) = 1 ! Формируем массив меток
do k = 2, n
Lbls(k) = Lbls(k-1) + n-k+1
enddo
case ("GB") ! --- Упакованная лента общего вида ---
nip = 3*nc+1; njp = n
case ("SB") ! --- Упакованная симметричная лента ---
nip = nc+1; njp = n
case ("SBL") ! --- Упакованная симметричная лента LAIPE2 ---
nip = (n-1)*nc+n; njp = 1
case ("SBLV") ! --- Упак. симм. лента перем. ширины LAIPE2 ---
nip = (n-1)*nc+n; njp = 1
allocate(Lbls(n))
Lbls(1) = 1 ! Формируем массив меток пока считая
do k = 2, n ! ширину ленты фактически постоянной
Lbls(k) = Lbls(k-1) + min(k-1, nc)
enddo
case ("SBI") ! --- Упакованная симметричная лента IMSL ---
nip = n+nc; njp = nc+2
case default
stop ("Ошибка. SlaeBeginAlloc: недопустимый формат упаковки.")
endselect
allocate(P(nip,njp), B(n))
#ifdef OpenBLAS
call OpenBLAS_set_num_threads(Threads)
#endif
#ifdef MKL
call MKL_set_num_threads(Threads)
#endif
#ifdef LAIPE2
call laipe$use(Threads)
#endif
endsubroutine
!> Освобождает память из упакованной матрицы P, вектора правых частей B,
!> и вспомогательных переменных модуля. Завершает работы библиотек
subroutine SlaeEndDealloc(P, B)
real(8), dimension(:,:), allocatable :: P !> Упакованная матрица (вх/вых)
real(8), dimension(:), allocatable :: B !> Вектор правых частей (вх/вых)
!-----------------------------------------------------------------------------
deallocate(P, B)
if (allocated(Lbls)) deallocate(Lbls)
if (allocated(work)) deallocate(work)
if (allocated(ipiv)) deallocate(ipiv)
#ifdef LAIPE2
call laipe$done()
#endif
endsubroutine
!> Засылает элемент в упакованную матрицу для места (i,j) обычной
subroutine PutToMatrix(P, n, nc, i, j, a, PackForm)
real(8), dimension(:,:) :: P !> Упакованная матрица (вх)
integer :: n !> Размер матрицы (вх)
integer :: nc !> Число нижн.кодиагоналей ленты (вх)
integer :: i, j !> Строка/столбец обычной формы (вх)
real(8) :: a !> Засылаемый элемент (вх)
character(*) :: PackForm !> Тип формата матрицы (вх)
integer :: i_, j_
!---------------------------------------------------------------------------
if (i>=j) then
i_ = i; j_ = j
else
i_ = j; j_ = i
endif
if (i_>j_+nc) stop ("Ошибка. PutToMatrix: выход за границы ленты.")
select case (trim(PackForm))
case ("G") ! --- Матрица общего вида ---
P(i_,j_) = a; P(j_,i_) = a
case ("ST") ! --- Упакованный треугольник ---
P(i_+(2*n-j_)*(j_-1)/2, 1) = a
case ("STL") ! --- Упакованный треугольник LAIPE2 ---
P(Lbls(j_)+i_-1, 1) = a
case ("GB") ! --- Упакованная лента общего вида ---
P(2*nc+i_-j_+1, j_) = a; P(2*nc+j_-i_+1, i_) = a
case ("SB") ! --- Упакованная симметричная лента ---
P(i_-j_+1, j_) = a
case ("SBL") ! --- Упакованная симметричная лента LAIPE2 ---
P((nc+1)*(j_-1)+i_-j_+1, 1) = a
case ("SBLV") ! --- Упак. симм. лента перем. ширины LAIPE2 ---
P(Lbls(i_)+j_-1, 1) = a
case ("SBI") ! --- Упакованная симметричная лента IMSL ---
P(nc+i_, i_-j_+1) = a
case default
stop ("Ошибка. PutToMatrix: недопустимый формат упаковки.")
endselect
endsubroutine
!> Забирает элемент из упакованной матрицы из мест (i,j) обычной
function GetFromMatrix(P, n, nc, i, j, PackForm)
real(8), dimension(:,:) :: P !> Упакованная матрица (вх)
integer :: n !> Размер матрицы
integer :: nc !> Число нижн.кодиагоналей ленты (вх)
integer :: i, j !> Строка/столбец обычной формы (вх)
character(*) :: PackForm !> Тип формата матрицы (вх)
real(8) :: GetFromMatrix !> Берущийся элемент (вых)
integer :: i_, j_ ! Индексы, где i_>=j_
!-----------------------------------------------------------------------------
if (i>=j) then
i_ = i; j_ = j
else
i_ = j; j_ = i
endif
if (i_>j_+nc) then
GetFromMatrix = 0.
else
select case (trim(PackForm))
case ("G") ! --- Матрица общего вида ---
GetFromMatrix = P(i_, j_)
case ("ST") ! --- Упакованный треугольник ---
GetFromMatrix = P(i_+(2*n-j_)*(j_-1)/2, 1)
case ("STL") ! --- Упакованный треугольник LAIPE2 ---
GetFromMatrix = P(Lbls(j_)+i_-1, 1)
case ("GB") ! --- Упакованная лента общего вида ---
GetFromMatrix = P(2*nc+i_-j_+1, j_)
case ("SB") ! --- Упакованная симметричная лента ---
GetFromMatrix = P(i_-j_+1, j_)
case ("SBL") ! --- Упакованная симметричная лента ---
GetFromMatrix = P((nc+1)*(j_-1)+i_-j_+1,1)
case ("SBLV") ! --- Упак.симм.лента пер.шир.LAIPE2 ---
if (j_ < Lbls(max(1,i_-1))-Lbls(i_)+i_) then
GetFromMatrix = 0.
else
GetFromMatrix = P(Lbls(i_)+j_-1, 1)
endif
case ("SBI") ! --- Упакованная симм. лента IMSL ---
GetFromMatrix = P(nc+i_, i_-j_+1)
case default
stop ("Ошибка. GetFromMatrix: недопустимый формат упаковки.")
endselect
endif
endfunction
!> Перепаковывает матрицу в формате SBLV из де факто постоянной ширины ленты
!> в матрицу с реально переменной шириной ленты на основе значений ее эл-тов
subroutine RepackSBLV(P, n, nc, nip, njp)
real(8), dimension(:,:),allocatable :: P !> Упакованная матрица (вх/вых)
integer :: n !> Размер матрицы (вх)
integer :: nc !> Макс.ч. нижн.кодиаг.ленты (вх)
integer :: nip, njp !> Новые размеры профиля (вых)
integer, dimension(:), allocatable :: Lnew ! Новый временный вектор меток
real(8), dimension(:,:),allocatable :: Pnew ! Новый временный профиль
integer :: i, j, ip
real(8) :: a
!-----------------------------------------------------------------------------
allocate(Lnew(n)) ! Формируем новый массив меток...
ip = 0; Lnew(1) = 1
do j = 1, n
do i = max(1,j-nc), j
a = GetFromMatrix(P, n, nc, i, j, "SBLV")
if (abs(a) < 1e-20 .and. i/=j) cycle ! если элемент – почти нулевой
if (j>1) Lnew(j) = Lnew(j-1) + j-i
ip = ip + j-i+1
exit
enddo
enddo
nip = ip; njp = 1
allocate(Pnew(nip, njp)) ! Заполняем новый профиль... лента
ip = 0
do j = 1, n
do i = j-(Lnew(j)-Lnew(max(1,j-1))), j
ip = ip + 1
Pnew(ip,1) = GetFromMatrix(P, n, nc, i, j, "SBLV")
enddo
enddo
deallocate(P); allocate(P(nip, njp)); P=Pnew; Lbls=Lnew
deallocate(Pnew, Lnew)
endsubroutine
!> Решает систему линейных уравнений
subroutine SlaeSolution(P, B, n, nc, Solver)
real(8), dimension(:,:) :: P !> Упакованная матрица коэфф. (вх/вых)
real(8), dimension(:) :: B !> Вектор правых частей / решение (вх/вых)
integer :: n !> Размер матрицы (неупакованной)
integer :: nc !> Число нижн.кодиагоналей ленты (вх)
character(*) :: Solver !> Решатель (см. case в теле)
integer :: info
!-------------------------------------------------------------------------
select case (Solver)
#ifdef LAPACK
case ("dgesv")
if (.not. allocated(ipiv)) allocate(ipiv(n))
call dgesv(n, 1, P, n, ipiv, B, n, info)
case ("dsysv")
if (.not. allocated(ipiv)) allocate(ipiv(n))
if (.not. allocated(work)) then
allocate(work(1))
call dsysv("L", n, 1, P, n, ipiv, B, n, work, -1, info)
lwork=Int(work(1))
deallocate(work); allocate(work(lwork))
endif
call dsysv("L", n, 1, P, n, ipiv, B, n, work, lwork, info)
case ("dposv")
call dposv("L", n, 1, P, n, B, n, info)
case ("dspsv")
if (.not. allocated(ipiv)) allocate(ipiv(n))
call dspsv ("L",n, 1, P, ipiv, B, n, info)
case ("dppsv")
call dppsv("L", n, 1, P, B, n, info )
case ("dgbsv")
if (.not. allocated(ipiv)) allocate(ipiv(n))
call dgbsv(n, nc, nc, 1, P, 3*nc+1, ipiv, B, n, info)
case ("dpbsv")
call dpbsv("L", n, nc, 1, P, nc+1, B, n, info)
#endif
#ifdef LAIPE2
case ("DAG_8")
call laipe$Solution_DAG_8(P, n, B, info)
case ("DSG_8")
call laipe$Solution_DSG_8(P, n, Lbls, B, info)
case ("DSP_8")
call laipe$Solution_DSP_8(P, n, Lbls, B, info)
case ("CSG_8")
call laipe$Solution_CSG_8(P, n, nc, B, info)
case ("CSP_8")
call laipe$solution_CSP_8(P, n, nc, B, info)
case ("VSG_8")
call laipe$Solution_VSG_8(P, n, Lbls, B, info)
case ("VSP_8")
call laipe$Solution_VSP_8(P, n, Lbls, B, info)
#endif
case default
stop ("Ошибка. SlaeSolution: недопустимый тип решателя.")
endselect
if (info/=0) stop ("Ошибка. SlaeSolution: не могу решить СЛАУ.")
endsubroutine
!> Устанавливает тип варианта решателя по условному номеру
subroutine SetSlaeSolution(is, S)
integer :: is !> Номер варианта решателя (вх)
type(SlaeSolType) :: S !> Тип решателя (вых)
select case (is)
case (1)
S%Library="LAPACK"; S%Solver="dgesv"; S%PackForm="G"; S%FactType="LU"
case (2)
S%Library="LAPACK"; S%Solver="dsysv"; S%PackForm="G"; S%FactType="LDLt"
case (3)
S%Library="LAPACK"; S%Solver="dposv"; S%PackForm="G"; S%FactType="LLt"
case (4)
S%Library="LAPACK"; S%Solver="dspsv"; S%PackForm="ST"; S%FactType="LDLt"
case (5)
S%Library="LAPACK"; S%Solver="dppsv"; S%PackForm="ST"; S%FactType="LLt"
case (6)
S%Library="LAPACK"; S%Solver="dgbsv"; S%PackForm="GB"; S%FactType="LU"
case (7)
S%Library="LAPACK"; S%Solver="dpbsv"; S%PackForm="SB"; S%FactType="LLt"
case (8)
S%Library="LAIPE2"; S%Solver="DAG_8"; S%PackForm="G"; S%FactType="LU"
case (9)
S%Library="LAIPE2"; S%Solver="DSG_8"; S%PackForm="STL"; S%FactType="LDLt"
case (10)
S%Library="LAIPE2"; S%Solver="DSP_8"; S%PackForm="STL"; S%FactType="LLt"
case (11)
S%Library="LAIPE2"; S%Solver="CSG_8"; S%PackForm="SBL"; S%FactType="LDLt"
case (12)
S%Library="LAIPE2"; S%Solver="CSP_8"; S%PackForm="SBL"; S%FactType="LLt"
case (13)
S%Library="LAIPE2"; S%Solver="VSG_8"; S%PackForm="SBLV";S%FactType="LDLt"
case (14)
S%Library="LAIPE2"; S%Solver="VSP_8"; S%PackForm="SBLV";S%FactType="LLt"
case default
stop ("Ошибка. SetSlaeSolution: недопустимый тип решателя.")
endselect
endsubroutine
endmodule
Оценка производительности решателей
Тестовые задачи
Мы сделали инструмент, позволяющий быстро выполнить сразу много расчетов. Решим разными решателями следующие 8 задач. В качестве самой первой повторим задачу, описанную в прошлом посте.
| 1 | 1000 | 999 | 500 500 (50%) | Плотная (неленточная) матрица |
| 2 | ––//–– | 250 | 219 625 (22%) | |
| 3 | ––//–– | Переменная: 250/25 | 28 451 (2.9%) | 10% ленты имеют полуширину 250, остальные 90% – 25 |
| 4 | ––//–– | 25 | 25 675 (2.6%) | Узкая |
| 5 | 10 000 | 250 | 2 478 625 (2.5%) | |
| 6 | ––//–– | Переменная: 250/25 | 45 3625 (0.45%) | 10% ленты имеют полуширину 250, остальные 90% – 25 |
| 7 | ––//–– | 25 | 259 675 (0.26%) | Узкая |
| 8 | 100 000 | 5 | 599 985 (0.006%) | Очень узкая |
* – число ненулевых уникальных элементов.
(Уточнение по условию задачи: значение диагональных элементов начиная с шага 2 увеличивалось до 100, с шага 5 – до 200).
В каждом случае СЛАУ решаем 100 раз, замеряем время этих 100 повторений.
Каждую задачу (из 100 шагов) будем запускать по 10 раз, запоминая наименьшее время работы текущего решателя.
Во всех программах будем использовать векторы и матрицы из вещественных чисел двойной точности (занимающие в памяти 8 байт).
Правильность и одинаковость решения разными средствами проверяем по контрольной сумме.
Время решения будем измерять по циклу решения СЛАУ – без учета предварительной подготовки матрицы коэффициентов.
Решение
С учетом нашего единого интерфейса для разных решателей и форматов упаковки матриц разого вида, решение на Фортране имеет почти «человеческое лицо». Здесь – решение первой задачи, для остальных нужно менять некоторые параметры:
program BandTest
use SlaeSolSym
implicit none
integer, parameter :: n = 1000 ! Число уравнений
integer, parameter :: nc = n-1 ! Число кодиагоналей ленты
integer, parameter :: nc2= nc ! Число кодиагоналей после 10% заполнения
integer, parameter :: m = 100 ! Число повторений решения
integer, parameter :: nr = 10 ! Число запусков каждого решателя
integer, parameter :: Threads = 8 ! Число потоков многопоточных решателей
integer :: nip, njp ! Размер упакованной матрицы
real(8),dimension(:,:),allocatable :: A, A0 ! Матрица коэффициентов (упакованная)
real(8),dimension(:), allocatable :: B ! Вектор правых частей / результатов
type(SlaeSolType), dimension(1:14) :: SS ! Массив типов решателей (от:до)
real(8) :: csum ! Контрольная сумма
integer :: i, j, k, is, ir, ncc, ne ! Вспомогательные
real(8) :: aij
integer :: t_0, t_1, t_rate, t_max ! Для измер. времени
real :: t_ss_min
!--------------------------------------------------------------------------------
print "(a3,a8,2a6,2a7,a15)","№","Библ.","Упак.","Разл.","Проц.","Время","К.сумма"
do is=lbound(SS,1), ubound(SS,1) ! Цикл по типам решателей
call SetSlaeSolution(is, SS(is)) ! Установка типа решателя по номеру
!if (SS(is)%PackForm(2:2) /= "B") cycle ! Раскомм.,если надо только лент. (!)
call SlaeBeginAlloc(A0, B, n, nc, SS(is)%PackForm, nip, njp, Threads)
A0=0.; B=0.; ncc=nc; ne=0
do j=1, n ! Заполняем матрицу А0
do i=1, j
if (j<=ncc+i .and. i<=ncc+j ) then
aij = 1-(dble(i-j)/n)**2; if (i==j) aij=10.
call PutToMatrix(A0, n, nc, i, j, aij, SS(is)%PackForm)
ne=ne+1
endif
if (j>n/10) ncc=nc2 ! После 10% меняем ширину ленты на nc2
enddo
enddo
if (SS(is)%PackForm=="SBLV") call RepackSBLV(A0, n, nc, nip, njp)
allocate(A(nip,njp))
t_ss_min = 1.e10
do ir=1, nr ! Цикл с запусками задачи на одном решателе
csum = 0.
call system_clock(t_0, t_rate, t_max)
do k = 1, m ! Цикл повторения решения
A = A0; B = k
call SlaeSolution(A, B, n, nc, SS(is)%Solver) ! Решение СЛАУ
csum = csum + sum(B)
enddo
call system_clock(t_1, t_rate, t_max)
t_ss_min = min(t_ss_min, real(t_1-t_0)/t_rate)
enddo
call SlaeEndDealloc(A0, B) ! Освобождение памяти и конечные установки
deallocate(A)
print "(i3,a8,a6,a6,a7,f7.2,f15.5)", is, SS(is)%Library, SS(is)%PackForm, &
SS(is)%FactType, SS(is)%Solver, t_ss_min, csum
enddo
print '("Ур-ний:",i8," Ненул.уник.эл-тов:",i8," (",f5.2,"%)")',n,ne,ne*100./(n*n)
endprogram
Компилируем (в данном случае для варианта c OpenBLAS и LAIPE2):
gfortran -c SlaeSolSym.F90 BandTest.F90 -static -fdollar-ok -O2 \
-lopenblas -llaipe2 -lneuloop4
gfortran *.o -static -fdollar-ok -O2 \
-lopenblas -llaipe2 -lneuloop4 -o BandTest.exe
Получаем:
№ Библ. Упак. Разл. Проц. Время К.сумма
1 LAPACK G LU dgesv 0.71 22439.19509
2 LAPACK G LDLt dsysv 2.51 22439.19509
3 LAPACK G LLt dposv 0.48 22439.19509
4 LAPACK ST LDLt dspsv 2.52 22439.19509
5 LAPACK ST LLt dppsv 2.47 22439.19509
6 LAPACK GB LU dgbsv 1.43 22439.19509
7 LAPACK SB LLt dpbsv 0.49 22439.19509
8 LAIPE2 G LU DAG_8 7.32 22439.19509
9 LAIPE2 STL LDLt DSG_8 5.14 22439.19509
10 LAIPE2 STL LLt DSP_8 4.95 22439.19509
11 LAIPE2 SBL LDLt CSG_8 4.78 22439.19509
12 LAIPE2 SBL LLt CSP_8 4.52 22439.19509
13 LAIPE2 SBLV LDLt VSG_8 4.18 22439.19509
14 LAIPE2 SBLV LLt VSP_8 4.11 22439.19509
Сравнительные результаты
Тесты проводились на обычном современном персональном компьютере с 4-ядерным 8-поточным процессором Intel i7 с частотой 4 ГГЦ с достаточным размером памяти.
Библиотеки, которые это позволяют (OpenBLAS, MKL и LAIPE2), тестировались в вариантах с настройками на 8, 4 и 1 потоков.
Вест софт, использовался в 64-разрядной реализации. ОС: Windows 10, 64.
n = 1000. Не лента (nc = 999). НУЭ = 500 500 (50%). 100 повторений
Что видно здесь интересного:
Колоссальная разница по времени между самым быстрым и самым медленным решателем (порядка 100 раз).
Треугольное представление матриц в основном не дает выигрыша во времени, скорее наоборот. Выигрыш – только по памяти.
В OpenBLAS есть хорошо оптимизированные решатели, а есть – откровенно плохо (LDLt).
MKL в этом смысле более ровный. Все в нем оптимизировано дальше некуда. Он – очень крут.
В целом даже если матрица не ленточная, ленточные решатели во многих случаях такие же быстрые как и решатели для обычного представления матрицы.
n = 1000. Лента (nc = 250). НУЭ = 219 625 (22%). 100 повторений
Это та же система, только матрица сделана ленточной. На диаграмме оставляем только ленточные решатели (время у остальных мало изменилось – с «их точки зрения» размер задачи остался прежним).
Самые быстрые решатели улучшили результат в 3-4 раза по сравнению с плотной матрицей.
MKL снова впереди всех.
Уменьшился эффект от распараллеливания.
n = 1000. Лента переменной ширины (nc = 250/25 для 10/90%). НУЭ = 28 451 (2.85%). 100 повторений
Для таких матриц форма представления с переменной шириной – эффективна. Благодаря этому в целом медлительный LAIPE2 практически догнал лидера MKL. При этом он быстр лишь на одном потоке, а многопоточность его замедляет. Как и OpenBLAS.
Классический LAPACK – в аутсайдерах.
n = 1000. Узкая лента (nc = 25). НУЭ = 25 675 (2.57%). 100 повторений
OpenBLAS для SB/LLt на 4 и 8 потоках явно неестественно тормозит. Баг? Он же в однопоточном варианте – на порядок быстрее.
Эффект от распараллеливания стал обратный: большинство решателей лучше работают на 1 потоке. Остальные – примерно одинаково что на одном, что на восьми.
Самое удивительное: ранее тормозной классический LAPACK (SB/LLt) вырывается вперед и наступает на пятки самому MKL!
n = 10 000. Лента (nc = 250). НУЭ = 2 478 625 (2.5%). 100 повторений
Относительно широкая лента. Все – предсказуемо.
OpenBLAS хорош только при малом числе потоков.
n = 10 000. Лента переменной ширины (nc = 250/25 для 10/90%). НУЭ = 453 625 (0.45%). 100 повторений
Да, решатели LAIPE для переменной ширины ленты вышли из аутсайдеров (особенно для одного потока), но до лидеров им далеко. В данном случае для этих решателей надо, чтобы разница в ширине ленты была куда больше.
n = 10 000. Узкая лента (nc = 25). НУЭ = 259 675 (0.26%). 100 повторений
Надо же! Какая же захватывающая получается гонка! Классический LAPACK, начав с самых последних мест, здесь уже в плотную подбирается к лидеру! Дышит ему в затылок! Это фантастика!
Очень заметно проявился уже отмеченный «баг» OpenBLAS. На 4 и 8 потоках его LLt неестественно тормозит. В одном потоке – он очень неплох.
n = 100 000. Очень узкая лента (nc = 5). НУЭ = 599 985 (0.006%). 100 повторений
Решатели LAIPE c переменной шириной ленты откинуты – они работали неприлично долго:
Та-дам! Йес! Мы с нетерпением ждали этого момента! Классический LAPACK вырывается вперед! Он – победитель! Он сделал MKL как стоячего! Это было невероятно захватывающее зрелище!
На больших и очень узких матрицах – распараллеливание и навороты, получается, неэффективны. (В протестированных реализациях). Классика побеждает!
Оле-е! Оле – оле – оле! LA-PACK – чемпи-он!
Заключение
Да уж. По скорости – все получилось очень запутано. Ясных закономерностей прослеживается немного.
Из того, что можно отметить:
Разница по времени между разными решателями – колоссальная, десятки и даже сотни раз, не смотря на то, что здесь мы вообще не рассматривали «тормозных» систем типа FreeMat. Только Фортран. Только хардкор. Только библиотеки – реальные претенденты на место «мотора» в наших реальных задачах.
Треугольное представление плотных матриц, похоже, почти не дает выигрыша во времени по отношению к полному плотному представлению. Во многих случаях – даже наоборот. Выигрыш – только по памяти. Может, конечно, на других n картина изменится.
Ленточные решатели, на ленточных матрицах – действительно сверх эффективны.
Ленточные решатели как правило даже на плотных матрицах мало уступают решателям для плотных матриц. Вывод: в принципе можно все матрицы мыслить как ленточные. Если они такие по факту – отлично, если нет – мы почти ни в чем не проигрываем.
На узких ленточных матрицах распараллеливание – не эффективно. (По крайней мере у рассмотренных вариантов). А в большинстве случаев – даже вредно. Иногда – очень вредно: решение на многих потоках замедлялось в десятки раз по сравнению с одноядерным вариантом. При этом на широких лентах – обычно все ровно наоборот.
На больших и очень узких ленточных матрицах классический LAPACK побеждает даже навороченный MKL. Если на плотных матрицах или матрицах с очень широкой лентой классический LAPACK был в аутсайдерах, то на больших матрицах с узенькой лентой он уверенно выходит вперед. Наблюдать за его рывком было чрезвычайно увлекательно – как будто Тур де Франс посмотрел.
Есть класс задач с очень сильно непостоянной шириной ленты, где учитывающие это решатели LAIPE2 работают хорошо, приближаясь к MKL. Но в большинстве случаев LAIPE2 – в аутсайдерах.
MKL – везде в лидерах. Самая молниеносная процедура – dpbsv (SB, LLt). Но поскольку свободность MKL мягко говоря довольно сомнительна, надо ориентироваться на аналоги. А они, все трое: OpenBLAS, классический LAPACK и LAIPE2, демонстрируют удивительное непостоянство: в каких-то вариантах проявляют рекордную тормознутость (в десятки, если не сотни раз медленнее, чем MKL), а в каких-то – вплотную приближаются к MKL. Причем так ведут себя все трое. :)
OpenBLAS – очень неровный. Где-то он отлично работает, где-то – явно недооптимизирован. Где-то странные тормоза при многопоточности, указываешь один поток – тормоза исчезают. Open Source, такой Open Source…
В рамках этой заметки был сделан единый интерфейс для целой плеяды решателей, работающих с различным форматом упакованных матриц. Теперь можно легко подключать их к реальным задачам и там экспериментальным путем подбирать наиболее подходящие. Как я уже писал выше, меня остро интересует не столько скорость (хотя при отличиях в скорости в десятки раз – это тоже критично), сколько устойчивость. Много копий с ней было сломано в практике: с одним решателем программа не работала хоть ты тресни, с другим, казалось бы по типу таким же, – ОК. До сих пор не было возможности протестировать этот феномен более вдумчиво.
Дальше есть планы изучить ScaLAPACK и несколько решателей для произвольно разряженных матриц. Их можно будет также замаскировать в рассмотренный интерфейс и подключать по необходимости. Как дойдут руки – сделаю тест и обзор и для них.
