Занятие 7. Программирование на CUDA C. Часть 3. - Лаборатория суперкомпьютерных и квантовых вычислений

Увеличение значение у элементов в матрице на единицу

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

Поэтому ниже будет приведен листинг всего кода с комментариями по новым вещам, что ранее в программах не было, а после каждый момент будет подробно разобран. Итак, листинг:

void incKernel ( float * data )
{
int idx = blockIdx.x * blockDim.x + threadIdx.x;
data [idx] = data [idx] + 1.0f;
}
 
int main ( int argc, char * argv [] )
{
int n = 16 * 1024 * 1024;
int numBytes = n * sizeof ( float );
 
// выделение памяти на хосте
float * a = new float [n];
for ( int i = 0; i < n; i++ )
a [i] = 0.0f;
 
// выделение памяти на девайсе
float * dev = NULL;
cudaMalloc ( (void**)&dev, numBytes );
 
// Установка конфигурации запуска ядра
 
dim3 threads = dim3(512, 1);
dim3 blocks = dim3(n / threads.x, 1);
 
// создание обработчиков событий cuda
 
cudaEvent_t start, stop;
float gpuTime = 0.0f;
 
cudaEventCreate ( &start );
cudaEventCreate ( &stop );
 
// асинхронно выдаем работу на GPU (все в поток 0)
 
cudaEventRecord ( start, 0 );
cudaMemcpy ( dev, a, numBytes, cudaMemcpyHostToDevice );
 
incKernel<<<blocks, threads>>>(dev);
 
cudaMemcpy ( a, dev, numBytes, cudaMemcpyDeviceToHost );
cudaEventRecord ( stop, 0 );
cudaEventSynchronize ( stop );
cudaEventElapsedTime ( &gpuTime, start, stop );
 
// Печатаем время работы на CPU и GPU
 
printf("time spent executing by the GPU: %.2f millseconds\n", gpuTime );
 
// проверка аутпута на корректность
 
printf("--------------------------------------------------------------\n");
 
for ( int i = 0; i < n; i++ )
if ( a [i] != 1.0f )
 
{
 
printf ( "Error in pos %d, %f\n", i, a [i] );
 
break;
 
}
 
// освобождение ресурсов
 
cudaEventDestroy ( start );
cudaEventDestroy ( stop );
cudaFree( dev );
 
delete a;
 
return 0;
 
}

В примере продемонстрирован обработчик событий, который не был ранее описан. Итак, обработчик событий CUDA cudaEvent_t start, stop задает начало и конец события. С помощью переменной float gpuTime = 0.0f будет производиться расчет времени выполнения программы в миллисекундах. Запись события производится с помощью специального объявления cudaEventRecord ( start/stop, 0 ). Вычисление времени производится с помощью cudaEventElapsedTime. Это время, которое потребовалось для вычисления только на стороне GPU. В последующих примерах работа с записью времени подробно описываться не будет, за исключением того, что не было описано здесь.

threadIdx и blockIdx определить какой именно элемент соответствует данному треду и увеличить именно его.

Так как и блоки и grid одномерны, то номер треда будет определятся как произведение номера блока на количество тредов в блоке, плюс номер треда внутри блока, т.е. blockIdx.x * blockDim.x + threadIdx.x.

cudaMalloc выделяется память под копию массива с данными в глобальной памяти (DRAM GPU). Далее данные копируются функцией cudaMemcpy из памяти CPU к глобальную память GPU.

По завершению копирования данных в глобальную память запускается ядро для обработки данных и после его вызова копируются результаты вычислений обратно из глобальной памяти GPU в память CPU.

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

Перемножение двух матриц

Следующий пример будет чуть сложнее, но не менее интереснее — перемножение двух квадратных матриц размерности N*N.

Допустим у нас есть две квадратные матрицы A и B размера N*N (будем считать, что N кратно 16). Простейший вариант использует по одному треду на каждый элемент получающейся матрицы C, при этом тред извлекает все необходимые элементы из глобальной памяти и производит требуемые вычисления.

Элемент c_i,jпроизведения двух матриц A и B вычисляется следующим фрагментом псевдокода:

Для каждого k-го элемента от 0 до N-1 считать:

c [i][j] += a [i][k] * b [k][j]

Таким образом для вычисления одного элемента произведения матриц нужно выполнить 2*N арифметических операций и 2*N чтений из глобальной памяти. Очевидно, что в основным лимитирующим фактором является скорость доступа к глобальной памяти, которая весьма низка. Опустим этот момент (его мы разберем в следующем примере), так как пример демонстрирует перемножение двух квадратных матриц, а не оптимизирует скорость доступа. Итак, листинг:

#define BLOCK_SIZE 16 // submatrix size
#define N 1024 // matrix size is N*N
 
__global__ void matMult ( float * a, float * b, int n, float * c )
 
{
int bx = blockIdx.x; // block index
int by = blockIdx.y;
int tx = threadIdx.x; // thread index
int ty = threadIdx.y;
 
float sum = 0.0f; // computed subelement
int ia = n * BLOCK_SIZE * by + n * ty; // a [i][0]
int ib = BLOCK_SIZE * bx + tx;
 
// Multiply the two matrices together;
 
for ( int k = 0; k < n; k++ )
sum += a [ia + k] * b [ib + k*n];
 
// Write the block sub-matrix to global memory;
 
// each thread writes one element
 
int ic = n * BLOCK_SIZE * by + BLOCK_SIZE * bx;
c [ic + n * ty + tx] = sum;
}
 
int main ( int argc, char * argv [] )
{
int numBytes = N * N * sizeof ( float );
 
// выделение памяти на хосте
 
float * a = new float [N*N];
float * b = new float [N*N];
float * c = new float [N*N];
 
for ( int i = 0; i < N; i++ )
for ( int j = 0; j < N; j++ )
{
int k = N*i + j;
a [k] = 0.0f;
b [k] = 1.0f;
}
 
// выделение памяти на девайсе
 
float * adev = NULL;
float * bdev = NULL;
float * cdev = NULL;
 
cudaMalloc ( (void**)&adev, numBytes );
cudaMalloc ( (void**)&bdev, numBytes );
cudaMalloc ( (void**)&cdev, numBytes );
 
// Установка конфигурации запуска ядра
 
dim3 threads ( BLOCK_SIZE, BLOCK_SIZE );
dim3 blocks ( N / threads.x, N / threads.y);
 
// Создание обработчика событий CUDA
 
cudaEvent_t start, stop;
float gpuTime = 0.0f;
cudaEventCreate ( &start );
cudaEventCreate ( &stop );
 
// асинхронно выдаваем работу на GPU (все в поток 0)
 
cudaEventRecord ( start, 0 );
cudaMemcpy ( adev, a, numBytes, cudaMemcpyHostToDevice );
cudaMemcpy ( bdev, b, numBytes, cudaMemcpyHostToDevice );
 
matMult<<<blocks, threads>>> ( adev, bdev, N, cdev );
 
cudaMemcpy ( c, cdev, numBytes, cudaMemcpyDeviceToHost );
cudaEventRecord ( stop, 0 );
cudaEventSynchronize ( stop );
cudaEventElapsedTime ( &gpuTime, start, stop );
 
// Печатаем время работы на GPU и CPU
 
printf("time spent executing by the GPU: %.2f millseconds\n", gpuTime );
 
// Освобождение ресурсов
 
cudaEventDestroy ( start );
cudaEventDestroy ( stop );
cudaFree ( adev );
cudaFree ( bdev );
cudaFree ( cdev );
 
delete a;
delete b;
delete c;
 
return 0;
 
}

Естественно, саму матрицу вы можете задать произвольно, в данном примере была указана лишь размерность без каких-либо значений элементов. Данный пример демонстрирует производительность перемножения двух матриц размерностей N*N:

Перемножение двух матриц с использованием shared-памяти

Поскольку используется разбиение бэндов на квадратные подматрицы, то псевдокод расчета элемент c_i,j изменится:

Для каждого шага от 0 до N/16 считать:

По каждому k-му элементу от 0..16 считать:

c [i][j] += a [i][k+step*16] * b [k+step*16][j]

Следует обратить внимание, что для каждого значения step значения из матриц A и B берутся из двух подматриц размером 16*16. Фактически полосы с рис. 4 просто поделены на квадратные подматрицы и каждому значению step соответствует по одной такой подматрице A и одной подматрице B.

На каждом шаге будем загружать в shared-память по одной 16*16 подматрице A и одной 16*16 подматрице B. Далее будем вычислять соответствующую им сумму для элементов произведения, потом загружаем следующие 16*16-подматрицы и т.д.

При этом на каждом шаге один тред загружает ровно по одному элементу из каждой из матриц A и B и вычисляет соответствующую им сумму членов. По окончании всех вычислений производится запись элемента в итоговую матрицу.

