среда, 26 апреля 2017 г.

GMSH 3.0

Вероятно, лучший в мире свободный мешер GMSH в апреле обновился до третьей версии. Почти все изменения связаны с появлением во встроенном сктриптовом языке, на котором можно описывать твердотельную геометрию системы, богатых возможностей технологии CSG (конструктивной блочной геометрии) на базе ядра OpenCASCADE:

  • появились команды для быстрого задания «готовых» примитивов: Sphere, Block, Torus, Cylinder, Cone, Wedge и другие;

  • появились полнофункциональные булевы операции над геометрическими объектами;

  • команды для «дополнительной обработки» типа скругления (Fillet) и т.п.;

  • возможность импорта геометрии и последующей работы с ней как с «родными» объектами.

В целом геометрический модуль внезапно стал очень мощным. В ряду «текстовых» CAD’ов для подготовки 3D-геометрии данный модуль GMSH занимает теперь достойное место (до третьей версии его вообще не приходилось серьезно рассматривать). Пожалуй, он вообще не уступает теперь знаменитому OpenSCAD, а во многом даже его превосходит, если учесть, что OpenSCAD работает только с сетками, а тут мы имеем «полноценное» описание на основе граничных поверхностей (BREP) на технологии OpenCASCADE.

В продолжение предыдущей заметки «Механические расчеты конструкций в свободном ПО. Часть 1» поупражняемся немного в языке GMSH. Тот же звездчатый многогранник «Соединение пяти тетраэдров», объединенный с шаром, на языке .geo GMSH можно построить так:

SetFactory("OpenCASCADE");

// Первый тетраэдр

Point(1) = { 1, 1, 1};
Point(2) = { 1,-1,-1};
Point(3) = {-1, 1,-1};
Point(4) = {-1,-1, 1};

Line(1) = {2,1};
Line(2) = {3,1};
Line(3) = {4,1};
Line(4) = {3,2};
Line(5) = {4,2};
Line(6) = {4,3};

Line Loop(1) = {1,5,3};
Line Loop(2) = {6,5,4};
Line Loop(3) = {2,3,6};
Line Loop(4) = {1,2,4};

Plane Surface(1) = {1};
Plane Surface(2) = {2};
Plane Surface(3) = {3};
Plane Surface(4) = {4};

Surface Loop(1) = {1,2,3,4};

Volume(1) = {1};

Dilate {{0,0,0}, 100/Sqrt(8)} {Volume{1};} // Масштабируем (длина ребра 100)
Rotate {{0,1,0}, {0,0,0}, -Atan(Sqrt(5)/2-0.5)} {Volume{1};} // Разворачиваем

// Остальные 4 тетраэдра

For i In {2:5}
    Rotate {{0,0,1}, {0,0,0}, (i-1)*2*Pi/5} {Duplicata{Volume{1};}}
EndFor

// Шар

Sphere(6) = {0,0,0, 64*Sqrt(6)/4};

// Объединение всех объектов

BooleanUnion(7) = {Volume{6}; Delete;} {Volume{1:5}; Delete;}; 

Код получился немного громоздким в силу того, что среди «готовых» геометрических примитивов нет тетраэдра: его пришлось строить по традиционной для GMSH идеологии «Точки -> Линии -> Замыкания линий -> Поверхности -> Замыкания поверхностей -> Тела». Если бы мы строили не тетраэдр, а, например, параллелепипед, призму или конус, все это заняло бы лишь одну строчку. Зато есть готовая команда для создания шара и булево объединение.

Язык, конечно, несколько «туповат». Лично мне больше всего не нравится именование примитивов по номерам, ну и общая бедность. С другой стороны, не надо относиться к нему как к языку программирования общего назначения. Для своих функций – описания геометрии не очень сложных объектов (а также последующей работы с ними) этот скриптовый язык вполне годный. Он на самом деле очень простой. На десятке маленьких учебных примеров изучается за пару часов.

В режиме отображения только точек и линий с их номерами получается следующая картина:


Забавно. Наш объект содержит 90 точек, 135 линий, 61 поверхность и одно тело.

После разбиения на конечные элементы второго порядка:




Здесь выполнено сечение объекта инструментом Clipping, чтобы посмотреть на элементы внутри.

