Порівняння продуктивності GPU-розрахунків на Python і C


Python має ряд привабливих переваг, до яких відноситься простота реалізації програмних рішень, наочність і лаконічність коду, наявність великої кількості бібліотек та численного активного ком'юніті. У той же час, відома всім повільність пітона часто обмежує його придатність для «важких» обчислень. Для ряду задач можна добитися істотного прискорення розрахунків шляхом використання технології CUDA для паралельних обчислень на GPU. Мета цього невеликого дослідження — аналіз можливостей ефективного використання python для розрахунків на GPU і порівняння продуктивності різних python-рішень з реалізацією на C.

Введення
За родом діяльності я часто займаюся завданнями чисельного моделювання. У багатьох випадках ми використовуємо технологію CUDA для прискорення розрахунків за рахунок використання можливостей паралельних обчислень на GPU, при цьому програми пишуться на мові C. У той же час, в деяких випадках хотілося б мати можливість реалізовувати розрахунки на пітоні, тому що це зручно, швидко, гнучко, лаконічно, молодіжно. При цьому, однак, дуже важливо не втрачати в швидкості виконання програм, оскільки деякі розрахунки можуть займати від декількох годин до доби. У статті продемонстровано використання python-бібліотек Numba і PyCUDA для реалізації паралельних обчислень на GPU і наведені результати порівняння їх продуктивності на тестовій задачі.

Тестова задача
Вибір тестового завдання і умов тестування був обумовлений характером реальних завдань, для вирішення яких надалі планується використання python. При цьому була обрана максимально просте завдання, а саме завдання рішення двовимірного рівняння теплопровідності з допомогою явній скінченно-різницевої схеми. Задача розглядалася на квадратній області з заданими значеннями температури на кордонах. Кількість вузлів розрахункової сітки по x і y однаково і дорівнює n. На малюнку на початку статті показано стале рішення тестової задачі.

Алгоритм і умови тестування
Алгоритм розв'язання задачі можна представити наступним псевдокодом:

Ініціалізувати масиви u0 і u, зберігають значення рішення на поточному та наступному часовому кроці
включити таймер
виділити пам'ять на девайсі і скопіювати масиви u0 і u з хоста на девайс
for(i = 0; i < nstp/2; i++)
обчислити значення u на наступному кроці за часом, знаючи значення u0 на поточному часовому кроці
обчислити значення u0 на наступному кроці за часом, знаючи значення u на поточному часовому кроці
скопіювати результат (u0) з девайса на хост
вимкнути таймер і вивести результат

У всіх тестах представлених нижче число вузлів квадратної сітки в кожному напрямку (n) варіювалося від 512 до 4096, а nstp = 5000.

Програмне забезпечення
Тестування проводилося на персональному комп'ютері:

Intel® Core(TM)2 Quad CPU Q9650 @ 3.00 GHz, 8 Gb ОЗУ
GPU: Nvidia GTX 580
Операційна система: Ubuntu 16.04 LTS з встановленою CUDA 7.5
Реалізація на C
Всі подальші python-реалізації порівнювалися з результатами, отриманими за допомогою програми на C, описаної в цьому розділі.

Код на C
#include < stdio.h>
#include <time.h>
#include <cuda_runtime.h>

#define BLOCK_SIZE 16

int N = 512; //1024; //2048;
double L = 1.0;
double h = L/N;
double h2 = h*h;
double tau = 0.1*h2;
double th2 = tau/h2;
int nstp = 5000;

__constant__ int dN;
__constant__ double dth2;

