понедельник, 27 февраля 2017 г.

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

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


Введение

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

От скорости решения СЛАУ тем или иным средством может зависеть сама возможность решения исходной задачи, ну и, конечно, удобство реализации. Если для проведения одного вычислительного эксперимента приходится ждать много часов или даже дней – это крайне досадно. Оптимизированная реализация той же задачи может считаться считанные минуты.

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

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

  • эффективность расчета разными средствами может сильно зависеть от типа матрицы СЛАУ (и выбора соответствующего метода решения);

  • у разных средств может сильно различаться масштабируемость относительно размера матрицы;

  • «под капотом» различных средств могут лежать разные (или одни и те же) решатели СЛАУ, обычно это скомпилированные библиотеки. Иногда в рамках одного средства можно выбрать различные решатели, и каждый решатель тоже можно настроить по-разному (например, на работу в одно- или многопоточном режиме). Сейчас мы разбираться с этим не будем, используем все «из коробки»;

  • на разных машинах (например, с разным количеством ядер) эффективность различных средств может очень сильно разниться;

  • с решением СЛАУ в инженерных задачах часто связана основная доля потери времени работы программы, но далеко не единственная. Скорость работы остального кода на разных языках и в разной реализации тоже может различаться в разы или десятки, сотни, тысячи… :) раз.

Тем не менее, для первичной оценки такой тест провести интересно.

В качестве матрицы коэффициентов СЛАУ выберем симметрическую положительно определенную матрицу с вещественными коэффициентами. К таким матрицам сводится множество физических задач. При этом матрицу коэффициентов будем считать (и сделаем) плотной. Вообще говоря, в реальной практике значительно чаще получаются не плотные, а сильно разряженные матрицы. Разработка эффективных алгоритмов для таких матриц – отдельное направление современной науки. Хорошо известны решатели PARDISO, MUMPS, SPOOLES, TAUCS, SuperLU, UMFPACK, PSPASES, WSMP, PaStiX и другие. При работе с разряженными матрицами на любом языке обычно вызывается одна из этих программ. Их настройка «в реальной жизни» и тестирование – очень важная и интересная для нас тема, но сейчас пока оставим ее за скобками. Разряженные матрицы иногда можно свести к ленточным, что чрезвычайно эффективно при малой ширине ленты. Ленточные матрицы тоже пока не рассматриваем.

Те средства, где доступны специальные методы решения СЛАУ с симметрической положительно определенной матрицей (на основе \(LL^T\) или \(U^TU\)-разложения Холецкого), протестируем два раза: с решением обычным методом для матриц общего вида (на основе классического \(LU\)-разложения) и методом с \(LL^T\)-разложением.

Будем тестировать следующие программные средства:

  • язык Fortran с различными математическими библиотеками (реализациями LAPACK/BLAS и другими);

  • язык Julia «из коробки», в котором доступны как «стандартные средства» решения СЛАУ (скрывающие используемый метод), так и вызовы конкретных процедур библиотеки LAPACK;

  • язык Python с математическими библиотеками NumPy и SciPy из дистрибутива Anaconda;

  • MatLAB-подобные открытые системы компьютерной математики (СКМ), установленные «из коробки»: Octave, SсiLAB и FreeLAB;

  • некоторые «прочие» открытые СКМ. Тут пока остановимся на СКА Maxima.

Все тесты проводились на обычном современном персональном компьютере с 4-ядерным 8-поточным процессором Intel i7 с достаточным размером памяти. Конкретную марку, частоты и т.п., думаю, описывать нет нужды, т.к. нас интересует лишь относительное сравнение средств.

ОС: Windows 10, 64. Вест софт (где это было доступно), использовался в 64-разрядной реализации.

Тестовая задача

Решим \(m\) раз систему из \(n\) уравнений

\[ \mathbf{A}\cdot\mathbf{X}=\mathbf{B}, \]

где \(\mathbf{A}\) – матрица коэффициентов размером \(n \times n\), \(\mathbf{X}\) – вектор неизвестных и \(\mathbf{B}\) – вектор правых частей.

Пусть коэффициенты матрицы \(\mathbf{A}\) вычисляются следующим образом:

\[ a_{i,j} =
\begin{cases}
1- \left(\displaystyle\frac{i-j}{n}\right)^2, &\text{если } i \ne j,\\
10, &\text{если } i = j.
\end{cases}
\]

При этом получается симметрическая положительно определенная матрица.

