Понимание и оптимизация потоков, блоков и сеток в pyCUDA

Question

Я очень новичок в программировании GPU и pyCUDA и имею довольно существенный пробел в своих знаниях. Я потратил немало времени на поиск SO, глядя на пример кода и чтение сопроводительной документации для CUDA/pyCUDA, но не нашел много разнообразия в объяснениях и не мог окунуться в некоторые моменты.

У меня возникли проблемы с правильным определением размеров блока и сетки. Код, я в настоящее время работает следующим образом, и стремится сделать поэлементное умножение массива a поплавком b:

from __future__ import division 
import pycuda.gpuarray as gpuarray 
import pycuda.driver as cuda 
import pycuda.autoinit 
from pycuda.compiler import SourceModule 
import numpy as np 

rows = 256 
cols = 10 
a = np.ones((rows, cols), dtype=np.float32) 
a_gpu = cuda.mem_alloc(a.nbytes) 
cuda.memcpy_htod(a_gpu, a) 

b = np.float32(2) 

mod = SourceModule(""" 
    __global__ void MatMult(float *a, float b) 
    { 
    const int i = threadIdx.x + blockDim.x * blockIdx.x; 
    const int j = threadIdx.y + blockDim.y * blockIdx.y; 
    int Idx = i + j*gridDim.x; 
    a[Idx] *= b; 
    } 
    """) 

func = mod.get_function("MatMult") 

xBlock = np.int32(np.floor(1024/rows)) 
yBlock = np.int32(cols) 

bdim = (xBlock, yBlock, 1) 
dx, mx = divmod(rows, bdim[0]) 
dy, my = divmod(cols, bdim[1]) 
gdim = ((dx + (mx>0)) * bdim[0], (dy + (my>0)) * bdim[1]) 
print "bdim=",bdim, ", gdim=", gdim 

func(a_gpu, b, block=bdim, grid=gdim) 
a_doubled = np.empty_like(a) 
cuda.memcpy_dtoh(a_doubled, a_gpu) 
print a_doubled - 2*a 

код должен напечатать размеры блоков bdim и размеры сетки gdim, а также как массив нулей.

Это работает для небольших размеров массива, например, если rows=256 и cols=10 (как в приведенном выше примере), выход следующим образом:

bdim= (4, 10, 1) , gdim= (256, 10) 
[[ 0. 0. 0. ..., 0. 0. 0.] 
[ 0. 0. 0. ..., 0. 0. 0.] 
[ 0. 0. 0. ..., 0. 0. 0.] 
..., 
[ 0. 0. 0. ..., 0. 0. 0.] 
[ 0. 0. 0. ..., 0. 0. 0.] 
[ 0. 0. 0. ..., 0. 0. 0.]] 

Однако, если я увеличить rows=512, я получаю следующий вывод :

bdim= (2, 10, 1) , gdim= (512, 10) 
[[ 0. 0. 0. ..., 0. 0. 0.] 
[ 0. 0. 0. ..., 0. 0. 0.] 
[ 0. 0. 0. ..., 0. 0. 0.] 
..., 
[ 2. 2. 2. ..., 2. 2. 2.] 
[ 2. 2. 2. ..., 2. 2. 2.] 
[ 2. 2. 2. ..., 2. 2. 2.]] 

Указывает, что умножение происходит дважды для некоторых элементов массива. не

Однако, если я заставляю размеры блоков для bdim = (1,1,1), проблема больше не возникает, и я получаю следующее (правильный) выход для большего размера массива:

bdim= (1, 1, 1) , gdim= (512, 10) 
[[ 0. 0. 0. ..., 0. 0. 0.] 
[ 0. 0. 0. ..., 0. 0. 0.] 
[ 0. 0. 0. ..., 0. 0. 0.] 
..., 
[ 0. 0. 0. ..., 0. 0. 0.] 
[ 0. 0. 0. ..., 0. 0. 0.] 
[ 0. 0. 0. ..., 0. 0. 0.]] 

Я не понимаю этого. Что здесь происходит, что означает, что этот метод определения размеров блока и сетки больше не подходит, поскольку размер массива увеличивается? Кроме того, если блок имеет размеры (1,1,1), значит ли это, что расчет выполняется серийно?

Заранее благодарим за любые указатели и помощь!

Answer 1

Работает на 2D-сетке 2D-блоков. В вашем ядре вы, кажется, считаете, что gridDim.x вернет число потоков в размере x размерности сетки.

__global__ void MatMult(float *a, float b) 
{ 
    const int i = threadIdx.x + blockDim.x * blockIdx.x; 
    const int j = threadIdx.y + blockDim.y * blockIdx.y; 
    int Idx = i + j*gridDim.x; 
    a[Idx] *= b; 
} 

gridDim.x возвращает число блоков г x направления сетки, а не количество потоков. Для того, чтобы получить число потоков в данном направлении следует умножать число потоков в блоке с числом блоков в сетке в том же направлении:

int Idx = i + j * blockDim.x * gridDim.x