ТЕХНОЛОГИЯ ПРОГРАММИРОВАНИЯ CUD. Учебное пособие казань 2017 2 удк 004. 4, 681. 3 Ббк 32. 81 Печатается по постановлению Редакционноиздательского совета
 Скачать 1.28 Mb.
|
|
Глава 4. Атомарные операции в CUDA CUDA для GPU, начиная с производительности 1.1 и выше, поддерживает атомарные операции над глобальной и разделяемой памятью. Атомарность заключается в том, что гарантируется корректность выполнения операции в случае, когда много нитей пытаются ее выполнить одновременно. Атомарные арифметические операции Наиболее часто используемыми атомарными арифметическими операциями являются atomicAdd и atomicSub, служащие для увеличения или уменьшения величины на заданное значение. При этом функция возвращает старое значение. Нужно обратить внимание, что atomicAdd поддерживает операции над 64-битовыми величинами, но только расположенными в глобальной памяти. int atomicAdd( int *addr, int value ); unsigned int atomicAdd( unsigned int *addr, unsigned int value ); unsigned long long atomicAdd(unsigned long long *addr, unsigned long long value ); int atomicSub( int *addr, int value ); unsigned int atomicSub(unsigned int *addr, unsigned int value ); Операция atomicExch осуществляет атомарный обмен значениями – передаваемое значение записывается по указанному адресу, а предыдущее значение возвращается. При этом подобный обмен происходит как одна транзакция, то есть ни одна нить не может «вклиниться» между шагами этого 35 обмена. Операции над 64-битовыми значениями поддерживаются только для глобальной памяти. int atomicExch( int *addr, int value ); unsigned int atomicExch(unsigned int *addr, unsigned int value ); unsigned long long atomicExch(unsigned long long *addr, unsigned long long value ); float atomicExch( float *addr, float value ); Следующие две операции сравнивают значение по адресу с переданным значением, записывают минимум/максимум из этих двух значений по заданному адресу и возвращают предыдущее значение, находившееся по адресу. Все эти шаги выполняются атомарно, как одна транзакция. int atomicMin( int *addr, int value ); unsigned int atomicMin(unsigned int *addr, unsigned int value ); int atomicMax( int *addr, int value ); unsigned int atomicMax(unsigned int *addr, unsigned int value ); Операция atomicInc читает слово по заданному адресу и сравнивает его с переданным значением. Если прочитанное слово больше, то по адресу записывается нуль, иначе значение по адресу увеличивается на единицу. Возвращается старое значение. 36 unsigned int atomicInc(unsigned int *addr, unsigned int value ); Операция atomicDec читает слово по переданному адресу. Если прочитанное значение равно нулю или больше переданного значения, то записывает по адресу переданное значение, иначе уменьшает значение по адресу на единицу. Возвращается старое значение. unsigned int atomicDec(unsigned int *addr, unsigned int value ); Следующая функция (CAS – Compare And Swap) читает старое 32- или 64-битовое значение по переданному адресу и сравнивает его с параметром compare. В случае совпадения по переданному адресу записывается значение параметра value, иначе значение по адресу не изменяется. Возвращается всегда старое прочитанное значение. int atomicCAS( int *addr, int compare, int value ); unsigned int atomicCAS(unsigned int *addr, unsigned int compare, unsigned int value ); unsigned long long atomicCAS(unsigned int *addr, unsigned long long compare, unsigned long long value ); Атомарные побитовые операции Побитовые атомарные операции читают слово по заданному адресу, применяют к нему заданную побитовую операцию с заданным параметром и записывают результат обратно. Возвращается всегда старое значение, находившееся по переданному адресу до начала операции. int atomicAnd( int *addr, int value ); 37 unsigned int atomicAnd(unsigned int *addr, unsigned int value ); int atomicOr( int *addr, int value ); unsigned int atomicOr(unsigned int *addr, unsigned int value ); int atomicXor( int *addr, int value ); unsigned int atomicXor(unsigned int *addr, unsigned int value ); Дополнительные возможности Compute Capability 6.x Начиная с производительности 6.x, существует ряд новвоведений. В частности, представлен новый тип атомарных операций. Например, atomicAdd_system гарантирует, что инструкция является атомарной относительно других CPU и GPU в системе. Функция atomicAdd_block подразумевает, что инструкция является атомарной только с учетом атомистики от других потоков в одном блоке потока. В следующем примере как CPU, так и GPU могут атомарно обновлять целочисленное значение по адресу addr: __global__ void add10(int *addr) { // only available on devices with compute // capability 6.x atomicAdd_system(addr, 10); } void foo() { int *addr; cudaMallocManaged(&addr, 4); *addr = 0; 38 add10<<<2, 5>>>(addr); // CPU atomic operation _sync_fetch_and_add(addr, 10); } Также в версиях производительности 6.x добавлено сложение 64 битных чисел. Пример использования атомарных операций Пусть задана двумерная сетка нитей, и требуется сосчитать количество нитей, индексы которых удовлетворяют некоторому условию. Эту задачу можно реализовать следующим способом: каждая нить вызывает некую булеву функцию, и если результат вызова этой функции равна true, добавляет единицу некой глобальной переменной res. Приведем ядро такой программы. __global__ void Kernel(int * res) { int i = blockIdx.x * blockDim.x + threadIdx.x; int j = blockIdx.y * blockDim.y + threadIdx.y; if ( Сondition(i, j) )// проверка условия { atomicAdd(res, 1); } } Здесь Сondition(i, j)некоторое условие, например, i>=j. Вопросы к разделу 1. Что такое атомарная операция? В чем ее ключевое отличие от обычных операций? 2. Какие атомарные операции вы знаете? 39 Лабораторная работа 1. Вычислите приблизительное значение числа π методом Монте Карло. Постройте равномерную сетку нитей в квадрате [0, 1]X[0, 1] и посчитайте количество нитей, которые попали в круг с центром в нуле и с радиусом, равным единице. Для вычисления случайного числа напишите __device__ функцию. Домашняя работа 1. Вызов атомарной функции, в целом, замедляет работу программы. Например, если очень много нитей должны добавить единицу одной и той же переменной, то они встанут в очередь, и будут простаивать, пока все по очереди не выполнят эту операцию. Одним из подходов для решения этой проблемы является уменьшение количества нитей, вызывающих атомарную операцию. Учитывая выше сказанное, подумайте, как можно уменьшить количество нитей, вызывающих атомарную операцию в методе Монте Карло. Оптимизируйте таким образом программу, и посмотрите, насколько уменьшается время ее работы. Также посмотрите, как меняется точность вычисления числа π. 40 Глава 5. Работа с векторами. Математические функции В языке Си (в частности, и в диалекте CUDA C) векторы обычно представляются массивами, как правило, динамическими (память выделяется с помощью функций семейства malloc() и освобождается вызовом free(), аналоги в CUDA см. ниже). В функции же передается не сам массив целиком, а лишь указатель на его начало (переменная фактически целочисленного типа, размер которого зависит от архитектуры (как правило, 32 или 64 разряда), содержащая адрес начала массива). Зная размер одного элемента, легко вычислить адрес любого элемента по его смещению. Сложение векторов Напишем программу, которая генерирует два вектора целых чисел, копирует их в память GPU, выполняет сложение на GPU, результат копируем обратно в память CPU. GPU обычно не работает на прямую с памятью CPU, хотя это технически возможно. Для начала необходимо выделить память под данные на device, для этого используется фукция cudaMalloc(), после этого копируем данные из памяти host в память device с помощью cudaMemcpy(). Далее для каждого элемента вектора запускается отдельный процесс сложения, т.е. в данном случае – 128 нитей. На CPU это могло быть накладно, но GPU как раз рассчитан на работу с большим количеством потоков, и всё происходит достаточно быстро. После того, как операция завершена с помощью cudaMemcpy(), копируем результат обратно на host и освобождаем память device с помощью cudaFree(). В тексте ядра есть переменная threadIdx, это системная переменная CUDA, содержащая номер нити. Приведем программу сложения двух векторов с помощью технологии CUDA. #define N 128 //количество элементов массива // Ядро 41 __global__ void add( int *a, int *b, int *c ) { //вычисление индекса элемента int tid = threadIdx.x; //проверка на выход за пределы массива if(tid > N-1) return; //поэлементное сложение массивов c[tid] = a[tid] + b[tid]; } int main() { // выделение памяти под массивы на CPU int host_a[N], host_b[N], host_c[N]; // выделение памяти под массивы для копирования // на GPU int *dev_a, *dev_b, *dev_c; // заполнение массивов for (int i=0; i } // выделение памяти под массивы на GPU cudaMalloc( (void**)&dev_a, N * sizeof(int) ); cudaMalloc( (void**)&dev_b, N * sizeof(int) ); cudaMalloc( (void**)&dev_c, N * sizeof(int) ); // копирование данных в память GPU cudaMemcpy( dev_a, host_a, N * sizeof(int), cudaMemcpyHostToDevice ) ; cudaMemcpy( dev_b, host_b, N * sizeof(int), cudaMemcpyHostToDevice ) ; // вызов ядра add<<<1, N>>>( dev_a, dev_b, dev_c ); // получаем результат расчета 42 cudaMemcpy( host_c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost ) ; // вывод результатов for (int i=0; i } // освобождение памяти cudaFree( dev_a ); cudaFree( dev_b ); cudaFree( dev_c ); return 0; } Вычисление математических функций В CUDA существуют аналоги ряда математических функций, включенных в стандартную библиотеку языка Си ( __device__ float __sinf(float x); Перечислим основные различия: - функции CUDA вычисляются непосредственно на GPU; - большинство функций работают только с float-аргументами (GPU, как правило, в основном работают с float, точности достаточно для основных задач, возлагаемых на GPU); - имеют модификатор __device__, вследствие чего могут вызываться только с GPU (т.е. из kernel-функций); Рассмотрим следующую задачу: 1. Выделить на GPU массив из n = 10 9 элементов типа float. 2. Инициализировать их в ядре следующим образом: 43 arr[i] = sin((i % 360) * pi / 180) где pi = 3.141592653. 3. Скопировать массив в память CPU, посчитать ошибку err = sum_i(abs(sin((i % 360) * pi / 180) - arr[i])) / n 4. Провести исследование при использовании массива типа double. Основной частью данной задачи является написание ядра, в котором будет инициализироваться массив (и вызываться синус). #define BASE_TYPE float #define MAX_GRIDSIZE 1 __global__ void sinMass(BASE_TYPE *A, int arraySize) { int index = (blockIdx.y * MAX_GRIDSIZE + blockIdx.x) * blockDim.x + threadIdx.x; if (index < arraySize) [index] = sin((BASE_TYPE)((index % 360) * M_PI / 180)); // другие варианты вычисления синуса // A[index] = sinf((index % 360) * M_PI / 180); // A[index] = __sinf((index % 360) * M_PI // / 180); } Приведем список встроенных математических функций для работы с типом float: 44 Функция Операция __fadd_[rn,rz,ru,rd](x,y), x + y __fsub_[rn,rz,ru,rd](x,y), x – y __fmul_[rn,rz,ru,rd](x,y), x * y __fmaf_[rn,rz,ru,rd](x,y,z) x * y + z __frcp_[rn,rz,ru,rd](x) 1 / x __fsqrt_[rn,rz,ru,rd](x) √x __frsqrt_rn(x) 1 / √x __fdiv_[rn,rz,ru,rd](x,y) x / y __fdividef(x,y) x / y; при y>2 126 результат = 0 __expf(x) e x __exp10f(x) 10 x __logf(x) ln x __log2f(x) log 2 x __log10f(x) lg x __sinf(x) sin x __cosf(x) cos x __sincosf(x,sptr,cptr) *sptr = sin x, cptr = cos x __tanf(x) tg x __powf(x, y) x y Здесь rn,rz,ru,rd - способы округления, например, __fadd_rn округляет сумму к ближайшему, __fadd_rz – к нулю, __fadd_ru – до большего, __fadd_rd – до меньшего. Список встроенных функций для работы с типом double несколько меньше: 45 Функция Операция __dadd_[rn,rz,ru,rd](x,y), x + y __dsub_[rn,rz,ru,rd](x,y), x – y __dmul_[rn,rz,ru,rd](x,y), x * y __fma_[rn,rz,ru,rd](x,y,z) x * y + z __drcp_[rn,rz,ru,rd](x) 1 / x __dsqrt_[rn,rz,ru,rd](x) √x __ddiv_[rn,rz,ru,rd](x,y) x / y Вопросы к разделу 1. Как организовать обращение к элементу массива в теле ядра? 2. Какие основные этапы работы с векторами можно выделить? 3. Какие встроенные (быстрые) математические функции есть на GPU? В чем отличие их от обычных функций? Лабораторная работа 1. Используя другие функции вместо синуса (__expf` (e^x), __exp10f` (10^x) и т.д.), рассчитайте относительную ошибку (отношение абсолютной ошибки к значению функции). 2. Напишите программу, реализующую скалярное произведение двух векторов. Домашняя работа 1. Сравните скорость и точность вычислений с типами float и double для функций из лабораторной работы. Какой тип для каких задач лучше использовать и почему? Сделайте вывод. 2. Напишите программу, реализующую процесс ортогонализации Грама- Шмидта. При этом на хосте создается набор из N векторов вида a 1 = (1,1,1,1,…), a 2 = (0,1,1,1,…), a 3 = (0,0,1,1,…), … , a N = (…,0,0,0,1). Далее эти вектора передаются на девайс, после чего производится 46 ортогонализация. Алгоритм процесса ортогонализации (получение ортогональных векторов b i по a i ) следующий: 2 , ) , ( ) , ( , 1 1 1 1 N i b b b b a a b a b k i k k k k i i i = − = = ∑ − = 47 Глава 6. Работа с матрицами Перед тем, как начать работать с матрицами, скажем несколько слов о типе dim3 – одном из базовых типов CUDA. Он представляет собой не что иное, как вектор из 3-х целочисленных значений. Создается он следующим образом: dim3 v1(512, 2, 4); dim3 v2(512, 128); dim3 v3(1024). Если какие-то значения опущены (как в примере с v2 и v3), то они по умолчанию равны 1 (а не 0). То есть, вышеприведенный код эквивалентен следующему: dim3 v1(512, 2, 4); dim3 v2(512, 128, 1); dim3 v3(1024, 1, 1). В CUDA есть несколько векторных типов для доступа к данным. Например, такой тип, как [u]int может быть 1-, 2-, 3-, 4-мерным вектором. Для создания переменной такого типа применяется функция специального вида, например: int3 a = make_int3 ( 1, 2, 3); Так же существуют типы (u)char, (u)short, (u)long, float, которые могут быть 1-, 2-, 3-, 4-мерными векторами, а вот типы long long и double могут быть только 1- и 2-мерными. На основе типа uint3, в CUDA существует специальный тип dim3, который имеет конструктор. Конструктор этого типа позволяет задавать не все компоненты 3-мерного вектора. Недостающие компоненты этого типа в конструкторе инициализируются единицами. Главное назначение типа dim3 – это задание параметров запуска ядра. Внутри каждого вычислительного ядра содержатся встроенные переменные, которые используют вышеописанные типы. С помощью этих 48 встроенных переменных можно оперировать отдельными нитями и различать их. Эти переменные были описаны выше, кроме переменной warpSize типа int – эта переменная содержит значение размера связки нитей (варп), которое на текущий момент всегда равно 32. Используя встроенные переменные, можно задавать разные параметры запуска ядра и различать отдельные нити в блоках или сами блоки. Создание матриц Создадим на CPU матрицу A размером 10x10 и зададим элементы матрицы следующим образом: = 99 92 91 90 29 22 21 20 19 12 11 10 9 2 1 0 A Создадим на GPU матрицу B размером 10x10. Элементы матрицы B зададим такими же, как и у матрицы А. При этом не будем использовать копирование данных с хоста на девайс. Необходимо скопировать матрицу B с девайса на хост и сравнить на хосте две матрицы – A и B. Нужно добиться, чтобы обе матрицы совпадали. Попробуйте запускать программу с разными размерами grid и block. Для решения этой задачи используем один блок с десятью нитями так, чтобы каждая нить создавала один элемент матрицы. Для удобства, блок нитей мы сделаем двумерным. Таким образом, индекс элемента матрицы вычислим как threadIdx.y * n + threadIdx.x. В итоге код ядра будет выглядеть так: __global__ void createMatrix(int *A, const int n) { 49 A[threadIdx.y * n + threadIdx.x] = 10 * threadIdx.y + threadIdx.x; } На CPU матрицу А будем создавать с помощью двойного цикла следующим образом: size_t size = n * n * sizeof(int); int *h_A = (int *)malloc(size); for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) h_A[j * n + i] = 10 * j + i; Вызов ядра, описанного выше, осуществляется таким образом: dim3 threadsPerBlock = dim3(10, 10); dim3 blocksPerGrid = dim3(1); createMatrix<< Проверку равенства двух матриц осуществляем следующим образом с помощью двойного цикла: for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) if (h_A[j * n + i] != h_B[j * n + i]) { printf("h_A[%d] != h_B[%d]\n", j * n + i, j * n + i); } Листинг кода. // Сравнение двух матриц #include "cuda_runtime.h" 50 #include "device_launch_parameters.h" #include #include #include // Ядро __global__ void createMatrix(int *A, const int n) { // Создание элементов матрицы на GPU A[threadIdx.y * n + threadIdx.x] = 10 * threadIdx.y + threadIdx.x; } int main() { // кол-во строк и столбцов матрицы const int n = 10; // размер матрицы size_t size = n * n * sizeof(int); // выделяем память для матрицы на CPU int *h_A = (int *)malloc(size); // инициализируем матрицу for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) h_A[j * n + i] = 10 * j + i; int *d_B = NULL; // выделяем память для матрицы на GPU cudaMalloc((void **)&d_B, size); // определение размеров сетки и блоков dim3 threadsPerBlock = dim3(10, 10); dim3 blocksPerGrid = dim3(1); // вызов ядра createMatrix<< 51 // выделяем память для матрицы B, чтобы // скопировать из GPU на CPU int *h_B = (int *)malloc(size); // копируем матрицу из GPU на CPU cudaMemcpy(h_B, d_B, size, cudaMemcpyDeviceToHost); // проверяем совпадение матрицы А и матрицы В for (int i = 0; i < n; i++) for (int j = 0; j < n; j++) if (h_A[j * n + i] != h_B[j * n + i]) { printf("h_A[%d] != h_B[%d]\n", j * n + i, j * n + i); } // освобождаем память на GPU cudaFree(d_B); // освобождаем память на CPU free(h_A); free(h_B); return 0; } Транспонирование матрицы Для решения этой задачи используем блоки размером 16х16. Каждая нить блока создает один элемент транспонированной матрицы. Код ядра будет выглядеть так: __global__ void matrixTranspose (const BASE_TYPE *A, BASE_TYPE *AT, int rows, int cols) { // Индекс элемента в исходной матрице 52 int iA = cols * (blockDim.y * blockIdx.y + threadIdx.y) + blockDim.x * blockIdx.x + threadIdx.x; // Индекс соответствующего элемента в // транспонированной матрице int iAT = rows * (blockDim.x * blockIdx.x + threadIdx.x) + blockDim.y * blockIdx.y + threadIdx.y; AT[iAT] = A[iA]; } Для удобства вычислений размеры исходной матрицы делаем кратными размеру блока. Для этого используем функцию toMultiple(int a, int b), которая возвращает число, кратное b и наиболее близкое к a. Матрицу инициализируем случайными числами. // Инициализация матрицы for (int i = 0; i < rows * cols; ++i) { h_A[i] = rand()/(BASE_TYPE)RAND_MAX; } Вызов ядра, описанного выше, осуществляется таким образом: dim3 threadsPerBlock = dim3(BLOCK_SIZE, BLOCK_SIZE); dim3 blocksPerGrid = dim3(cols / BLOCK_SIZE, rows / BLOCK_SIZE); //вызов ядра matrixTranspose<< Осуществляем проверку правильности работы ядра. 53 for (int i = 0; i < rows; i++) for (int j = 0; j < cols; j++) if (h_A[i * cols + j] != h_AT[j * rows + i]) { fprintf(stderr, "Result verification failed at element [%d, %d]!\n", i, j); exit(EXIT_FAILURE); } Листинг кода: #include "cuda_runtime.h" #include "device_launch_parameters.h" #include #include // размер блока #define BLOCK_SIZE 16 // тип, который будут иметь элементы матриц #define BASE_TYPE float // Ядро // Функция транспонирования матрицы __global__ void matrixTranspose(const BASE_TYPE *A, BASE_TYPE *AT, int rows, int cols) { // Индекс элемента в исходной матрице int iA = cols * (blockDim.y * blockIdx.y + threadIdx.y) + blockDim.x * blockIdx.x + threadIdx.x; // Индекс соответствующего элемента в // транспонированной матрице int iAT = rows * (blockDim.x * blockIdx.x + threadIdx.x) + blockDim.y * blockIdx.y + threadIdx.y; 54 AT[iAT] = A[iA]; } // Функция вычисления числа, которое больше числа а // и кратное числу b int toMultiple(int a, int b) { int mod = a % b; if (mod != 0) { mod = b - mod; return a + mod; } return a; } int main() { // Объекты событий cudaEvent_t start, stop; cudaEventCreate(&start); cudaEventCreate(&stop); // Количество строк и столбцов матрицы int rows = 1000; int cols = 2000; // Меняем количество строк и столбцов матрицы // на число, кратное размеру блока (16) rows = toMultiple(rows, BLOCK_SIZE); printf("rows = %d\n", rows); cols = toMultiple(cols, BLOCK_SIZE); printf("cols = %d\n", cols); size_t size = rows * cols * sizeof(BASE_TYPE); // Выделение памяти под матрицы на хосте // Исходная матрица BASE_TYPE *h_A = (BASE_TYPE *)malloc(size); 55 // Транспонированная матрица BASE_TYPE *h_AT = (BASE_TYPE *)malloc(size); // Инициализация матрицы for (int i = 0; i < rows * cols; ++i) { h_A[i] = rand()/(BASE_TYPE)RAND_MAX; } // Выделение глобальной памяти на девайсе // для исходной матрицы BASE_TYPE *d_A = NULL; cudaMalloc((void **)&d_A, size); // Выделение глобальной памяти на девайсе для // транспонированной матрицы BASE_TYPE *d_AT = NULL; cudaMalloc((void **)&d_AT, size); // Копируем матрицу из CPU на GPU cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice); // Определяем размер блока и сетки dim3 threadsPerBlock = dim3(BLOCK_SIZE, BLOCK_SIZE); dim3 blocksPerGrid = dim3(cols / BLOCK_SIZE, rows / BLOCK_SIZE); // Начать отсчета времени cudaEventRecord( start, 0); // Запуск ядра matrixTranspose<< // Окончание работы ядра, остановка времени cudaEventRecord( stop, 0); cudaEventSynchronize( stop ); 56 float KernelTime; cudaEventElapsedTime( &KernelTime, start, stop); printf("KernelTime: %.2f milliseconds\n", KernelTime); // Копируем матрицу из GPU на CPU cudaMemcpy(h_AT, d_AT, size, cudaMemcpyDeviceToHost); // Проверка правильности работы ядра for (int i = 0; i < rows; i++) for (int j = 0; j < cols; j++) { if (h_A[i * cols + j] != h_AT[j * rows + i]) fprintf(stderr, "Result verification failed at element [%d, %d]!\n", i, j); exit(EXIT_FAILURE); } printf("Test PASSED\n"); // Освобождаем память на GPU cudaFree(d_A); // Освобождаем память на GPU cudaFree(d_AT); // Освобождаем память на CPU free(h_A); free(h_AT); // Удаляем объекты событий cudaEventDestroy( start ); cudaEventDestroy( stop ); printf("Done\n"); return 0; } 57 Сложение матриц Поскольку сложение матриц выполняется поэлементно (еще проще, чем умножение), то все сказанное в следующем разделе (по поводу параллельного вычисления и представления массивов в памяти) справедливо и для сложения матриц. Для решения задачи используем функцию __global__ void matrixAdd. // Ядро __global__ void matrixAdd(const BASE_TYPE *A, const BASE_TYPE *B, BASE_TYPE *C, int cols) { // Вычисление индекса элемента матрицы на GPU int ind = cols * (blockDim.y * blockIdx.y + threadIdx.y) + blockDim.x * blockIdx.x + threadIdx.x; C[ind] = A[ind] + B[ind]; } // Функция вычисление числа, которое больше // числа а и кратное числу b int toMultiple(int a, int b) { int mod = a % b; if (mod != 0) { mod = b - mod; return a + mod; } return a; } 58 Умножение матриц Пусть даны матрицы 𝐴𝐴 = (𝑎𝑎 𝑖𝑖𝑖𝑖 ) и 𝐵𝐵 = (𝑏𝑏 𝑖𝑖𝑖𝑖 ) размерами 𝑚𝑚 × 𝑛𝑛 и 𝑛𝑛 × 𝑝𝑝 соответственно. Тогда произведение 𝐶𝐶 = (𝑐𝑐 𝑖𝑖𝑖𝑖 ) = 𝐴𝐴𝐵𝐵 размерностью 𝑚𝑚 × 𝑝𝑝, как известно, определяется по формуле 𝑐𝑐 𝑖𝑖𝑖𝑖 = ∑ 𝑎𝑎 𝑖𝑖𝑖𝑖 𝑏𝑏 𝑖𝑖𝑖𝑖 𝑛𝑛 𝑖𝑖=1 Элементы матрицы-произведения зависят только от элементов исходных матриц , поэтому их вычисление может выполняться параллельно, не требуя дополнительных "ухищрений". Поэтому рассмотрим лишь вопрос о вариантах представления матрицы в памяти – в виде двухмерного (например, float** a, затем обращение вида a[i][j]) или одномерного массива (float* a и затем обращение вида a[i * N + j], где N - длина строки (количество столбцов)). Рассмотрим вопрос о том, как работает операция индексирования: 1) в случае с одномерным массивом: a) вычисляется выражение i * N + j (в современных CPU это выполняется за одну операцию, так называемую fused multiply-add); b) полученное число (смещение) умножается на размер типа и складывается с базовым адресом (адресом начала массива); c) как правило, переменные i, j и базовый адрес будут храниться в регистрах CPU, поэтому вышеприведенные операции выполняются быстро. 2) в случае с двухмерным массивом: a) значение i умножается на размер указателя и складывается с базовым адресом, вычисляется адрес a[i] (начала i-й строки); b) Значение j умножается на размер элемента и складывается с адресом a[i], вычисляется адрес a[i][j]; c) Адрес a[i], очевидно, не хранится в регистре, поэтому операция с ним будет медленнее. С другой стороны, одномерный массив требует большого последовательного участка памяти. Для решения задачи опишем функцию __global__ void matrixMult. 59 __global__ void matrixMult(const BASE_TYPE *A, const BASE_TYPE *B, BASE_TYPE *C, int Acols, int Bcols) { int i0 = Acols * (blockDim.y * blockIdx.y + threadIdx.y); int j0 = blockDim.x * blockIdx.x + threadIdx.x; BASE_TYPE sum = 0; for (int k = 0; k < Acols; k++) sum += A[i0 + k] * B[k * Bcols + j0]; int ind = Bcols * (blockDim.y * blockIdx.y + threadIdx.y) + blockDim.x * blockIdx.x + threadIdx.x; C[ind] = sum; } Листинг кода: #include "cuda_runtime.h" #include "device_launch_parameters.h" #include #include #define BLOCK_SIZE 16 // тип, который будут иметь элементы матриц #define BASE_TYPE double // функция перемножения матриц __global__ void matrixMult(const BASE_TYPE *A, const BASE_TYPE *B, BASE_TYPE *C, int Acols, int Bcols) { int i0 = Acols * (blockDim.y * blockIdx.y + threadIdx.y); int j0 = blockDim.x * blockIdx.x + threadIdx.x; BASE_TYPE sum = 0; 60 for (int k = 0; k < Acols; k++) sum += A[i0 + k] * B[k * Bcols + j0]; int ind = Bcols * (blockDim.y * blockIdx.y + threadIdx.y) + blockDim.x * blockIdx.x + threadIdx.x; C[ind] = sum; } int toMultiple(int a, int b) { int mod = a % b; if (mod != 0) { mod = b - mod; return a + mod; } return a; } int main() { //start, stop - for Kernel time cudaEvent_t start, stop; cudaEventCreate(&start); cudaEventCreate(&stop); // количество строк и столбцов матрицы int Arows = 100; int Acols = 200; int Brows = Acols; int Bcols = 150; Arows = toMultiple(Arows, BLOCK_SIZE); printf("Arows = %d\n", Arows); Acols = toMultiple(Acols, BLOCK_SIZE); printf("Acols = %d\n", Acols); 61 Brows = toMultiple(Brows, BLOCK_SIZE); printf("Brows = %d\n", Brows); Bcols = toMultiple(Bcols, BLOCK_SIZE); printf("Bcols = %d\n", Bcols); size_t Asize = Arows * Acols * sizeof(BASE_TYPE); size_t Bsize = Brows * Bcols * sizeof(BASE_TYPE); size_t Csize = Arows * Bcols * sizeof(BASE_TYPE); BASE_TYPE *h_A = (BASE_TYPE *)malloc(Asize); BASE_TYPE *h_B = (BASE_TYPE *)malloc(Bsize); BASE_TYPE *h_C = (BASE_TYPE *)malloc(Csize); for (int i = 0; i < Arows * Acols; ++i) { h_A[i] = rand()/(BASE_TYPE)RAND_MAX; } for (int i = 0; i < Brows * Bcols; ++i) { h_B[i] = rand()/(BASE_TYPE)RAND_MAX; } BASE_TYPE *d_A = NULL; cudaMalloc((void **)&d_A, Asize); BASE_TYPE *d_B = NULL; cudaMalloc((void **)&d_B, Bsize); BASE_TYPE * d_C = NULL; cudaMalloc((void **)&d_C, Csize); cudaMemcpy(d_A, h_A, Asize, cudaMemcpyHostToDevice); cudaMemcpy(d_B, h_B, Bsize, cudaMemcpyHostToDevice); dim3 threadsPerBlock = dim3(BLOCK_SIZE, BLOCK_SIZE); dim3 blocksPerGrid = dim3(Bcols / BLOCK_SIZE, Arows / BLOCK_SIZE); 62 cudaEventRecord( start, 0); matrixMult<< KernelTime); cudaMemcpy(h_C, d_C, Csize, cudaMemcpyDeviceToHost); printf("Test STARTED\n"); for (int i = 0; i < Arows; i++) { for (int j = 0; j < Bcols; j++) { BASE_TYPE sum = 0; for (int k = 0; k < Acols; k++) sum += h_A[i * Acols + k] * h_B[k * Bcols + j]; if (fabs(h_C[i * Bcols + j] - sum) > 1e-3) { fprintf(stderr, "Result verification failed at element [%d, %d]!\n", i, j); printf("sum = %f, h_C[i * Bcols + j] = %f\n", sum, h_C[i * Bcols + j]); exit(EXIT_FAILURE); } } } printf("Test PASSED\n"); cudaFree(d_A); cudaFree(d_B); cudaFree(d_C); free(h_A); 63 free(h_B); free(h_C); cudaEventDestroy( start ); cudaEventDestroy( stop ); return 0; } Вопросы к разделу 1. Чем отличается создание матрицы на GPU от создания матрицы на CPU и последующего копирования ее элементов на GPU? 2. В каких участках кода при операциях над матрицами возникают потери в производительности? Лабораторная работа 1. Проверьте, является ли данная квадратная матрица ортогональной, путем умножения на транспонированную (не выделяя под нее дополнительной памяти) и сравнения результата с единичной матрицей. Домашняя работа 1. Проверьте, являются ли две заданные квадратные матрицы A и B коммутирующими (т.е. их произведение не зависит от порядка умножения: AB = BA). 2. Напишите функцию для сложения двух матриц M×N, заданных как вектор векторов (float a[M][N], b[M][N]). Матрицы задаются на хосте, после копируются на девайс, и там производится сложение. |