Например, для \(n=5\) матрица \(\mathbf{A}\) будет такой:

  10.00   0.96   0.84   0.64   0.36
   0.96  10.00   0.96   0.84   0.64
   0.84   0.96  10.00   0.96   0.84
   0.64   0.84   0.96  10.00   0.96
   0.36   0.64   0.84   0.96  10.00

Пусть на каждом шаге решения \( k = 1 \dots m \) коэффициенты вектора правых частей равны номеру шага: \( b_{i} = k \).

Чтобы контролировать правильность решения каждым программным средством будем подсчитывать сумму коэффициентов вектора \(\mathbf{X}\) на каждом шаге решения и суммировать ее по шагам:

\[ S = \sum_{k=1}^{m} \sum_{i=1}^{n} x_i^{(k)}. \]

Все программы будем тестировать при размере матрицы \(n=1000\) и числе повторений \(m=100\) шагов. При этих параметрах должно получаться \( S \approx 22 439,195\). Во всех программах будем использовать векторы и матрицы из вещественных чисел двойной точности (занимающие в памяти 8 байт).

Договоримся, что разложение (факторизацию) матрицы \(\mathbf{A}\) будем делать на каждом шаге, т.е. каждый раз систему уравнений будем решать полностью.

Время решения будем измерять по циклу решения СЛАУ – без учета предварительной подготовки матрицы \(\mathbf{A}\).

Разные запуски реализаций решения на ПК дают несколько разное время т.к. ОС периодически отбирает ресурсы от нашей программы для своих нужд. Из нескольких запусков будем записывать минимальное время.

Решение различными средствами

Fortran

GFortran + BLAS/LAPACK (LU)

Для начала опробуем решение задачи на GNU GFortran (версия 5.3) с классическими библиотеками BLAS и LAPACK (версия 3.7.0). Установка описана в прошлом посте.

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

program LinTest
    implicit none
    integer, parameter      :: n = 1000 ! Число уравнений
    integer, parameter      :: m = 100  ! Число шагов
    real(8), dimension(n,n) :: A, A0    ! Матрица коэффициентов
    real(8), dimension(n)   :: B        ! Вектор правых частей / результатов
    real(8)                 :: s        ! Контрольная сумма
    integer                 :: i, j, k  ! Вспомог. индексы
    integer, dimension(n)   :: ipiv     ! Вспомог. матрица для процед. LAPACK
    integer                 :: info     ! Инфо об успешности решения
    integer                 :: t_0, t_1, t_rate, t_max ! Для измер. времени
    !------------------------------------------------------------------------
    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.
    s = 0.
    call system_clock(t_0, t_rate, t_max) 
    do k = 1, m
        A = A0
        B = k   ! Формируем вектор правых частей
        call dgesv(n, 1, A, n, ipiv, B, n, info) ! Решаем СЛАУ 
        if (info/=0) stop ("Ошибка")
        s = s + sum(B)
    enddo
    call system_clock(t_1, t_rate, t_max)
    print *, "Время:", real(t_1-t_0)/t_rate
    print *, "Сумма:", s
endprogram

Здесь использована процедура LAPACK dgesv, реализующая решение СЛАУ с \(LU\)-разложением и подходящая для матриц общего вида.

Матрицу A приходится восстанавливать из эталонной A0 на каждом шаге цикла т.к. dgesv ее изменяет (как и вектор B, куда записывается решение).

Компиляция:

gfortran LinTest.f90 -static -llapack -lrefblas -O2 -o LinTest.exe

Время решения задачи в таком варианте составило 14,6 с. При этом задействовался только один процессор. Библиотеки LAPACK/BLAS в классическом варианте – непараллельные. Видно, что в данном варианте они не оптимизированы – время расчета вызывает разочарование.

Размер exe-файла получился 821 КБ.

Компиляция с более агрессивными ключами оптимизации ускорения почти не дала.

Если вынести факторизацию матрицы

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

Для решения СЛАУ с матрицей общего вида в два приема в LAPACK имеются функции dgetrf и dgetrs. Первая реализует \(LU\)-разложение, вторая на его основе решает СЛАУ. Фрагмент кода в зоне главного цикла станет таким:

call dgetrf(n, n, A0, n, ipiv, info)  ! LU-разложение
do k = 1, m
    B = k
    call dgetrs("N", n, 1, A0, n, ipiv, B, n, info)  ! Решение СЛАУ
    if (info/=0) stop ("Ошибка") 
    s = s + sum(B)
enddo

Время решения задачи стало всего лишь 0,21 с, т.е. в 70 раз меньше (!). Оно и понятно. Трудоемкая часть в решении СЛАУ – факторизация.