__global__ void NextStpGPU(double *u0, double *u)
{
int i = blockDim.x * blockIdx.x + threadIdx.x;
int j = blockDim.y * blockIdx.y + threadIdx.y;

double uim1, uip1, ujm1, ujp1, u00, d2x, d2y;

if (i > 0)
uim1 = u0[(i - 1) + dN*j];
else
uim1 = 0.0;

if (i < dN - 1)
uip1 = u0[(i + 1) + dN*j];
else
uip1 = 0.0;

if (j > 0)
ujm1 = u0[i + dN*(j - 1)];
else
ujm1 = 0.0;

if (j < dN - 1)
ujp1 = u0[i + dN*(j + 1)];
else
ujp1 = 1.0;

u00 = u0[i + dN*j];
d2x = (uim1 - 2.0*u00 + uip1);
d2y = (ujm1 - 2.0*u00 + ujp1);
u[i + dN*j] = u00 + dth2*(d2x + d2y);

}

int main(void)
{
size_t size = N * N * sizeof(double);
double pStart, pEnd, pD;
int i;

double *h_u0 = (double *)malloc(size);
double *h_u1 = (double *)malloc(size);

for (i = 0; i < N*N; ++i)
{
h_u0[i] = 0.0;
h_u1[i] = 0.0;
}

pStart = (double) clock();

double *d_u0 = NULL;
cudaMalloc((void **)&d_u0, size);
double *d_u1 = NULL;
cudaMalloc((void **)&d_u1, size);

cudaMemcpy(d_u0, h_u0, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_u1, h_u1, size, cudaMemcpyHostToDevice);

cudaMemcpyToSymbol(dN, &N, sizeof(int));
cudaMemcpyToSymbol(dth2, &th2, sizeof(double));

dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE, 1);
dim3 dimGrid(N/BLOCK_SIZE, N/BLOCK_SIZE, 1);

for(i = 0; i < int(nstp/2); i++)
{
NextStpGPU<<<dimGrid,dimBlock>>>(d_u0, d_u1);
NextStpGPU<<<dimGrid,dimBlock>>>(d_u1, d_u0);
}

cudaMemcpy(h_u0, d_u0, size, cudaMemcpyDeviceToHost);
cudaFree(d_u0);
cudaFree(d_u1);

pEnd = (double) clock();
pD = (pEnd-pStart)/CLOCKS_PER_SEC;
printf("Calculation time on GPU = %f sec\n", pD);
}

free(h_u0);
free(h_u1);

return 0;
}


Розрахунки показують, що для N = 512 час виконання C-програми розпаралелиного на GPU становить 0.27 секунд проти 33.06 секунд для послідовної реалізації алгоритму на CPU. Тобто прискорення CPU/GPU становить близько 120 разів. З ростом N величина прискорення не убуває.

Python with Numba
Бібліотека Numba надає можливість jit (just-in-time) компіляції коду на пітоні в байт-код порівнянний по продуктивності з C або Fortran кодом. Numba підтримує компіляцію і запуск python-коду не лише на CPU, але і на GPU, при цьому стиль і вигляд програми, що використовує бібліотеку Numba, залишається чисто питоновским.

Ось як виглядає в цьому випадку наша програма

from numba import cuda
import numpy as np
import matplotlib.pyplot as plt
from time import time

n = 512
blockdim = 16, 16
griddim = int(n/blockdim[0]), int(n/blockdim[1])

L = 1.
h = L/n
dt = 0.1*h**2
nstp = 5000

@cuda.jit("void(float64[:,:], float64[:,:])")
def nextstp_gpu(u0, u):
i,j = cuda.grid(2)

if i > 0:
uim1 = u0[i-1,j]
else:
uim1 = 0.
if i < n-1:
uip1 = u0[i+1,j]
else:
uip1 = 0.
if j > 0:
ujm1 = u0[i,j-1]
else:
ujm1 = 0.
if j < n-1:
ujp1 = u0[i,j+1]
else:
ujp1 = 1.

d2x = (uim1 - 2.*u0[i,j] + uip1)
d2y = (ujm1 - 2.*u0[i,j] + ujp1)
u[i,j] = u0[i,j] + (dt/h**2)*(d2x + d2y)

