Русский
Русский
English
Статистика
Реклама

Чистый Cython VS nvc жжем металлические пластины на GPU для сравнения скорости

image

image
Будем греть металлические пластины на GPU

Все знают, что Python не блещет скоростью сам по себе. На мой взгляд язык прекрасен своей читабельностью, но основная ниша его применения там, где вы большую часть времени ожидаете ввода/вывода каких-то данных. Условно, вы можете написать суперпроизводительный код на Rust или С, но 99% времени он будет просто ждать.

Тем не менее, Python прекрасен еще и как высокоуровневый синтаксический клей. В этом случае, его неторопливая интерпретируемая часть вызывает быстродействующий код, написанный на компилируемых языках программирования. Обычно для этого используются такие традиционные библиотеки как NumPy.

Но мы пойдем чуть дальше попробуем распараллелить вычисления на CUDA и задействуем странный, но работающий гибрид C++, stdpar и компилятора nvc++ от Nvidia. Ну и заодно попробуем оценить быстродействие. Возьмем две задачи: сортировку чисел и метод Якоби, которым будем рассчитывать нагрев металлической пластины.

Вызываем C++ из Python


Наш код сортировки будет иметь следующий вид:

# distutils: language=c++from libcpp.algorithm cimport sortdef cppsort(int[:] x):    sort(&x[0], &x[-1] + 1)

В первой строчке мы явно указываем, что Cython должен использовать C++, а не дефолтный C. Во второй мы импортируем функцию сортировки из C++, а дальше следует сама логика. Помещаем код в файл cppsort.pyx. Обратите внимание, что расширение отличается от привычного py, так как мы будем его компилировать или выполнять cythonize в терминологии Cython.

Компиляцию можно выполнить вручную или включить в setup.py, где мы полноценно описываем подготовку нашего окружения.

В setup.py это выглядит примерно так:

from setuptools import setupfrom Cython.Build import cythonizesetup(    ext_modules = cythonize("cppsort.pyx"))

Но мы можем и просто выполнить это в командной строке:

cythonize -i cppsort.pyx

Под капотом произойдет примерно следующее:

image

  1. Cython транслирует python код в C++ и сгенерирует валидный cppsort.cpp.
  2. C++ компилятор (в данном случае g++) компилирует C++ код в Python extension module
  3. Python extension module импортируется в код как обычный питоновский модуль.

После компиляции можем импортировать и сразу протестировать сортировку:

from cppsort import cppsortimport numpy as npx = np.array([4, 3, 2, 1], dtype="int32")print(x)cppsort(x)print(x)

Массив [4, 3, 2, 1] успешно отсортируется в [1, 2, 3, 4] с помощью C++ std::sort.

А давайте на GPU?


Стандартные библиотечные алгоритмы C++ могут вызываться с указанием аргумента parallel execution policy. Этот аргумент говорит компилятору о том, что вы хотите раскидать алгоритм на параллельные процессы.

При этом C++ имеет несколько вариантов этой политики:

  1. std::execution::seq: Последовательное выполнение. Параллельность запрещена.
  2. std::execution::unseq: Векторизированное выполнение в рамках вызвавшего потока.
  3. std::execution::par: Параллельное выполнение в одном и более потоках.
  4. std::execution::par_unseq: Параллельное выполнение в одном и более потоках. Каждый поток будет по возможности векторизирован.

При этом вы сами должны следить за race condition и deadlock. Стандартный компилятор g++ постарается распределить вычисления на ядра CPU. Но мы можем взять проприетарный компилятор от Nvidia nvc++ и скормить ему опцию "-stdpar". stdpar это C++ Standard Parallelism от Nvidia с выполнением параллельного кода на GPU.

Перепишем код, с учетом необходимости создавать локальную копию массива, так как GPU не сможет получить доступ к массиву, расположенному в рамках NumPy.