Результат решения, разумеется, не изменился т.к. матрица \(\mathbf{A}\) на разных шагах решения по условию задачи у нас одинаковая.

GFortran + BLAS/LAPACK (LLt)

Для решения СЛАУ с симметрическими положительно определенными матрицами в LAPACK имеются процедуры на основе \(LL^T\)-разложения Холецкого. В исходной программе достаточно заменить одну строчку:

call dposv("L", n, 1, A, n, B, n, info) ! Решаем СЛАУ 

Время решения составило 8,6 с.

Размер exe-файла – примерно тот же.

Если вынести факторизацию матрицы

Фрагмент кода в зоне главного цикла:

call dpotrf("L", n, A0, n, info)  ! LU-разложение 
do k = 1, m 
    B = k
    call dpotrs("L", n, 1, A0, n, B, n, info)  ! Решение СЛАУ
    if (info/=0) stop ("Ошибка") 
    s = s + sum(B)
enddo

Время решения стало 0,16 с – в 55 раз меньше.

GFortran + OpenBLAS (LU)

OpenBLAS (www.openblas.net/) – свободная библиотека с функциями BLAS+LAPACK, оптимизированная для современных процессоров.

Устанавливается в GFortran OpenBLAS очень просто. Достаточно скачать версию со скомпилированными бинарниками и положить библиотеку libopenblas.a туда, где GCC хранит библиотеки .a. Текущая версия 0.2.19.

Код программы идентичен предыдущим примерам.

Компиляция:

gfortran LinTest.f90 -static -lopenblas -O2 -o LinTest.exe

Время решения составило 0,79 с. При этом задействовались все процессоры.

Слинкованный статически exe-файл получился очень большим: 26 (!) мегабайт. Можно ли его уменьшить я не разобрался. В дистрибутиве OpenBLAS прилагается также dll-библиотека. Но динамическая линковка с ним у меня не заработала.

Вынос факторизации за тело цикла уменьшил время до 0,03 с (!).

GFortran + OpenBLAS (LLt)

Время решения составило 0,51 с.

Это – шикарно!

Если бы еще размер exe-файла получался поменьше…

Вынос факторизации за тело цикла уменьшил время до 0,07 с.

GFortran + LAIPE2 (LU)

Библиотека Laipe2 (www.equation.com) и способ ее установки были упомянуты в предыдущем посте.

Это оригинальная реализация параллельных решателей СЛАУ с отличной от LAPACK идеологией. Для матриц общего вида вызов совмещенного решателя осуществляется так:

call laipe$Solution_DAG_8(A, n, B, info) ! Вызов процедуры LAIPE2

(Доступны также процедуры с разложениями и последующим решением).

Предварительно нужно инициализировать используемое число ядер (потоков) процессора (в моей случае – максимум 8). Без этого работать не будет:

call laipe$use(8)

Компиляция:

gfortran LinTest.f90 -static -fdollar-ok -llaipe2 -lneuloop4 -O2 -o LinTest.exe

Время решения – 7,4 с. Не густо. При этом все ядра загружаются максимально.

Размер exe-файла – 910 КБ. Это – приятно.

GFortran + LAIPE2 (LLt)

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

Приведу полный код – он отличается от прошлых реализаций:

program LinTest
    implicit none
    integer, parameter      :: n = 1000 ! Число уравнений
    integer, parameter      :: m = 100  ! Число шагов
    real(8), dimension(((n+1)*n)/2) :: A, A0 ! Профиль матрицы коэффициентов
    real(8), dimension(n)   :: B        ! Вектор правых частей / результатов
    integer, dimension(n)   :: Label    ! Массив меток
    real(8)                 :: s        ! Контрольная сумма      
    integer                 :: i, j, k  ! Вспомог. индексы
    integer                 :: info     ! Инфо об успешности решения
    integer                 :: t_0, t_1, t_rate, t_max ! Для измер. времени
    !------------------------------------------------------------------------
    Label(1) = 1      ! Формируем массив меток
    do j = 2, n
        Label(j) = Label(j-1) + n-j+1
    enddo
    do j = 1, n       ! Формируем профиль матрицы коэффициентов
        do i = j, n
            A0(Label(j)+i-1) = 1-(dble(i-j)/n)**2  ! j-столбец, i-строка
            if (i==j) A0(Label(j)+i-1) = 10.
        enddo
    enddo
    s = 0.
    call laipe$use(8) ! Число используемых потоков
    call system_clock(t_0, t_rate, t_max)
    do k = 1, m       ! Решение задачи
        A = A0
        B = k
        call laipe$Solution_DSP_8(A, n, Label, B, info)  ! Вызов LAIPE2
        if (info==1) stop ("Ошибка") 
        s = s + sum(B)
    enddo
    call system_clock(t_1, t_rate, t_max)
    print *, "Время:", real(t_1-t_0)/t_rate
    print *, "Сумма:", s
