Тест решателей СЛАУ с разряженными симметрическими матрицами в плотном и ленточном представлении

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

Содержание

• Введение
• Форматы представления симметрических матриц
  • G – общий формат
  • ST – упакованный треугольник
  • STL – упакованный треугольник для LAIPE2
  • GB – упакованная ленточная матрица общего вида
  • SB – упакованная симметрическая ленточная матрица
  • SBL – упакованная симметрическая ленточная матрица для LAIPE2
  • SBLV – упакованная симметрическая ленточная матрица с переменной шириной ленты для LAIPE2
  • SBI – упакованная симметрическая ленточная матрица для IMSL
• Интерфейс решателей для симметрических матриц
  • Описание
  • Код модуля
• Оценка производительности решателей
  • Тестовые задачи
  • Решение
  • Сравнительные результаты
    • n = 1000. Не лента (nc = 999). НУЭ = 500 500 (50%). 100 повторений
    • n = 1000. Лента (nc = 250). НУЭ = 219 625 (22%). 100 повторений
    • n = 1000. Лента переменной ширины (nc = 250/25 для 10/90%). НУЭ = 28 451 (2.85%). 100 повторений
    • n = 1000. Узкая лента (nc = 25). НУЭ = 25 675 (2.57%). 100 повторений
    • n = 10 000. Лента (nc = 250). НУЭ = 2 478 625 (2.5%). 100 повторений
    • n = 10 000. Лента переменной ширины (nc = 250/25 для 10/90%). НУЭ = 453 625 (0.45%). 100 повторений
    • n = 10 000. Узкая лента (nc = 25). НУЭ = 259 675 (0.26%). 100 повторений
    • n = 100 000. Очень узкая лента (nc = 5). НУЭ = 599 985 (0.006%). 100 повторений
• Заключение

Введение

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

Чтобы эффективно решать СЛАУ с сильно разреженными матрицами можно пойти двумя путями:

1. Использовать специальные библиотеки для разреженных матриц – такие как PARDISO, MUMPS, SPOOLES, TAUCS, SuperLU, UMFPACK, PSPASES, WSMP, PaStiX и другие.

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

В обоих случаях в компьютерной памяти хранится только «содержательная» часть матрицы. Для варианта 1 – ненулевые компоненты и их «адреса», для варианта 2 – лента.

Получить ленточную матрицу можно тоже двумя путями:

1. В алгоритме решения базовой задачи позаботиться, чтобы матрица получалась ленточной «автоматом». Для МКЭ это означает, что нужно правильно организовать нумерацию узлов конечных элементов в системе. Это может быть чрезвычайно эффективно для всяких «вытянутых» и «длинных» систем. Скажем, мы на работе обсчитываем контактную подвеску железных дорог. Она – очень-преочень длинная, для нее (если позаботиться о правильной нумерации узлов), ленточная матрица получается естественным образом.

2. Второй путь – сначала сделать разреженную матрицу «как придется», а затем преобразовать ее в ленточную посредством перенумерации строк и столбцов. Для этого существуют эффективные алгоритмы, такие как GPS, GGPG, RCM, Sloan, JCL и другие. В случае с МКЭ можно автоматически преренумеровывать не столбцы и строки, а номера узлов конечных элементов т.е. номера блоков в блочном представлении матрицы. В одной из задач мы использовали для этого алгоритм JCL (сайт). Результаты превзошли всякие ожидания: программа из типа «Запустил – пойду пообедаю» стала интерактивной – расчет стал выполняться секунды.

Проиллюстрируем перенумерацию узлов алгоритмом JCL для простейшей системы и ее блочной матрицы. Слева – до, справа – после. Цифры в матрице – ненулевые элементы у блоков.

Различных подпрограмм для решения СЛАУ с ленточными матрицами – тоже множество. Для нас во-первых интересно оценить скорость их работы, а во-вторых – реализовать некоторый единый интерфейс для работы с ними, скрывающий различные форматы представления матриц. Это позволит пробовать разные решатели в наших «околонаучных» программах и менять их «на лету».

Почему это важно? Дело в том, что разные решатели отличаются не только скоростью. Они также различаются параметрами устойчивости алгоритмов и имеют разные «пределы терпимости» к плохо обусловленным матрицам, которые иногда получаются в жизни. В некоторых задачах мы сильно намучились с устойчивостью решений. Хочется опробовать разные решатели в реальных программах, чтобы выбрать самые «некапризные» в плане устойчивости и приемлемые по скорости. Сделать это без единого интерфейса – очень проблемно.

Как оценить устойчивость в синтетическом тесте – не очевидно. Устойчивость – скорее характеристика решения реальной задачи, поэтому «тест на устойчивость» в рамках данной заметки мы делать не будем. Ограничимся разбирательством с форматами ленточных матриц, формированием интерфейса решателей и «бытовым» тестом их скоростных качеств.