u0 = np.full((n,n), 0., dtype = np.float64)
u = np.full((n,n), 0., dtype = np.float64)

st = time()

d_u0 = cuda.to_device(u0)
d_u = cuda.to_device(u)
for i in xrange(0, int(nstp/2)):
nextstp_gpu[griddim,blockdim](d_u0, d_u)
nextstp_gpu[griddim,blockdim](d_u, d_u0)

cuda.synchronize()
u0 = d_u0.copy_to_host()
print 'time on GPU = ', time() - st


Зазначимо тут декілька приємних особливостей. По-перше, ця реалізація набагато коротше і наочніше. Тут ми використовували двовимірні масиви, що робить код набагато більш зручним для читання. По-друге, якщо C-реалізації від нас було потрібно передати всі константи (наприклад N) шляхом виконання функцій на зразок
cudaMemcpyToSymbol(dN, &N, sizeof(int));
то тут ми просто користуємося глобальними змінними як у звичайній python-функції. Ну і нарешті, реалізація не вимагає жодних знань мови C і архітектури GPU.

Цей код легко переписати і для випадку використання одновимірних масивів розміру n*n, як буде показано далі, це істотно впливає на швидкість виконання.

Тут представлений такий код

from numba import cuda
import numpy as np
import matplotlib.pyplot as plt
from time import time

n = 512
blockdim = 16, 16
griddim = int(n/blockdim[0]), int(n/blockdim[1])

L = 1.
h = L/n
dt = 0.1*h**2
nstp = 5000

@cuda.jit("void(float64[:], float64[:])")
def nextstp_gpu(u0, u):
i,j = cuda.grid(2)

u00 = u0[i + n*j]
if i > 0:
uim1 = u0[i-1 + n*j]
else:
uim1 = 0.
if i < n-1:
uip1 = u0[i+1 + n*j]
else:
uip1 = 0.
if j > 0:
ujm1 = u0[i + n*(j-1)]
else:
ujm1 = 0.
if j < n-1:
ujp1 = u0[i + n*(j+1)]
else:
ujp1 = 1.

d2x = (uim1 - 2.*u00 + uip1)
d2y = (ujm1 - 2.*u00 + ujp1)
u[i + n*j] = u00 + (dt/h/h)*(d2x + d2y)

u0 = np.full(n*n, 0., dtype = np.float64)
u = np.full(n*n, 0., dtype = np.float64)

st = time()

d_u0 = cuda.to_device(u0)
d_u = cuda.to_device(u)
for i in xrange(0, int(nstp/2)):
nextstp_gpu[griddim,blockdim](d_u0, d_u)
nextstp_gpu[griddim,blockdim](d_u, d_u0)

cuda.synchronize()
u0 = d_u0.copy_to_host()
print 'time on GPU = ', time() - st


PyCUDA
Другий з протестованих python-бібліотек була бібліотека PyCUDA. На відміну від Numba тут від розробника потрібно написати код ядра на C, тому без знання цієї мови не обійтися. З іншого боку, крім власне ядра на C нічого писати не треба.

Код з використанням PyCUDA виходить таким

import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
import numpy as np
import matplotlib.pyplot as plt
from time import time

n = 512
blockdim = 16, 16
griddim = int(n/blockdim[0]), int(n/blockdim[1])

L = 1.
h = L/n
dt = 0.1*h**2
nstp = 5000

mod = SourceModule("""
__global__ void NextStpGPU(int* dN, double* dth2, double *u0, double *u)
{
int i = blockDim.x * blockIdx.x + threadIdx.x;
int j = blockDim.y * blockIdx.y + threadIdx.y;

double uim1, uip1, ujm1, ujp1, u00, d2x, d2y;
uim1 = exp(-10.0);
if (i > 0)
uim1 = u0[(i - 1) + dN[0]*j];
else
uim1 = 0.0;

if (i < dN[0] - 1)
uip1 = u0[(i + 1) + dN[0]*j];
else
uip1 = 0.0;

if (j > 0)
ujm1 = u0[i + dN[0]*(j - 1)];
else
ujm1 = 0.0;

if (j < dN[0] - 1)
ujp1 = u0[i + dN[0]*(j + 1)];
else
ujp1 = 1.0;

u00 = u0[i + dN[0]*j];
d2x = (uim1 - 2.0*u00 + uip1);
d2y = (ujm1 - 2.0*u00 + ujp1);
u[i + dN[0]*j] = u00 + dth2[0]*(d2x + d2y);
}
""")