from libcpp.algorithm cimport sort, copy_nfrom libcpp.vector cimport vectorfrom libcpp.execution cimport pardef cppsort(int[:] x):    cdef vector[int] temp    temp.resize(len(x))    copy_n(&x[0], len(x), temp.begin())    sort(par, temp.begin(), temp.end())    copy_n(temp.begin(), len(x), &x[0])

image

Теперь это нужно снова скомпилировать, но уже с использованием nvc++. В этот раз напишем нормальный setup.py и вызовем его:

CC=nvc++ python setup.py build_ext --inplace

Импортируем в код и пробуем вызвать:

x = np.array([4, 3, 2, 1], dtype="int32")cppsort(x) # этот кусок выполняется на GPU

Производительность


Традиционно GPU хороши там, где есть много однотипных легковесных вычислений. Тяжелые задачи одиночные задачи GPU не подойдут. Более того, стоит учитывать объем ваших данных. Если у вас немного данных, то вы получите большой оверхед на процесс распараллеливания, ввод/вывод на между CPU и GPU. В итоге, такой код скорее всего наиболее эффективно будет выполняться на чистом CPU, иногда даже в пределах одного ядра, если данных совсем немного. Но на больших массивах GPU однозначно будет впереди.

Вот тут есть отличное сравнение сортировки. За единицу брали скорость NumPy, а затем считали кратность прироста скорости в каждом методе относительно него.


image

Как мы видим, гипотеза о том, что выигрыш в использовании GPU на большом количестве данных растет, подтверждается.

Вычисляем нагрев пластины


Возьмем задачу более приближенную к реальному инженерному моделированию вычислениям по методу Якоби. В частности, они отлично подходят для вычисления температурных процессов в 2D-пространстве.

image

Код для симуляции
"""simulates heat equation on rectangle returning a heat map at a number of timesboundary and initial conditions are 0, source represents burner on a stoveThis program is based on the script FEniCS tutorial demo program: Diffusion of a Gaussian hill.       u'= Laplace(u) + f  in a square domain  u = u_D = 0            on the boundary  u = u_0 = 0            at t = 0  u_D = f = stove burner flameThis program succesfully runs in the fenics docker image, see the book Solving PDEs in Python.to animate: convert -delay 4 -loop 100 heatequation10*.png heatstovelinn.gifto crop:convert heatstovelinn.gif -coalesce -repage 0x0 -crop 810x810+95+15 +repage heatstovelin.gif"""from fenics import *import timeimport matplotlib.pyplot as pltfrom matplotlib import cm# Create mesh and define function spacenx = ny = 100mesh = RectangleMesh(Point(-2, -2), Point(2, 2), nx, ny)V = FunctionSpace(mesh, 'P', 1)# Define boundary, source, initialdef boundary(x, on_boundary):    return on_boundarybc = DirichletBC(V, Constant(0), boundary)u_0 = interpolate(Constant(0), V)f = Expression('exp(-sqrt(pow((a*pow(x[0], 2) + a*pow(x[1], 2)-a*1),2)))', degree=2, a=5) #steep guassian centred on the unit spherefinal_time = 0.035num_pics = 72for i in range(num_pics):    T =   final_time*(i+1.0)/(num_pics+1)      #solve time even space    #T = final_time*1.1**(i-num_pics+1)        #solve time log  space    num_steps = 30    dt = T / num_steps # time step size    # Define variational problem    u = TrialFunction(V)    v = TestFunction(V)    F = u*v*dx + dt*dot(grad(u), grad(v))*dx - (u_0 + dt*f)*v*dx    a, L = lhs(F), rhs(F)    # Time-stepping    u = Function(V)    t = 0    for n in range(num_steps):        t += dt              #step        solve(a == L, u, bc) #solve        u_0.assign(u)        #update    #plot solution    plot(u,cmap=cm.hot,vmin=0,vmax=0.07)    plt.axis('off')    plt.savefig('heatequation10%s.png'%(i+10),figsize=(8, 8), dpi=220,bbox_inches='tight', pad_inches=0,transparent=True)