endprogram

Время решения в данной реализации составило 4,9 с. Ядра загружаются максимально.

Размер exe-файла – 910 КБ.

Intel Fortran + MKL (LU)

В интернете оказался доступен триал коммерческого компилятора Intel Fortran с математической библиотекой MKL, содержащей все те же BLAS и LAPACK, но хорошо оптимизированные для современных процессоров. Он стоит немалых денег. Хотя настоящий блог посвящен свободному ПО, я не мог не попробовать. Очень интересно сравнить.

Код программы остается все тем же, что и для классических BLAS/LAPACK.

Компиляция:

ifort /Qmkl /static /O2 LinTest.f90

Время решения составило 0,58 с. Конечно, великолепно.

Размер exe-файла – 6,3 МБ. Тоже немало, но не 26 все же. :)

Intel Fortran + MKL (LLt)

Время решения составило всего 0,33 с. Круто, что скажешь. Забегая вперед, это абсолютный рекорд из всех протестированных решений. Программисты Intel едят хлеб не зря. Intel Fortran, определенно, своих денег стоит.

Julia

У меня установлена Julia 0.5.0. У нее «под капотом» – реализация LAPACK OpenBlas, очевидно, аналогичная той, которую мы уже «щупали» под Fortran (может быть версия постарее).

Julia, стандартное средство (LU)

Код решения задачи на Джулии – лаконичен и понятен без комментариев:

n = 1000; m = 100
s = 0.
A = [1-((i-j)/n)^2 for i in 1:n, j in 1:n] + 9*eye(n,n)
  @time for k=1:m
    B = fill(Float64(k), n)
    X = A \ B
    s += sum(X)
  end
println(s)

По сравнению с кодом на Fortran, что называется, «почувствуйте разницу».

Время решения в таком варианте составило 1,2 с. Задействовались все ядра.

Макрос @time выводит как время выполнения цикла, а также объем выделенной памяти. В данном случае 765 МБ.

Функция «\» скрывает фактический метод решения. На самом деле он выбирается в зависимости от типа матрицы А.

Например, можно явно указать, что матрица симметричная:

A = Symmetric(A)

После этого время решения почему-то несколько возрастает… :)

Скорее всего в исходном варианте использовалась LU-факторизация. Ее можно выполнить предварительно явно. Заменим строку с решением уравнения:

X = lufact(A) \ B

При этом функция «\» так же выберет необходимый метод с учетом типа аргументов, поступающих на вход. Первый аргумент будет иметь тип уже не матрицы, а Base.LinAlg.LU{Float64,Array{Float64,2}}, содержащий в себе факторизованную матрицу. Функция «\» это поймет.

Время решения в таком варианте не изменилось.

Если вынести факторизацию матрицы

Если lufact(A) вынести за тело цикла программа будет выглядеть так:

n = 1000; m = 100
s = 0.
A = [1-((i-j)/n)^2 for i in 1:n, j in 1:n] + 9*eye(n,n)
tic()
  luA = lufact(A)
  for k=1:m
    B = fill(Float64(k), n)
    X = luA \ B
    s += sum(X)
  end
toc()
println(s)

(Тут изменен также способ измерения времени на tic/toc).

Время решения уменьшилось до 0,04 с – в 3 раза.

Julia, стандартное средство (LLt)

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

X = cholfact(A) \ B

Время получается 0,98 с.

Julia, LAPACK (LU)

Из Джулии также можно напрямую вызвать процедуры LAPACK. Вариант с LU-разложением и вызовом процедуры gesv! будет выглядеть так:

import Base.LinAlg.LAPACK: gesv!
n = 1000; m = 100
s = 0.
A0 = [1-((i-j)/n)^2 for i in 1:n, j in 1:n] + 9*eye(n,n)
  @time for i = k:m
    B = fill(Float64(k), n)
    A = copy(A0)
    gesv!(A, B)
    s += sum(B)
  end
println(s)

