Половинкин А.Н.

Постановка задачи Алгоритм умножения матриц на GPU Программная реализация

C = A*B A – прямоугольная матрица число строк – hA число столбцов – wA B – прямоугольная матрица число строк – hB число столбцов - wB

A sub B sub

Каждый блок потоков занимается вычислением одной подматрицы C sub матрицы С Каждый поток внутри блока потоков занимается вычислением одного элемента подматрицы C sub

C sub = A sub * B sub A sub размерности (block_size, wA) B sub размерности (wA, block_size) Матрицы A sub и B sub в общем случае могут не помещаться в общей памяти устройства, что приведёт к потере производительности Решение: разбить A sub и B sub на блоки A sub,i и B sub,i размерности (block_size, block_size), вычислять C sub как сумму произведений этих блоков:

for i = 1 to wA/block_size загрузить A sub,i и B sub,i из глобальной (device) памяти в общую (shared) память блока потоков. Каждый поток загружает «свой» элемент! каждый поток вычисляет «свой» элемент в произведении A sub,i *B sub,i и сохраняет аккумулированное значение end каждый поток загружает «свой» вычисленный элемент в глобальную (device) память

matrixMul.h – содержит определения (через define) размера блока и размеров матриц matrixMul_gold.cpp computeGold matrixMul.cu main randomInit printDiff runMultiplication matrixMul_kernel.cu matrixMul (kernel)

Используется схема хранения матриц по строкам в виде одномерного массива __global__ void matrixMul( float* C, float* A, float* B, int wA, int wB) вычисляем координаты текущего блока потоков и сохраняем их в переменные bx, by int bx = blockIdx.x int by = blockIdx.y вычисляем координаты текущего потока в блоке потоков и сохраняем их в переменные tx, ty

вычисляем индекс элемента в массиве, хранящем исходную матрицу A, который соответствует первому элементу первой обрабатываемой подматрицы A sub,1 int aBegin = wA * BLOCK_SIZE * by; вычисляем шаг для перехода к первому элементу следующей обрабатываемой подматрицы int aStep = BLOCK_SIZE; вычисляем индекс (условие останова перебора подматриц A sub,i ) int aEnd = aBegin + wA - 1;

вычисляем индекс элемента в массиве, хранящем исходную матрицу B, который соответствует первому элементу первой обрабатываемой подматрицы B sub,1 int bBegin = … вычисляем шаг для перехода к первому элементу следующей обрабатываемой подматрицы int bStep = … объявляем переменную, в которой хранится элемент произведения, вычисляемый текущим потоком float Csub = 0.f;

В цикле по всем подматрицам A sub,i, B sub,i выполняем следующие действия for (int a = aBegin, b = bBegin; a

каждый поток вычисляет свой элемент Csub в произведении подматриц A sub,i *B sub,i синхронизируем потоки в блоке потоков конец цикла загрузить вычисленный элемент Csub в соответствующий элемент матрицы С (в глобальную память) int c = wB * BLOCK_SIZE * by + BLOCK_SIZE * bx; C[c + wB * ty + tx] = Csub;

void runMultiplication(int argc, char** argv) инициализируем устройство (device) CUT_DEVICE_INIT(argc, argv); выделяем память на хосте для хранения матриц A и B unsigned int size_A = WA * HA; unsigned int mem_size_A = sizeof(float) * size_A; float* h_A = (float*)malloc(mem_size_A); unsigned int size_B = WB * HB; unsigned int mem_size_B = sizeof(float) * size_B; float* h_B = (float*)malloc(mem_size_B);

инициализируем матрицы A и B случайными значениями randomInit(h_A, size_A); randomInit(h_B, size_B); выделяем память под матрицы A и B на устройстве, копируем данные с хоста на устройство float* d_A; CUDA_SAFE_CALL(cudaMalloc((void**)&d_A, mem_size_A)); float* d_B; CUDA_SAFE_CALL(cudaMalloc((void**)&d_B, mem_size_B)); CUDA_SAFE_CALL(cudaMemcpy(d_A, h_A, mem_size_A, cudaMemcpyHostToDevice) ); CUDA_SAFE_CALL(cudaMemcpy(d_B, h_B, mem_size_B, cudaMemcpyHostToDevice) ); выделяем память под матрицу C на хосте и на устройстве (имена соответствующих переменных h_C, d_C, size_C, mem_size_C)

создаем и инициализируем таймер unsigned int timer = 0; CUT_SAFE_CALL(cutCreateTimer(&timer)); CUT_SAFE_CALL(cutStartTimer(timer)); определяем конфигурацию выполнения ядра (размер решетки блоков потоков и блока потоков) dim3 threads(BLOCK_SIZE, BLOCK_SIZE); dim3 grid(WC / threads.x, HC / threads.y); запускаем ядро копируем вычисленную матрицу С с устройства на хост

останавливаем таймер, выводим время вычислений, освобождаем ресурсы таймера CUT_SAFE_CALL(cutStopTimer(timer)); printf("Processing time: %f (ms) \n", cutGetTimerValue(timer)); CUT_SAFE_CALL(cutDeleteTimer(timer)); вычисляем то же самое произведение на CPU float* reference = (float*) malloc(mem_size_C); computeGold(reference, h_A, h_B, HA, WA, WB);

сравниваем результат, полученный на GPU, с результатом, полученным на CPU (по евклидовой норме) CUTBoolean res = cutCompareL2fe(reference, h_C, size_C, 1e-6f); printf("Test %s \n", (1 == res) ? "PASSED" : "FAILED"); if (res!=1) printDiff(reference, h_C, WC, HC); освобождаем память

Nvidia CUDA Programming Guide Многочисленные курсы по CUDA: (на русском языке)

?