Напишем аналогичный солвер на Cython для последующей компиляции по CUDA:

# distutils: language=c++# cython: cdivision=Truefrom libcpp.algorithm cimport swapfrom libcpp.vector cimport vectorfrom libcpp cimport bool, floatfrom libc.math cimport fabsfrom algorithm cimport for_each, any_of, copyfrom execution cimport par, seq cdef cppclass avg:    float *T1    float *T2    int M, N     avg(float* T1, float *T2, int M, int N):        this.T1, this.T2, this.M, this.N = T1, T2, M, N    inline void call "operator()"(int i):        if (i % this.N != 0 and i % this.N != this.N-1):            this.T2[i] = (                this.T1[i-this.N] + this.T1[i+this.N] + this.T1[i-1] + this.T1[i+1]) / 4.0cdef cppclass converged:    float *T1    float *T2    float max_diff     converged(float* T1, float *T2, float max_diff):        this.T1, this.T2, this.max_diff = T1, T2, max_diff     inline bool call "operator()"(int i):        return fabs(this.T2[i] - this.T1[i]) > this.max_diff def jacobi_solver(float[:, :] data, float max_diff, int max_iter=10_000):    M, N  = data.shape[0], data.shape[1]    cdef vector[float] local    cdef vector[float] temp    local.resize(M*N)    temp.resize(M*N)    cdef vector[int] indices = range(N+1, (M-1)*N-1)    copy(seq, &data[0, 0], &data[-1, -1], local.begin())    copy(par, local.begin(), local.end(), temp.begin())    cdef int iterations = 0    cdef float* T1 = local.data()    cdef float* T2 = temp.data()     keep_going = True    while keep_going and iterations < max_iter:        iterations += 1        for_each(par, indices.begin(), indices.end(), avg(T1, T2, M, N))        keep_going = any_of(par, indices.begin(), indices.end(), converged(T1, T2, max_diff))        swap(T1, T2)     if (T2 == local.data()):        copy(seq, local.begin(), local.end(), &data[0, 0])    else:        copy(seq, temp.begin(), temp.end(), &data[0, 0])    return iterations

image

image

В итоге отрыв GPU получается еще более существенным.

Минусы


  1. Написание такого кода несколько сложнее, чем чистого варианта на Python и требует понимания принципов работы параллельных вычислений на GPU.
  2. Требуется копирование данных в отдельный массив для передачи на GPU, куда видеокарта не имеет доступа. Это может быть проблемой при работе с очень большими массивами.

Источник: habr.com
К списку статей
Опубликовано: 11.01.2021 20:23:30
0

Сейчас читают

Комментариев (0)
Имя
Электронная почта

Блог компании ruvds.com

Python

Программирование

Cyton

Nvc++

Ruvds_статьи

Gpu

Nvidia

Категории

Последние комментарии

  • Имя: Макс
    24.08.2022 | 11:28
    Я разраб в IT компании, работаю на арбитражную команду. Мы работаем с приламы и сайтами, при работе замечаются постоянные баны и лаги. Пацаны посоветовали сервис по анализу исходного кода,https://app Подробнее..
  • Имя: 9055410337
    20.08.2022 | 17:41
    поможем пишите в телеграм Подробнее..
  • Имя: sabbat
    17.08.2022 | 20:42
    Охренеть.. это просто шикарная статья, феноменально круто. Большое спасибо за разбор! Надеюсь как-нибудь с тобой связаться для обсуждений чего-либо) Подробнее..
  • Имя: Мария
    09.08.2022 | 14:44
    Добрый день. Если обладаете такой информацией, то подскажите, пожалуйста, где можно найти много-много материала по Yggdrasil и его уязвимостях для написания диплома? Благодарю. Подробнее..
© 2006-2024, personeltest.ru