Здесь матрицу А приходится восстанавливать из эталонной A0, т.к. gesv! ее изменяет. Отметим особенности:

  • написать A = A0 нельзя, т.к. в этом случае будет присвоен только указатель на матрицу;

  • функции, меняющие значения аргументов, в Джулии принято называть с восклицательным знаком.

Время решения – 1,3 с.

Julia, LAPACK (LLt)

Вариант с LU-разложением и вызовом процедуры LAPACK posv!:

import Base.LinAlg.LAPACK: posv!
n = 1000; m = 100
s = 0.
A0 = [1-((i-j)/n)^2 for i in 1:n, j in 1:n] + 9*eye(n,n)
  @time for k = 1:m
    B = fill(Float64(k), n)
    A = copy(A0)
    posv!('L', A, B)
    s += sum(B)
  end
println(s)

Время решения – 0,78 с.

Python

У меня установлен Python 3.5.2 из дистрибутива Anaconda 4.2.0 (64-bit), в который «из коробки» входят библиотеки NumPy 1.11.1 и ScyPy 0.18.1.

«Под капотом» у NumPy/SciPy – Intel MKL (интересно, почему это остается бесплатным?). Ожидаем отличную производительность.

Python + NumPy (LU)

Код:

import numpy as np
import time
n = 1000; m = 100
s = 0.
A = np.fromfunction(lambda i,j: 1-((i-j)/n)**2, (n,n)) + 9*np.eye(n)
t = time.time()
for k in range(1, m+1):
    B = np.full((n), float(k))
    X = np.linalg.solve(A, B)
    s += np.sum(X)
print("Время:", time.time() - t)
print("Сумма:", s)

Время работы – 0,93 с. Отличный результат для матрицы общего вида.

Модуль линейной алгебры у NumPy скромный. Готовая функция для решения СЛАУ здесь только для матриц общего вида, вызывающая ту же функцию LAPACK gesv.

Библиотека SciPy – намного богаче. В частности, в нее входят функции решения СЛАУ с матрицами различного вида.

Python + SciPy (LU)

Изменения:

    ...
    from scipy import linalg
    ...
    X = linalg.solve(A, B)
    ...

Время – 1,2 с. Странно, почему, больше чем у NumPy? Видимо там немного разные версии MKL вызываются.

Python + SciPy (LLt)

Изменения:

    ...
    from scipy import linalg
    ...
    X = linalg.solve(A, B, sym_pos = True)
    ...

Время – 0,98 с.

MatLAB-подобные СКМ

Octave (LU)

У меня установлен Octave 4.2.1. Для решения СЛАУ он использует OpenBLAS.

Код:

n = 1000; m = 100;
tic
A = zeros(n,n);
for i = 1:n
    for j = 1:n
        A(i,j) = 1-((i-j)/n)**2;
        if (i == j)
            A(i,j) = 10.;
        endif
    endfor
endfor
toc
s = 0.;
tic
for k = 1:m
    B = k * ones(n, 1);
    X = A \ B;
    s = s + sum(X);
endfor
toc
s

Цикл с решением СЛАУ выполнился за достойное время 1,4 с (все-таки здесь работает OpenBLAS), а вот предшествующая часть с заполнением матрицы отняла неприличные 7 секунд (!). Я в ступоре и прострации – что я делаю не так??? Во всех остальных «приличных» языках эта часть занимает пренебрежимо малое время.

SciLab (LU)

У меня установлен SciLab 6.0.0. Для решения СЛАУ он использует Intel MKL.

Код:

n = 1000; m = 100
A = zeros(n,n)
for i = 1:n
    for j = 1:n
        A(i,j) = 1-((i-j)/n)**2
        if i == j then
            A(i,j) = 10.
        end
    end
end
s = 0.
tic()
for k = 1:m
    B = k * ones(n, 1)
    X = A \ B
    s = s + sum(X)
end
disp(toc(), s)

Время – 1,2 с.

FreeMat (LU)

FreeMat – MatLAB-подобная СКМ, свободная и относительно легковесная, за что многие коллеги ее очень любят. К сожалению, проект не развивается с 2013 года. Доступна только в 32-битной версии.

tic;
n = 1000; m = 100;
A=zeros(n,n);
for i = 1:n
    for j = 1:n 
        A(i,j) = 1-((i-j)/n)^2;
        if (i == j) 
            A(i,j) = 10.;
        end
    end
end
disp(toc)
s=0.;
tic;
for k = 1:m
    B = k * ones(n, 1);
    X = A \ B;
    s = s + sum(X);
end
disp(toc, s)

Время решения задачи составило 96 с.

Время заполнения матрицы А – 15 с.