В новой версии также улучшена и ускорена оптимизация разбиения.

Собственно, об основных возможностях GMSG, связанных с мешенгом, я планировал написать в продолжениях заметки «Механические расчеты конструкций в свободном ПО». Надеюсь, что когда-нибудь дойдут руки. Здесь попытался отметить именно нововведения в версии 3. В целом GMSH растет. Из «просто мешера» он уверено превращается в полноценную среду для различных конечноэлементных исследований. В сочетании с решателем GetDP из него теперь вообще можно не вылезать.

среда, 22 марта 2017 г.

Мой Git-конспект

Введение

В работе с Git я пользуюсь прекрасной оболочкой GitExtensions, а также расширениями к текстовым редакторам и IDE.

Поэтому помнить команды Git обычно нет никакой нужды – достаточно нажимать на кнопочки и горячие клавиши.

Есть лишь несколько команд и возможностей Git, выходящих за рамки GitExtensions, которые я постоянно забываю и каждый ищу заново и интернете. Чтобы сэкономить себе время, буду собирать их здесь.

Создание bare-репозитория

Bare-репозиторий – это репозиторий, в который не предполагается делать коммиты, а только пуши. Он предназначен не для работы, а для хранения изменений, сделанных в других (локальных) репозиториях. Обычно bare-репозиторий расположен на сервере, а локальные – на клиентах.

Для своих небольших проектов я создаю bare-репозитории в папке на DropBox. А локальные репозитории для работы – в специальном месте на каждой машине (иногда даже на RAM-диске, чтобы все просто летало).

Создаем bare-репозиторий:

  1. Создаем папку для него, например:
C:\Users\I\Dropbox\Rep\MyProject.git

Имя папки оканчивается на .git в силу соглашения об именовании bare-репозиториев. Из его имени сразу видно, что это не рабочий, а bare-репозиторий.

  1. Идем в консоль Git:
cd C:/Users/I/Dropbox/Rep/MyProject.git
git init --bare
exit

Обращаю внимание на прямые слеши вместо обратных в пути.

Все, bare-репозиторий – создан.

  1. Задаем в настройках GitExtensions «Внешний репозиторий» путь к нашему bare-репозиторию (имя указываем origin) и пушим туда наш проект.

Удаление из индекса ранее проиндексированных файлов

Если какое-то время мы работали и сохраняли в Git историю изменений некоторых файлов, а затем решили записать их в .gitignore, то они игнорироваться не станут, т.к. являются проиндексированными. Их надо убрать из индекса:

git rm --cached `git ls-files -i --exclude-from=.gitignore`

Или так, (это лучше работает):

git ls-files -i -z --exclude-from=.gitignore | xargs -0 git rm --cached

суббота, 11 марта 2017 г.

GFortran: конвертируем библиотеки .a в DLL

Бывает, что у нас имеются библиотеки .a только для статической линковки, а мы бы хотели вместо них иметь DLL и подключать ее динамически.

Примером может быть библиотека линейной алгебры LAIPE2 с стайта www.equation.com.

Оказывается, из .a очень просто сделать DLL автоматически.

Нам потребуется файл-заглушка Dllmain.f90 с пустой процедурой:

subroutine my_dll_main()
endsubroutine 

Откомпилируем его в DLL с подключением соответствующих библиотек (в данном примере – liblaipe2.a и нужная ей для работы libneuloop4.a):

gfortran dllmain.f90  -shared -o laipe2.dll -Wl,-export-all-symbols 
-Wl,-enable-auto-import -Wl,-no-undefined -Wl,--enable-runtime-pseudo-reloc 
-Wl,-whole-archive liblaipe2.a libneuloop4.a -Wl,-no-whole-archive

Тут довольно много ключей для линковщика – я их взял из примера, найденного в интернете, в смысл особо не погружался. Главное, что это работает.

Протестировать полученную DLL можно дополнив пример из прошлого поста:

