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

[miekrudakov]

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

1. Скорость Numpy зависит от реализации

2. Скорость Numba зависит от реализации

3. Скорость C++ зависит от реализации

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

[tenzink]

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

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

[Anton_Menshov]

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

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

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

[Pand5461]

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

[Anton_Menshov]

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

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

[laptev_am]

Я уже лет 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  

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