Ясно, что FreeMat для больших вычислений не предназначен.

Maxima

Maxima – свободная, мощная и очень старая система компьютерной алгебры (СКА). Не смотря на преклонный возраст проект развивается, периодически появляются новые версии. Текущая – 5.39.0. Я привык работать с графической оболочной wxMaxima. Когда нужно что-то быстро посчитать, сделать символьные преобразования – отличная вещь. Грешен: старушку Максиму я по-прежнему нежно люблю.

Конечно, Maxima не предназначена для вычислений большого объема. Однако, система – универсальная. В частности, в ней есть модуль линейной алгебры с вызовом библиотек LAPACK (в какой реализации – непонятно, видимо несовременной).

Решим нашу задачу и в Maxima, хотя, это скорей «по приколу» (ну и чтобы не забывать ее несколько специфичный язык):

load(lapack)$

n: 1000$  m: 100$
t0: elapsed_run_time()$
A: zeromatrix (n, n)$
for i:1 thru n do
  for j:1 thru n do
    A[i,j]: if (i#j) then float(1-((i-j)/n)**2) 
                     else float(10)$
B: zeromatrix (n,1)$
s: 0.$
t1: elapsed_run_time()$
t1-t0;
for k:1 thru m do (
    B: B + 1.,
    X: dgesv (A, B),
    s: s + sum(X[i,1], i, 1, n)
)$
elapsed_run_time()-t1;
s;

Время заполнения матрицы А составило аж 70 с.

Время решения задачи… спустя несколько часов я не дождался. :) На малых объемах данных приведенное решение работает правильно, но при n=1000 и m=100 – все о-очень долго.

Сравнение и заключение

Мы решали систему уравнений с матрицей вещественных чисел двойной точности, размером 1000 х 1000, повторяя решение в цикле 100 раз.

Общие сравнительные результаты выглядят следующим образом:




Это иллюстрация очевидного факта: скорость решения СЛАУ разными средствами может колоссально разниться – в разы, на порядки…

Не предназначенная для решения подобных задач Maxima работала так долго, что я решения не дождался. Ей можно простить. Неоптимизированный FreeMat, позиционируемый как аналог MatLAB, тоже рекордно забуксовал. Разочаровал GFortran с классическими BLAS/LAPACK без распараллеливания и оптимизации.

Если откинуть перечисленных «аутсайдеров» получается следующая картина:




Абсолютный лидер – коммерческий Intel Fortran + MKL (LLt). Время решения задачи у него заняло всего 0,33 с. Инженеры Intel едят хлеб не зря. Если надо создать по-настоящему быструю реализацию на Fortran чего-либо серьезного, покупка лицензии Intel – правильное решение.

Догоняющий лидера свободный аналог – GFortran + OpenBLAS (LLt) с временем 0,51 с. Единственный недостаток этой реализации – огромный exe-файл (26 МБ). Это несколько напрягает. Надо будет попробовать с Фортраном другие оптимизированные реализации BLAS/LAPACK: например, GotoBLAS2 или ScaLAPACK.

Библиотека LAIPE2 от профессора Ляо как-то разочаровала.

В протестированных современных высокоуровневых средствах – языках Julia, Python и СКМ Octave и SciLAB в качестве решателей СЛАУ используются скомпилированные оптимизированные библиотеки Intel MKL или OpenBLAS. Поэтому результаты у них получились весьма достойными, а между собой – близкими. Выигрышней всего смотрится Julia + LAPACK (LLt) – 0,78 с. Единственное, неприятно удивил Octave, который собирал исходную матрицу неприличные 7 секунд. Может, я что-то сделал не так?

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

Очевидная истина – решатель СЛАУ надо выбирать соответственно типа матрицы. В нашем примере была симметрическая положительно определенная матрица. Все решения с \(LL^T\)-разложением для нее были заметно быстрее аналогов с \(LU\)-разложением, подходящим для матриц общего вида. Для разряженных матриц следует использовать соответствующие решатели. Если можно сделать разряженную матрицу ленточной с малой шириной ленты – будет совсем хорошо. Аналогичный тест для разряженных матриц – надеюсь, еще проведем. Там свой «мир» решателей, изучить и освоить его – было бы очень полезно и интересно.

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

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

2 комментария:

  1. Благодарю, Евгений! Весьма полезный сравнительный анализ. ИМХО, уже на подходе задачи, в которых результаты вашей работы будут определять выбор инструментария программирования. С уважением, И.Б.

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

    ОтветитьУдалить