mariakasimceva19 мар в 08:15

Почему Python + Numba обгоняет C? Эксперимент с алгоритмом прогонки

Средний

4 мин

6.7K

C * Python * Алгоритмы *

Аналитика

Из песочницы

Комментарии 12

miekrudakov 19 мар в 08:25

Выводы не совсем верные

Скорость Numpy зависит от реализации
Скорость Numba зависит от реализации
Скорость C++ зависит от реализации
Правильный вывод: измерять скорость (правильно измерять!) и не бояться JIT компиляции

tenzink 19 мар в 09:23

прогрев (для исключения влияния кэша и JIT-компиляции)

Не кажется, что это немного тепличные условия? В реальном приложении никакого прогрева вероятней всего не будет и приближенным к реальности будет вариант первого запуска с загрузкой всего и вся. Интересно, насколько там большая дисперсия получается, если не делать прогрев

Anton_Menshov 19 мар в 14:00

Если вы пишете линейную алгебру на C/C++ и НЕ пользуетесь высокопроизводительной BLAS/LAPACK библиотекой - есть только три варианта:

1. линейная алгебра составляет менее 1% вычислений в вашей задаче - и тогда все равно какая у нее производительность.
2. вам [пока] не стоит писать программы для линейной алгебры ни для чего, кроме как для обучения себя.
3. вы сами пишете свой BLAS\LAPACK по какой-то экзотической причине. вполне возможно обоснованной.

А так, статья хорошая: лучше использовать хорошие примитивы из numba, чем писать линейку на С\С++ вручную без понимания того, как и когда это нужно делать хорошо.

Pand5461 19 мар в 15:48

Для трёхдиагональных систем библиотеки мало что дадут и могут быть даже медленнее ручного кода, если, например, дают гарантию устойчивости для любой системы. Прогонка всё равно в скорость доступа к памяти упирается, там же алгоритм последовательный и не векторизуемый.
Бонусом - если писать вручную, то можно соптимизировать для специальных типов матриц (например, если на диагонали все значения одинаковые, можно их не хранить в явном виде). А в задачах, где тридиагональные системы возникают, они как раз имеют специальный вид часто.

Anton_Menshov 19 мар в 17:13

В такой постановке проблемы вы правы. Релевантно:
- How does LAPACK solve tridiagonal systems and why - partial pivoting - это о гарантии устойчивости и то, что LAPACK делает, а вам может и не надо.
- Об оптимальности Томаса для трехдиагональной системы
- Parallel TDMA - для оооооочень больших трехдиагональных систем

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

Pand5461 20 мар в 10:03

А я вот так и узнал, что в LAPACK алгоритм чуть сложнее стандартной прогонки, как в первой ссылке - обнаружил, что своя реализация по учебнику быстрее библиотечной процентов на 10-20. Но в реальных задачах часто гарантируется, например, диагональное преобладание - а этого достаточно для устойчивости даже без выбора главного элемента.
С параллелизмом есть вопросы. Не так часто возникают задачи, где в основе будет гигантская трёхдиагональная система. По ссылке, впрочем, говорится также про оптимизированный солвер для большого числа трёхдиагональных систем, это как раз реально.

или уж действительно выжимать все соки из своей имплементации (например, не аллоцировать рабочие массивы каждый раз, как в статье)

Тут полностью согласен. В Julia тут была бы стандартная идиома - функция thomas_solver!(rhp, ld, d, ud), которая перезаписывает d и ud, а в rhp записывает решение, и thomas_solver(rhp, ld, d, ud), которая создаёт временные массивы и с ними запускает предыдущую.

laptev_am 19 мар в 18:32

Я уже лет 15 математикой не занимался, так что мне пришлось прибегнуть к помощи нейронки. Но результаты от этого не сильно пострадали.

Железо и компилятор

$ lscpu | grep "Model name"
Model name: AMD Ryzen 7 3700U with Radeon Vega Mobile Gfx
$ sudo dmidecode --type memory | grep -i speed
	Speed: 2400 MT/s
	Configured Memory Speed: 2400 MT/s
	Speed: 2400 MT/s
	Configured Memory Speed: 2400 MT/s
$ gfortran --version
GNU Fortran (Ubuntu 13.3.0-6ubuntu2~24.04.1) 13.3.0
Copyright (C) 2023 Free Software Foundation, Inc.

листинг

program thomas_test
    implicit none
    integer, parameter :: n = 10**7  ! 10 миллионов элементов
    real(8), allocatable :: a(:), b(:), c(:), d(:), x(:), cp(:), dp(:)
    real(8) :: m, t1, t2
    integer :: i

    call cpu_time(t1)

    print *, "Выделение памяти для n =", n
    allocate(a(n-1), b(n), c(n-1), d(n), x(n), cp(n-1), dp(n))

    ! Наполним данными, чтобы процессор не скучал
    a = 1.0d0; b = 4.0d0; c = 1.0d0; d = 10.0d0

    print *, "Поехали! Считаем прогонку..."

    ! --- Тот самый алгоритм ---
    cp(1) = c(1) / b(1)
    dp(1) = d(1) / b(1)
    do i = 2, n
        m = 1.0d0 / (b(i) - a(i-1) * cp(i-1))
        if (i < n) cp(i) = c(i) * m
        dp(i) = (d(i) - a(i-1) * dp(i-1)) * m
    end do

    x(n) = dp(n)
    do i = n-1, 1, -1
        x(i) = dp(i) - cp(i) * x(i+1)
    end do
    ! -------------------------

    call cpu_time(t2)
    print "(A, F10.6, A)", " Время выполнения: ", t2 - t1, " сек."
    print *, "Контрольное значение x(1):", x(1)

    deallocate(a, b, c, d, x, cp, dp)
end program thomas_test

результаты

$ gfortran -O3 -march=native -ffast-math thomas_test.f90 -o thomas_test
$ ./thomas_test 
 Выделение памяти для n =    10000000
 Поехали! Считаем прогонку...
 Время выполнения:   0.371780 сек.
 Контрольное значение x(1):   2.1132486540518713

Вывода не будет...

muwrec 20 мар в 04:59

Если C код медленнее питона, то к кого-то руки из ж...

1kovand1 20 мар в 07:49

Numba подстраивается под конкретный процессор, тогда как бинарник C обычно компилируется с некоторым усреднённым набором инструкций.

Но в данном случае gcc же явно вызывался с параметром -march=native. Это буквально по определению значит использовать все доступные наборы инструкций:

Using -march=native enables all instruction subsets supported by the local machine (hence the result might not run on different machines).

Pand5461 20 мар в 14:08

Намеряны какие-то очень не те времена. У меня получается с numpy примерно 20 секунд на 10М уравнений, с numba примерно 200 мс, на Си с -Ofast и -O3 130 мс (AMD Ryzen 5 3500U). И объяснений, почему вдруг Numba быстрее Си, уже не требуется.

Anton_Menshov 20 мар в 16:03

несмотря на то, что вы правы, все же нумбу в этом случае нужно сравнивать со включенной fastmath=True опцией. Или не использовать -Ofast

Pand5461 20 мар в 16:26

Видимо, непонятно написал - имеется в виду, что -O3 и -Ofast не отличаются по производительности.
В Julia ещё проверял, там есть эффект в 10-15% от отключения проверок границ массива (fastmath, как и в Си, на скорость не повлиял). Если для numba есть такой флажок, то он, скорее всего, тоже с 200 до 170 мс время скинет.

Зарегистрируйтесь на Хабре, чтобы оставить комментарий