ТЕХНОЛОГИЯ ПРОГРАММИРОВАНИЯ CUD. Учебное пособие казань 2017 2 удк 004. 4, 681. 3 Ббк 32. 81 Печатается по постановлению Редакционноиздательского совета
 Скачать 1.28 Mb.
|
|
Глава 7. Типы памяти. Разделяемая память В CUDA предусмотрены несколько типов памяти, каждый из которых предназначен под определенные нужды. Ниже приведена таблица типов памяти GPU. Таблица 3. Типы памяти GPU Тип памяти Доступ Выделяется Скорость работы Регистры чтение/запись на поток Высокая Локальная чтение/запись на поток Низкая Разделяемая чтение/запись на блок Высокая Глобальная чтение/запись на все решетки Низкая Константная чтение на все решетки Низкая Текстурная чтение на все решетки Низкая Регистры используется для хранения локальных переменных. Скорость доступа для данного типа памяти самая быстрая, и выделение происходит на каждую нить. Локальная память используется для хранения локальных переменных, когда регистров не хватает, скорость доступа низкая. Память выделяется так же, как и для регистров. Разделяемая память используется для хранения массивов данных, используемых совместно всеми нитями в блоке. Эта память имеет чуть меньшую скорость доступа, чем регистры, и выделяется на блок. В глобальной памяти хранятся исходно все данные, это общая память для всех блоков. Константная память используется для передачи параметров в ядро, превышающих допустимые размеры для параметров ядра, выделение целиком на решетку. Текстурная память используется для хранения больших массивов данных, выделяется целиком на решетку. Обращение к памяти в CUDA осуществляется одновременно из половины количества нитей в блоке. Все типы памяти, кроме регистровой и локальной, 65 можно использовать как эффективно, так и не эффективно. Поэтому необходимо четко понимать их особенности и области применения. Рассмотрим более подробно разделяемую память. Основной особенностью разделяемой памяти, является то, что она оптимизирована под обращение к ней одновременно половины количества нитей в блоке (16 нитей). Вся разделяемая память разбита на 32 банка, состоящих из 32 битовых слов, причем, последовательно идущие слова попадают в последовательно идущие блоки. Обращение к каждому банку происходит независимо. Второй важной особенностью разделяемой памяти является то, что она выделяется на блок. Для задания разделяемой памяти используется спецификатор __shared__. Часто использование разделяемой памяти требует явной синхронизации внутри блока. При работе с разделяемой памятью необходимо обеспечить синхронизацию потоков, что в CUDA довольно просто – она предоставляет примитивный барьер, функцию __syncthreads(). При ее вызове потоки приостанавливаются до того момента, пока все потоки не достигнут этой точки. Важно, чтобы вызов __syncthreads() у всех потоков находился в одном и том же месте (в частности, был безусловным, т. е. не находился внутри оператора if и т.п.), в противном случае результаты могут быть непредсказуемыми, в частности, может возникнуть "мертвая блокировка" (deadlock). Приведем фрагмент кода создания массива данных в разделяемой памяти. #define BS 256 //Размер блока __global__ void Kernel (float* data) { // Создание массива в разделяемой памяти __shared__ float a[BS]; int i = blockIdx.x * BS + threadIDx.x; // Копирование из глобальной памяти в разделяемую a[i] = data[i]; // Синхронизируем нити. 66 // Перед использованием нужно быть уверенным, что все данные скопированились __syncthreads(); // Использонание данных разделяемой памяти } Данный пример демонстрирует только работу с разделяемой памятью. Для того, чтобы понять, как же разделяемая память помогает в ускорении работы напишем программу для вычисления скалярного произведения векторов. Попробуем написать эту программу с использованием глобальной памяти. Используя механизм событий, замерим время работы ядра. Рассмотрим, как изменится ядро, если, вместо глобальной памяти, использовать разделяемую память. __global__ void scalMult(const BASE_TYPE *A, const BASE_TYPE *B, BASE_TYPE *C, int numElem) { // Переменная для хранения суммы элементов BASE_TYPE sum = 0.0; // Создание массивов в разделяемой памяти __shared__ BASE_TYPE ash[BLOCK_SIZE]; __shared__ BASE_TYPE bsh[BLOCK_SIZE]; // Копирование из глобальной памяти ash[threadIdx.x] = A[blockIdx.x * blockDim.x + threadIdx.x]; bsh[threadIdx.x] = B[blockIdx.x * blockDim.x + threadIdx.x]; // Синхронизация нитей __syncthreads(); // Вычисление скалярного произведения if (threadIdx.x == 0) { sum = 0.0; for (int j = 0; j < blockDim.x; j++) sum += ash[j] * bsh[j]; 67 C[blockIdx.x] = sum; } } Заменим ядро программы и оценим время работы ядра, использующего разделяемую память. Так как каждая нить теперь работает с памятью, чья скорость работы более высокая, большую часть времени работы программы занимают вычисления, а не доступ к памяти. Более наглядно это можно наблюдать в случае перемножения матриц. Банки данных Доступ к памяти осуществляется через 32 банка данных. Причем номер банка можно посчитать по формуле Номер банка = (Адрес в байтах / 4) % 32. Каждый банк в определенный момент времени обрабатывает одну транзакцию (чтение или запись). Оптимальный доступ к данным обеспечивается, когда каждая нить последовательно соотносится со своим банком (см. Рис. 6а). Также конфликта между банками нет, если каждой нити в любом порядке соответствует свой банк. Пример такого доступа изображен на Рис. 6б. Конфликт имеет место, когда несколько нитей пытаются обратиться к одному банку данных, как изображено на Рис. 6в. Однако, в случае, когда все нити записывают или считывают с одного банка, то конфликта не происходит. Рис. 6а. Оптимальный доступ к банкам данных 68 Рис. 6б. Бесконфликтный доступ к банкам данных Рис. 6в. Конфликтный доступ к банкам данных Рассмотрим наиболее распространенные случаи, когда могут иметь место конфликты банков данных. В следующем примере __shared__ char s[N]; char data = s[threadIdx.x]; имеем однобайтовый массив, и, так как в одном банке содержится четыре байта, то одновременно в один банк будут обращаться четыре нити. Это приведет к конфликту четвертого порядка (4 нити на 1 банк). Отметим, что на Рис. 6в приведен пример конфликтов шестого порядка. Один из способов решения такого конфликта состоит в добавлении избыточных данных: __shared__ char s[4 * N]; char data = s[threadIdx.x * 4]; В следующем примере будет конфликт второго порядка (short – двухбайтовый тип) 69 __shared__ short s[N]; short data = s[threadIdx.x]; Чаще возникают конфликты, когда используются четырехбайтовые числа, и номер в массиве вычисляется по формуле. Например, в следующем коде: __shared__ float s[N]; float data = s[threadIdx.x * L]; при четных L возникают конфликты. Действительно, при L = 2 нулевая нить и 16-натая нить обращаются в нулевой банк, а при L = 4 уже нити с номерами 0, 8, 16 и 24 обращаются в нулевой банк, создавая конфликт четвертого порядка. Пример. Оптимизация программы по перемножению двух матриц Пусть даны две матрицы размером 𝑁𝑁 ∗ 𝑁𝑁 (для удобства будем считать, что 𝑁𝑁 кратно 16). Если решать данную задачу, получая данные постоянно из глобальной памяти, то понадобится 2𝑁𝑁 чтений из глобальной памяти. Становится ясно, что основным тормозящим программу фактором является не что иное, как чтение из глобальной памяти, которая имеет низкую скорость передачи. Заметно ускорить программу поможет правильная организация вычислений и разделяемая память. Обратим внимание на то, что для вычисления элемента новой матрицы необходимы лишь небольшие части исходных матриц. Для повышения быстродействия программы разобьем новую матрицу на подматрицы таким образом, чтобы обработкой данной подматрицы занимался только один блок. В данной задаче разобьем на подматрицы (𝐶𝐶′) размером 16*16 элементов (см. Рис. 7). 70 Рис. 7. Вычисление подматрицы новой матрицы Однако, из-за малого объема разделяемой памяти, подматрицы 𝐴𝐴′ и 𝐵𝐵′ не могут быть целиком скопированы. Поступим с 𝐴𝐴′ и 𝐵𝐵′ так же, как поступили с подматрицей 𝐶𝐶′, разобьем каждую матрицу на подматрицы размером 16*16. В итоге вычисление 𝐶𝐶′ будем проводить в 𝑁𝑁/16 шагов (см. Рис. 8). Вычисление подматрицы C' можно записать следующим образом: 𝐶𝐶 ′ = 𝐴𝐴 1 ′ ∗ 𝐵𝐵 1 ′ + 𝐴𝐴 2 ′ ∗ 𝐵𝐵 2 ′ + ⋯ + 𝐴𝐴 𝑁𝑁 16 ′ ∗ 𝐵𝐵 𝑁𝑁 16 ′ Теперь на каждом шаге будем подгружать в разделяемую память по одной подматрице 𝐴𝐴 𝑖𝑖 ′ и 𝐵𝐵 𝑖𝑖 ′ , вычислять соответствующую им сумму для элементов произведения и по окончанию всех вычислений записывать в результирующую матрицу. Укажем в начале программы размер блока и тип элементов массива: // размер блока или размер подматрицы #define BLOCK_SIZE 16 //тип, который будут иметь элементы матриц #define BASE_TYPE double 71 Рис. 8. Разбиение исходных матриц Учитывая все вышесказанное об организации вычислений, ядро будет выглядеть следующим образом: __global__ void matrixMult(const BASE_TYPE *A, const BASE_TYPE *B, BASE_TYPE *C, int Acols, int Bcols) { // индекс начала первой подматрицы А, которую // обрабатывает блок int aBegin = Acols * blockDim.y * blockIdx.y; // индекс конца подматрицы А, которую обрабатывает блок int aEnd = aBegin + Acols - 1; // шаг для перебора подматриц А int aStep = blockDim.x; // индекс начала первой подматрицы В, которую // обрабатывает блок int bBegin = blockDim.x * blockIdx.x; // шаг для перебора подматриц В int bStep = blockDim.y * Bcols; // Выделение разделяемой памяти для подматриц __shared__ BASE_TYPE as [BLOCK_SIZE][BLOCK_SIZE]; __shared__ BASE_TYPE bs [BLOCK_SIZE][BLOCK_SIZE]; 72 // переменная для вычисления элемента подматрицы BASE_TYPE sum = 0.0; for (int ia = aBegin, ib = bBegin; ia < aEnd; ia += aStep, ib += bStep) { // загрузка подматриц А и В из глобальной памяти в // разделяемую as[threadIdx.y][threadIdx.x] = A[ia + Acols * threadIdx.y + threadIdx.x]; bs[threadIdx.y][threadIdx.x] = B[ib + Bcols * threadIdx.y + threadIdx.x]; // синхронизация нитей __syncthreads(); // перемножение двух матриц for (int k = 0; k < blockDim.x; k++) sum += as[threadIdx.y][k] * bs[k][threadIdx.x]; // синхронизация нитей __syncthreads(); } // индекс результирующего элемента в глобальной памяти int ind = Bcols * (blockDim.y * blockIdx.y + threadIdx.y) + blockDim.x * blockIdx.x + threadIdx.x; // запись элемента в глобальную память C[ind] = sum; } После любой работы с разделяемой памятью необходимо синхронизировать нити. В данном примере важно произвести синхронизацию до момента начала расчетов, чтобы убедиться в том, что все данные загружены, 73 а также после расчетов, чтобы убедиться в том, что текущие подматрицы не нужны, и можно подгружать новые подматрицы. Вызов ядра из главной функции: // параметры запуска ядра dim3 threadsPerBlock = dim3(BLOCK_SIZE, BLOCK_SIZE); dim3 blocksPerGrid = dim3(Bcols / BLOCK_SIZE, Arows / BLOCK_SIZE); // старт времени работы ядра cudaEventRecord( start, 0); // вызов ядра matrixMult<< // стоп времени работы ядра cudaEventRecord( stop, 0); cudaEventSynchronize( stop ); // вычисление времени работы float KernelTime; cudaEventElapsedTime( &KernelTime, start, stop); Использование такой алгоритм для перемножении двух матриц требует 2 * N / 16 чтений из глобальной памяти. Понятно, что правильное использование разделяемой памяти может поднять быстродействие программы даже более чем на порядок. Также метод вычисления суммы произведений элементов по блокам, который был применен в ядре для вычисления произведения двух матриц, ускоряет вычисления и часто применяется в подобных задачах. Если размер разделяемой памяти не задан изначально, и определяется в процессе работы программы, то необходимо его указывать при вызове ядра 74 myKernelFunc<< Вопросы к разделу 1. Какие типы памяти выделяют в технологии CUDA? 2. Какой спецификатор применяется для задания разделяемой памяти? Перечислите основные особенности разделяемой памяти. Лабораторная работа 1. Напишите программу для скалярного умножения двух векторов, используя разделяемую память. Разделяемая память выделяется с учетом ограничения ее размеров на блок. Количество нитей при вызове ядра равно числу элементов массива с разделяемой памятью. Замерьте время работы ядра программы. 2. Напишите программу вычисления определенного интеграла функции одной переменной по квадратурной формуле центральных прямоугольников. Домашняя работа 1. Доработайте программу вычисления определенного интеграла, используя другие квадратурные формулы: формулы трапеций и Симпсона. Сравните точность результатов. 2. Напишите программу вычисления евклидовой нормы вектора (путем скалярного умножения его на себя и вычисления корня из произведения). 75 Глава 8. Константая и текстурная память Помимо использования разделяемой памяти существует еще два вида памяти, которые можно использовать в CUDA – это константная и текстурная память. Эти виды памяти являются медленнее разделяемой, но за счет кэширования они превосходят по скорости глобальную память. Константная память Константная память используется тогда, когда в ядро необходимо передать много различных данных, которые будут одинаково использоваться всеми нитями ядра. Эта память выделяется в глобальной памяти в пределах 64 кБ. При этом кэширование происходит отдельного для каждого мультипроцессора, размер кэша равен 8 кБ для всех версий производительности кроме 6.0, для 6.0 объем кэша составляет 4 кБ. Объявление переменной делается глобально, используя слово __constant__. Пример использования константного массива: __device__ __constant__ float contsData[4]; float hostData[4]; …… // присваивание значений в hostData // копирование данных с центрального процессора в // константную память cudaMemcpyToSymbol (constData, hostData, 4*sizeof(float), 0, cudaMemcpyHostToDevice); Слово __device__ перед __constant__ можно опускать, так как по умолчанию константная память выделяется на графическом устройстве. Использование внутри ядра константной памяти ничем не отличается от использования любой глобальной переменной. Рассмотрим следующий пример использования константы: 76 __constant__ int add; // функция сложения двух векторов с прибавлением // постоянной add __global__ void addKernel(int *c, const int *a, const int *b) { int i = threadIdx.x; c[i] = a[i] + b[i] + add; } …… // инициализация и вызов функции сложения в основной // программе, size – размер массивов int r = 2; cudaMemcpyToSymbol(add, &r, sizeof(int), 0, cudaMemcpyHostToDevice); addKernel<<<1, size>>>(dev_c, dev_a, dev_b); Текстурная память При рассмотрении разделяемой памяти, мы предполагали, что использование памяти внутри блока сильно локализовано, то есть что блоку для его работы требуется лишь какой-то определенный участок массива, если же локализовать обращение не удается, а скорость обращения к памяти является узким местом программы, можно попробовать использовать текстурную память. Основная особенность текстурной памяти – это использования кэша. Однако существуют дополнительные стадии конвейера (преобразование адресов, фильтрация, преобразование данных), которые снижают скорость первого обращения. Поэтому текстурную память разумно использовать с следующих случаях: • Объем данных не влезает в shared память; 77 • Паттерн доступа хаотичный; • Данные переиспользуются разными потоками. Для использования текстурной памяти необходимо задать объявление текстуры как глобальную переменную: texture < type, dim, tex_type> g_TexRef ; Type – тип хранимых переменных, dim – размерность текстуры (1, 2, 3) Tex_type – тип возвращаемых значений • cudaReadModeNormalizedFloat, • cudaReadModeElementType. Кроме того, для более полного использования возможностей текстурной памяти можно задать описание канала: struct cudaChannelFormatDesc { int x , y , z , w; enum cudaChannelFormatKind f ; } ; Задает формат возвращаемого значения int x, y, z, w – числа в интервале [0, 32], проекция исходного значения по битам cudaChannelFormatKind – тип возвращаемого значения • cudaChannelFormatKindSigned – знаковые int • cudaChannelFormatKindUnsigned – беззнаковые int • cudaChannelFormatKindFloat – float Текстурная память, как и константная, кэшируется. Причем, размер кэша зависит от производительности. До производительности 6.0 размер кэша изменяется от 12 до 48 КБ, для 6.x – от 24 до 48 КБ, для 7.0 – от 32 до 128 КБ. 78 Примеры использования одномерной текстуры Приведем пример программы для работы с одномерной текстурой, привязанной к CUDA массиву. #include "cuda_runtime.h" #include "device_launch_parameters.h" #include // объявляем одномерный текстурный объект texture __global__ void readTexels(int n, float *d_out) { int idx = blockIdx.x * blockDim.x + threadIdx.x; if (idx < n) { float x = tex1D(texRef, float(idx)); d_out[idx] = x; } } #define NUM_THREADS 256 int main() { int N = 10; int nBlocks = N / NUM_THREADS + ((N % NUM_THREADS) ? 1 : 0); float *d_out; // выделяем память на GPU cudaMalloc((void**)&d_out, sizeof(float) * N); // выделяем память на CPU float *h_out = (float*)malloc(sizeof(float)*N); // заполняем массив data float *data = (float*)malloc(N * sizeof(float)); 79 for (int i = 0; i < N; i++) data[i] = float(2 * i); // создаем CUDA массив на device cudaArray* cuArray; cudaMallocArray(&cuArray, &texRef.channelDesc, N, 1); cudaMemcpyToArray(cuArray, 0, 0, data, sizeof(float)*N, cudaMemcpyHostToDevice); // привязываем текстуру к CUDA массиву cudaBindTextureToArray(texRef, cuArray); // устанавливаем атрибуты текстуры texRef.normalized = false; texRef.filterMode = cudaFilterModePoint; readTexels << // копируем с девайса на хост cudaMemcpy(h_out, d_out, sizeof(float)*N, cudaMemcpyDeviceToHost); for (int i = 0; i < N; i++) printf("%f\n", h_out[i]); free(h_out); cudaFree(d_out); cudaFreeArray(cuArray); cudaUnbindTexture(texRef); getchar(); } Программа выделяет память на GPU, привязывает к ней текстуру. Затем копирует из текстуры в массив. Ответ должен быть следующим: 0.000000 2.000000 80 4.000000 6.000000 8.000000 10.000000 12.000000 14.000000 16.000000 18.000000 Следующий пример кода связывает текстуру с линейной памятью. Эта программа заполняет на девайсе обыкновенной массив из текстуры, после этого копирует его на хост и сравнивает с исходным массивом. Программа должна выдать сообщение «Массивы совпали». #include "cuda_runtime.h" #include "device_launch_parameters.h" #include #include { cudaError_t err = cudaGetLastError(); if (cudaSuccess != err) { fprintf(stderr, "Cuda error: %s: %s.\n", msg, cudaGetErrorString(err)); exit(EXIT_FAILURE); } } texture __global__ void kernel(int n, float *d_out) { int idx = blockIdx.x*blockDim.x + threadIdx.x; if (idx < n) 81 d_out[idx] = tex1Dfetch(texRef, idx); } #define NUM_THREADS 256 int main() { int N = 5000; int nBlocks = N / NUM_THREADS + ((N % NUM_THREADS) ? 1 : 0); int memSize = N * sizeof(float); // data fill array with increasing values float *data; data = (float*)malloc(memSize); for (int i = 0; i < N; i++) data[i] = float(i); float *d_a; cudaMalloc((void **)&d_a, memSize); cudaMemcpy(d_a, data, memSize, cudaMemcpyHostToDevice); cudaBindTexture(0, texRef, d_a, memSize); checkCUDAError("bind"); kernel << 82 free(h_out); free(data); getchar(); } Вопросы к разделу 1. Где выделяется константная и текстурная память? 2. Где расположен кэш для константной и текстурной памяти? Какие размеры кэша? Лабораторная работа 1. Напишите программу для скалярного умножения двух векторов, используя константную память. Замерьте время работы ядра программы. 2. Напишите программу вычисления определенного интеграла функции одной переменной по квадратурной формуле центральных прямоугольников. Используйте тектстурную память, привязанную к линейной памяти. Домашняя работа 1. Напишите программу для скалярного умножения двух векторов, расположенных в текстурной памяти. Испольуйте три алгоритма: 1. В обоих случаях привязка к CUDA массиву; 2. В одном – к CUDA массиву, в другом – к линейной памяти; 3. В обоих случаях привязаны к линейной памяти. Сравните время работы для различных размеров массивов. |