u0 = np.full(n*n, 0., dtype = np.float64)
u = np.full(n*n, 0., dtype = np.float64)
nn = np.full(1, n, dtype = np.int64)
th2 = np.full(1, dt/h/h, dtype = np.float64)

st = time()

u0_gpu = cuda.to_device(u0)
u_gpu = cuda.to_device(u)
n_gpu = cuda.to_device(nn)
th2_gpu = cuda.to_device(th2)

func = mod.get_function("NextStpGPU")
for i in xrange(0, int(nstp/2)):
func(n_gpu, th2_gpu, u0_gpu, u_gpu, 
block=(blockdim[0],blockdim[1],1),grid=(griddim[0],griddim[1],1))
func(n_gpu, th2_gpu, u_gpu, u0_gpu, 
block=(blockdim[0],blockdim[1],1),grid=(griddim[0],griddim[1],1))

u0 = cuda.from_device_like(u0_gpu, u0)

print 'time on GPU = ', time() - st


Виглядає це все як чистий пітон за винятком локальної вставки C-коду.

Порівняння продуктивності
На Рисунку 1 і в Таблиці 1 наведено залежності часу виконання тестової програми (у секундах) від розміру сітки n, отримані при запуску C-коду (крива CUDA C) і python-реалізацій з бібліотекою Numba і двовимірними масивами (Numba 2DArr), з бібліотекою Numba і одновимірними масивами (Numba 1DArr), з бібліотекою PyCUDA (крива PyCUDA).

Рисунок 1
Малюнок 1
Таблиця 1
n Cuda C Numba 2DArr Numba 1DArr PyCUDA
512 0.25 0.8 0.66 0.216
1024 0.77 3.26 1.03 0.808
2048 2.73 12.23 4.07 2.87
3073 6.1 27.3 9.12 6.6
4096 11.05 55.88 16.27 12.02
На Рисунку 2 наведено відношення часу виконання різних python-реалізацій до часу виконання, C-коду. Як видно з малюнків, самою повільною, з розглянутих, є реалізація з допомогою бібліотеки Numba з використанням двовимірних масивів. При цьому, цей підхід є найбільш наочним і простим. Цікаво, що розгортка двовимірних масивів в одномірні призводить до приблизно триразового прискоренню коду. Найшвидшим рішенням виявилося використання бібліотеки PyCUDA. У той же час, як зазначалося вище, використання цієї бібліотеки декілька більш трудомістким, оскільки вимагає написання ядра на C. Однак витрати окупаються і швидкість виконання такої python-програми лише на 5-8% менше ніж програми, написаної повністю на C.

Рисунок 2
Малюнок 2

Висновки
Чудес не буває і самі прості і наочні рішення виявляються одночасно і самими повільними. У той же час, наявні бібліотеки дозволяють досягти швидкості виконання python-програм, порівнянної зі швидкістю виконання чистого C-коду. Існуючі бібліотеки дають розробнику вибір між більш і менш високорівневими рішеннями. Однак цей вибір завжди є компроміс між швидкістю розробки і швидкістю виконання програми.

Посилання
» Документація бібліотеки Numba
» Приклади використання Numba
» Документація бібліотеки PyCUDA
» Приклади використання PyCUDA
Джерело: Хабрахабр

0 коментарів

Тільки зареєстровані та авторизовані користувачі можуть залишати коментарі.