! Вставка в файл LapackFromDll.F90:
...
    !> Интерфейсы функций LAIPE2
    abstract interface
        subroutine laipe_use_(N)
            integer  N
        endsubroutine
        subroutine laipe_done_()
        endsubroutine
        subroutine laipe_DAG_8_(A_io, N_i, X_io, NoGood_o)
            integer           NoGood_o, N_i
            double precision  A_io(*), X_io(*) 
        endsubroutine
    endinterface
    !> Процедуры-указатели на функции LAIPE2
    procedure(laipe_use_)  , pointer  :: laipe_use
    procedure(laipe_done_) , pointer  :: laipe_done
    procedure(laipe_DAG_8_), pointer  :: laipe_DAG_8
...
! Вставка в файл LinTest.F90:
...
    call LoadDll(LinAlg_dll, "laipe2.dll")
    call DllProcAssoc(LinAlg_dll, laipe_use,   "laipe$use_")
    call DllProcAssoc(LinAlg_dll, laipe_done,  "laipe$done_")
    call DllProcAssoc(LinAlg_dll, laipe_DAG_8, "laipe$solution_dag_8_")
    call laipe_use(8)
    A = A0 ; X = B 
    call laipe_DAG_8(A, n, X, info)
    print "(' *** Решение с помощью laipe2.dll (info = ',i1,'):')", info
    print FRM, X
    call laipe_done()
    call FreeDll(LinAlg_dll) 
...

пятница, 10 марта 2017 г.

GFortran: вызываем функции LAPACK из DLL динамически

Как известно, существует два способа связывания DLL-библиотек с остальной частью программы:

  1. Статическое, или неявное связывание. Все заботы берет на себя компилятор (точнее, линковщик), вы просто вызываете функцию, содержащуюся в DLL. Для этого ему нужно иметь соответствующую (сопоставленную) нашей DLL библиотеку *.a (или *.lib), которая содержит необходимую для линковщика информацию. Все связанные так DLL подключаются сразу после запуска программа.

  2. Динамическое или явное связывание. Может быть выполнено в любое время работы программы и выполняется программистом «руками». Библиотеки *.a (или *.lib) – не требуются.

Недостаток первого метода для нас заключается в том, что все DLL должны иметься на момент запуска программы. Если какой-то DLL нет, то программа не запустится. А нам иногда очень хочется подключить ту или иную DLL в зависимости от настроек. Кроме того, первый способ не позволяет указать путь к нужным DLL, из-за чего может получиться, что подключатся какие-то не те версии.

Например, как мы уже отмечали, библиотека линейной алгебры LAPACK имеется как минимум в трех реализациях: классической, ObenBLAS и MKL (вообще-то, их куда больше). Для экспериментов нам хотелось бы иметь возможность переключатся между ними «на лету». Пусть эти библиотеки имеются в отдельных DLL. (На самом деле, в каждом случае нужна не одна, а несколько связанных библиотек). Напишем модуль, реализующий их динамическое подключение (на примере трех функций LAPACK: dgesv, dposv и dsysv). На Фортране он получается чуть-чуть мудреным, но в принципе, все не так страшно:

!> ******   Модуль с интерфейсом подключения процедур LAPACK из DLL   ******
module LapackFromDll
    use, intrinsic :: iso_c_binding, only:  c_f_procpointer, &
         c_funptr, c_intptr_t, c_null_char, c_char, c_associated, c_bool

    implicit none

    !> Интерфейсы функции WinAPI, работающих с DLL
    interface
        ! Интерфейс загрузки DLL (из Kernel32.lib)
        function LoadLibrary(lpFileName) bind(c, name='LoadLibraryA')
            use, intrinsic :: iso_c_binding, only: c_intptr_t, c_char
            implicit none 
            character(kind=c_char)          :: lpFileName(*) 
            !GCC$ attributes stdcall        :: LoadLibrary 
            integer(c_intptr_t)             :: LoadLibrary 
        endfunction 
        ! Интерфейс получения адреса функции DLL (из Kernel32.dll)
        function GetProcAddress(hModule, lpProcName) bind(c, name='GetProcAddress')
            use, intrinsic :: iso_c_binding, only: c_funptr, c_intptr_t, c_char
            implicit none
            !GCC$ attributes stdcall        :: GetProcAddress
            type(c_funptr)                  :: GetProcAddress
            integer(c_intptr_t), value      :: hModule
            character(kind=c_char)          :: lpProcName(*)
        endfunction 
        ! Интерфейс освобождения DLL (из Kernel32.dll)
        function FreeLibrary(hModule) bind(c, name='FreeLibrary')
            use, intrinsic :: iso_c_binding, only: c_funptr, c_intptr_t, c_bool
            implicit none
            !GCC$ attributes stdcall        :: FreeLibrary
            logical(c_bool)                 :: FreeLibrary
            integer(c_intptr_t), value      :: hModule
        endfunction  
    endinterface 

    !> Интерфейсы функций LAPACK
    abstract interface
        subroutine dgesv_(N, NRHS, A, LDA, IPIV, B, LDB, INFO) 
            integer            INFO, LDA, LDB, N, NRHS
            integer            IPIV(*)
            double precision   A(lda,*), B(ldb,*)
        endsubroutine
        subroutine dposv_(UPLO, N, NRHS, A, LDA, B, LDB, INFO) 
            character          UPLO
            integer            INFO, LDA, LDB, N, NRHS
            double precision   A(LDA,*), B(LDB,*)
        endsubroutine   
        subroutine dsysv_(UPLO, N, NRHS, A, LDA, IPIV, B, LDB, WORK, LWORK, INFO)
            character          UPLO
            integer            INFO, LDA, LDB, LWORK, N, NRHS
            integer            IPIV(*)
            double precision   A(LDA,*), B(LDB,*), WORK(*)
        endsubroutine  
    endinterface

    !> Процедуры-указатели на функции LAPACK
    procedure(dgesv_), pointer  :: dgesv        
    procedure(dposv_), pointer  :: dposv
    procedure(dsysv_), pointer  :: dsysv        

 contains  ! ============== Обертки функций для работы с DLL =================

    !> Загружает dll и связывает с ним Dll_handle
    subroutine LoadDll(Dll_handle, Dll_name)
        integer(c_intptr_t)     :: Dll_handle
        character(*, c_char)    :: Dll_name
        !------------------------------------------------------
        Dll_handle = LoadLibrary(Dll_name // c_null_char)
        if (Dll_handle == 0) stop 'Не удалось загрузить ' // Dll_name
    endsubroutine

    !> Ассоциирует процедуру-указатель Proc с функцией из Dll
    subroutine DllProcAssoc(Dll_handl, Proc, Proc_name)
        integer(c_intptr_t)     :: Dll_handl
        character(*, c_char)    :: Proc_name
        procedure(), pointer    :: Proc 
        type(c_funptr)          :: proc_address
        !------------------------------------------------------
        proc_address = GetProcAddress(Dll_handl, Proc_name // c_null_char)
        if (.not. c_associated(proc_address)) stop 'Не могу найти адрес процедуры '&
                                                    // Proc_name // " в DLL" 
        call c_f_procpointer(proc_address, Proc)
    endsubroutine

    !> Освобождает связанные с DLL ресурсы
    subroutine FreeDll(Dll_handle)
        integer(c_intptr_t)       :: Dll_handle
        logical(c_bool)           :: FrDll
        !------------------------------------------------------
        FrDll = FreeLibrary( Dll_handle )
        if (.not. FrDll) stop 'Не удалось выгрузить DLL'
    endsubroutine    

endmodule

Интерфейсы функций LAPACK можно копировать прямо из исходников с сайта www.netlib.org/lapack/.

Проверяем. Вызовем функцию dgesv «на лету» последовательно из трех реализаций LAPACK – классической, ObenBLAS и MKL:

program LinTest
    use, intrinsic :: iso_c_binding, only: c_intptr_t
    use  LapackFromDll
    implicit none
    integer(c_intptr_t)      :: LAPACK_dll ! Хендл LAPACK_dll
    integer, parameter       :: n = 10     ! Размерность матриц
    real(8), dimension(n,n)  :: A, A0      ! Матрица коэффициентов
    real(8), dimension(n,1)  :: B, X       ! Векторы правых частей / рез.
    integer                  :: i,j        ! Вспомогательный индекс
    integer, dimension(n)    :: ipiv       ! Вспомогательные для LAPACK
    integer                  :: info       ! Инфо об успешности решения
    character(100)           :: FRM        ! Формат для вывода матриц
    !---------------------------------------------------------------------
    print *, "Проверка динамической загрузки процедур LAPACK из DLL"
    forall(i=1:n, j=1:n) A0(i,j) = 1-(dble(i-j)/n)**2  ! Формируем матрицу
    forall(i=1:n, j=1:n, i==j) A0(i,j) = 10.
    forall(i = 1:n) B(i,1) = 100*i        ! Заполняем вектор B
    write(FRM, "(a,i10,a)") "(",n,"f8.2)" ! Формат для вывода матриц
    print *, "Матрица коэффициентов A:"
    print FRM, transpose(A0)
    print *, "Вектор правых частей B:"
    print FRM, B

    call LoadDll(LAPACK_dll, "liblapack.dll")
    call DllProcAssoc(LAPACK_dll, dgesv, "dgesv_")
    A = A0 ; X = B 
    call dgesv(n, 1, A, n, ipiv, X, n, info)  
    print "(' *** Решение с помощью liblapack.dll (info=',i1,'):')", info
    print FRM, X
    call FreeDll(LAPACK_dll)

    call LoadDll(LAPACK_dll, "libopenblas.dll")
    call DllProcAssoc(LAPACK_dll, dgesv, "dgesv_")
    A = A0 ; X = B 
    call dgesv(n, 1, A, n, ipiv, X, n, info)  
    print "(' *** Решение с помощью libopenblas.dll (info=',i1,'):')", info
    print FRM, X
    call FreeDll(LAPACK_dll)

    call LoadDll(LAPACK_dll, "mkl_rt.dll")
    call DllProcAssoc(LAPACK_dll, dgesv, "dgesv_")
    A = A0 ; X = B 
    call dgesv(n, 1, A, n, ipiv, X, n, info)  
    print "(' *** Решение с помощью mkl_rt.dll (info=',i1,'):')", info
    print FRM, X
    call FreeDll(LAPACK_dll)    
endprogram 

Компилируем:

gfortran -c LapackFromDll.F90
gfortran -c LinTest.F90
gfortran *.o -o LinTest.exe

Получаем:

 Проверка динамической загрузки процедур LAPACK из DLL
 Матрица коэффициентов A:
  10.00    0.99    0.96    0.91    0.84    0.75    0.64    0.51    0.36    0.19
   0.99   10.00    0.99    0.96    0.91    0.84    0.75    0.64    0.51    0.36
   0.96    0.99   10.00    0.99    0.96    0.91    0.84    0.75    0.64    0.51
   0.91    0.96    0.99   10.00    0.99    0.96    0.91    0.84    0.75    0.64
   0.84    0.91    0.96    0.99   10.00    0.99    0.96    0.91    0.84    0.75
   0.75    0.84    0.91    0.96    0.99   10.00    0.99    0.96    0.91    0.84
   0.64    0.75    0.84    0.91    0.96    0.99   10.00    0.99    0.96    0.91
   0.51    0.64    0.75    0.84    0.91    0.96    0.99   10.00    0.99    0.96
   0.36    0.51    0.64    0.75    0.84    0.91    0.96    0.99   10.00    0.99
   0.19    0.36    0.51    0.64    0.75    0.84    0.91    0.96    0.99   10.00
 Вектор правых частей B:
 100.00  200.00  300.00  400.00  500.00  600.00  700.00  800.00  900.00 1000.00
 *** Решение с помощью liblapack.dll (info=0):
  -6.20    0.36    7.63   15.60   24.29   33.68   43.77   54.58   66.09   78.30
 *** Решение с помощью libopenblas.dll (info=0):
  -6.20    0.36    7.63   15.60   24.29   33.68   43.77   54.58   66.09   78.30
 *** Решение с помощью mkl_rt.dll (info=0):
  -6.20    0.36    7.63   15.60   24.29   33.68   43.77   54.58   66.09   78.30

понедельник, 6 марта 2017 г.

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

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


Содержание

Введение

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

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

  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 и несколько решателей для произвольно разряженных матриц. Их можно будет также замаскировать в рассмотренный интерфейс и подключать по необходимости. Как дойдут руки – сделаю тест и обзор и для них.