Все «интерфейсы» и тесты будут сделаны на Фортране с компилятором GFortran. В прошлый раз мы убедились, что в разных языках решатели СЛАУ вызываются одни и те же (в виде скомпилированных библиотек). Значит – переписывать одно и то же на разных языках нет смысла.

Будут рассмотрены решатели и соответствующие форматы матриц из следующих библиотек:

1. LAPACK (+BLAS) в трех реализациях: классической (версия 3.7.0), OpenBLAS (0.2.19) и Intel MKL (11.3 – он поставился у меня из дистрибутива Python Anaconda).

2. 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(:)	Вход / выход	Вектор правых частей или решений. На входе – размещенный в памяти, на выходе – освобожденный.

Особенности:

1. Между вызовами SlaeBeginAlloc и SlaeEndDealloc можно работать только с одними экземплярами матриц (в модуле имеются внутренние данные, имеющиеся в единственном числе).

2. Формат SBI поддерживается только на уровне выделения памяти, упаковки и адресации к элементам.

3. Какие библиотеки применяются с модулем (классический LAPACK, OpenBLAS, MKL, LAIPE2) – управляется директивами компилятора.

4. В модулей введен тип данных, определяющий вариант решателя СЛАУ и имеется вспомогательная процедура его установки по условному номеру SetSlaeSolution.

5. При работе с форматом ленты переменной ширины 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 задач. В качестве самой первой повторим задачу, описанную в прошлом посте.

№	Число уравнений n	Полуширина ленты nc	НУЭ*	Примечание
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 повторений

Что видно здесь интересного:

1. Колоссальная разница по времени между самым быстрым и самым медленным решателем (порядка 100 раз).

2. Треугольное представление матриц в основном не дает выигрыша во времени, скорее наоборот. Выигрыш – только по памяти.

3. В OpenBLAS есть хорошо оптимизированные решатели, а есть – откровенно плохо (LDLt).

4. MKL в этом смысле более ровный. Все в нем оптимизировано дальше некуда. Он – очень крут.

5. В целом даже если матрица не ленточная, ленточные решатели во многих случаях такие же быстрые как и решатели для обычного представления матрицы.

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 – чемпи-он!

Заключение

Да уж. По скорости – все получилось очень запутано. Ясных закономерностей прослеживается немного.

Из того, что можно отметить:

1. Разница по времени между разными решателями – колоссальная, десятки и даже сотни раз, не смотря на то, что здесь мы вообще не рассматривали «тормозных» систем типа FreeMat. Только Фортран. Только хардкор. Только библиотеки – реальные претенденты на место «мотора» в наших реальных задачах.

2. Треугольное представление плотных матриц, похоже, почти не дает выигрыша во времени по отношению к полному плотному представлению. Во многих случаях – даже наоборот. Выигрыш – только по памяти. Может, конечно, на других n картина изменится.

3. Ленточные решатели, на ленточных матрицах – действительно сверх эффективны.

4. Ленточные решатели как правило даже на плотных матрицах мало уступают решателям для плотных матриц. Вывод: в принципе можно все матрицы мыслить как ленточные. Если они такие по факту – отлично, если нет – мы почти ни в чем не проигрываем.

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

6. На больших и очень узких ленточных матрицах классический LAPACK побеждает даже навороченный MKL. Если на плотных матрицах или матрицах с очень широкой лентой классический LAPACK был в аутсайдерах, то на больших матрицах с узенькой лентой он уверенно выходит вперед. Наблюдать за его рывком было чрезвычайно увлекательно – как будто Тур де Франс посмотрел.

7. Есть класс задач с очень сильно непостоянной шириной ленты, где учитывающие это решатели LAIPE2 работают хорошо, приближаясь к MKL. Но в большинстве случаев LAIPE2 – в аутсайдерах.

8. MKL – везде в лидерах. Самая молниеносная процедура – dpbsv (SB, LLt). Но поскольку свободность MKL мягко говоря довольно сомнительна, надо ориентироваться на аналоги. А они, все трое: OpenBLAS, классический LAPACK и LAIPE2, демонстрируют удивительное непостоянство: в каких-то вариантах проявляют рекордную тормознутость (в десятки, если не сотни раз медленнее, чем MKL), а в каких-то – вплотную приближаются к MKL. Причем так ведут себя все трое. :)

9. OpenBLAS – очень неровный. Где-то он отлично работает, где-то – явно недооптимизирован. Где-то странные тормоза при многопоточности, указываешь один поток – тормоза исчезают. Open Source, такой Open Source…

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

Дальше есть планы изучить ScaLAPACK и несколько решателей для произвольно разряженных матриц. Их можно будет также замаскировать в рассмотренный интерфейс и подключать по необходимости. Как дойдут руки – сделаю тест и обзор и для них.