Следует обратить внимание, что после загрузки элементов из A и B нужно выполнить синхронизацию тредов при помощи вызова __synchronize для того, чтобы к моменту начала расчетов все нужные элементы (загружаемые остальными тредами блока) были бы уже загружены. Точно также по окончании обработки загруженных подматриц и перед загрузкой следующих также нужна синхронизация.

Ниже приводится соответствующий исходный код:

#define BLOCK_SIZE 16 // submatrix size
#define N 1024 // matrix size is N*N
 
__global__ void matMult ( float * a, float * b, int n, float * c )
{
int bx = blockIdx.x; // block index
int by = blockIdx.y;
int tx = threadIdx.x; // thread index
int ty = threadIdx.y;
 
// Index of the first sub-matrix of A processed by the block
int aBegin = n * BLOCK_SIZE * by;
int aEnd = aBegin + n - 1;
 
// Step size used to iterate through the sub-matrices of A
int aStep = BLOCK_SIZE;
 
// Index of the first sub-matrix of B processed by the block
int bBegin = BLOCK_SIZE * bx;
 
// Step size used to iterate through the sub-matrices of B
int bStep = BLOCK_SIZE * n;
float sum = 0.0f; // computed subelement
for ( int ia = aBegin, ib = bBegin; ia <= aEnd; ia += aStep, ib += bStep )
{
 
// Shared memory for the sub-matrix of A
__shared__ float as [BLOCK_SIZE][BLOCK_SIZE];
 
// Shared memory for the sub-matrix of B
__shared__ float bs [BLOCK_SIZE][BLOCK_SIZE];
 
// Load the matrices from global memory to shared memory;
as [ty][tx] = a [ia + n * ty + tx];
bs [ty][tx] = b [ib + n * ty + tx];
 
__syncthreads(); // Synchronize to make sure the matrices are loaded
 
// Multiply the two matrices together;
for ( int k = 0; k < BLOCK_SIZE; k++ )
sum += as [ty][k] * bs [k][tx];
 
// Synchronize to make sure that the preceding
// computation is done before loading two new
// sub-matrices of A and B in the next iteration
 
__syncthreads();
}
 
// Write the block sub-matrix to global memory;
// each thread writes one element
 
int ic = n * BLOCK_SIZE * by + BLOCK_SIZE * bx;
c [ic + n * ty + tx] = sum;
 
}
int main ( int argc, char * argv [] )
{
int numBytes = N * N * sizeof ( float );
 
// allocate host memory
 
float * a = new float [N*N];
float * b = new float [N*N];
float * c = new float [N*N];
 
for ( int i = 0; i < N; i++ )
for ( int j = 0; j < N; j++ )
{
a [i] = 0.0f;
b [i] = 1.0f;
}
 
// allocate device memory
 
float * adev = NULL;
float * bdev = NULL;
float * cdev = NULL;
 
cudaMalloc ( (void**)&adev, numBytes );
cudaMalloc ( (void**)&bdev, numBytes );
cudaMalloc ( (void**)&cdev, numBytes );
 
// set kernel launch configuration
 
dim3 threads ( BLOCK_SIZE, BLOCK_SIZE );
dim3 blocks ( N / threads.x, N / threads.y);
 
// create cuda event handles
 
cudaEvent_t start, stop;
float gpuTime = 0.0f;
cudaEventCreate ( &start );
cudaEventCreate ( &stop );
 
// asynchronously issue work to the GPU (all to stream 0)
 
cudaEventRecord ( start, 0 );
cudaMemcpy ( adev, a, numBytes, cudaMemcpyHostToDevice );
cudaMemcpy ( bdev, b, numBytes, cudaMemcpyHostToDevice );
 
matMult<<<blocks, threads>>> ( adev, bdev, N, cdev );
 
cudaMemcpy ( c, cdev, numBytes, cudaMemcpyDeviceToHost );
cudaEventRecord ( stop, 0 );
cudaEventSynchronize ( stop );
cudaEventElapsedTime ( &gpuTime, start, stop );
 
// print the cpu and gpu times
 
printf("time spent executing by the GPU: %.2f millseconds\n", gpuTime );
 
// release resources
 
cudaEventDestroy ( start );
cudaEventDestroy ( stop );
 
cudaFree ( adev );
cudaFree ( bdev );
cudaFree ( cdev );
 
delete a;
delete b;
delete c;
 
return 0;
 
}

Теперь для вычисления одного элемента произведения матриц нам нужно всего 2*N/16 чтений из глобальной памяти. И по результатам сразу видно за счет использования shared-памяти нам удалось поднять быстродействие более чем на порядок – 3.66 милисекунд: