Лекторы: Боресков А.В. (ВМиК МГУ) Харламов А.А. (NVidia) Иерархия памяти CUDA. Разделяемая и константная память. Основные алгоритмы и решение СЛАУ на CUDA.

2011 Типы памяти в CUDA Тип памятиДоступУровень выделения Скорость работы Регистры R/WPer-threadВысокая(on-chip) Локальная R/WPer-threadНизкая (DRAM) Shared R/WPer-blockВысокая(on-chip) Глобальная R/WPer-gridНизкая (DRAM) Constant R/OPer-gridВысокая(L1 cache) Texture R/OPer-gridВысокая(L1 cache)

2011 Типы памяти в CUDA Самая быстрая – shared (on-chip) Самая медленная – глобальная (DRAM) Для ряда случаев можно использовать кэшируемую константную и текстурную память Доступ к памяти в CUDA идет отдельно для каждой половины warpа (half-warp)

Работа с shared-памятью zСамая быстрая (on-chip) zСейчас всего 16 Кбайт на один мультипроцессор zСовместно используется всеми нитями блока zОтдельное обращение для каждой половины warpа (half-warp) zКак правило, требует явной синхронизации

Типичное использование 1.Загрузить необходимые данные в shared-память (из глобальной) 2.__syncthreads () 3.Выполнить вычисления над загруженными данными 4.__syncthreads () 5.Записать результат в глобальную память

Умножение матриц (2) z При вычислении C постоянно используются одни и те же элементы из А и В z По много раз считываются из глобальной памяти z Эти многократно используемые элементы формируют полосы в матрицах А и В z Размер такой полосы N*16 и для реальных задач даже одна такая полоса не помещается в shared-память

Умножение матриц (2) z«Разделяй и властвуй» zРазбиваем каждую полосу на квадратные матрицы (16*16) zТогда требуемая подматрица произведения C может быть представлена как сумма произведений таких матриц 16*16 zДля работы нужно только две матрицы 16*16 в shared-памяти

Эффективная реализация __global__ void matMult ( float * a, float * b, int n, float * c ) { int bx = blockIdx.x, by = blockIdx.y; int tx = threadIdx.x, ty = threadIdx.y; int aBegin = n * BLOCK_SIZE * by; int aEnd = aBegin + n - 1; int bBegin = BLOCK_SIZE * bx; int aStep = BLOCK_SIZE, bStep = BLOCK_SIZE * n; float sum = 0.0f; for ( int ia = aBegin, ib = bBegin; ia

Эффективная реализация zНа каждый элемент z2*N арифметических операций z2*N/16 обращений к глобальной памяти zПоскольку разные warpы могут выполнять разные команды нужна явная синхронизация всех нитей блока zБыстродействие выросло более чем на порядок (2578 vs 132 миллисекунд)

Эффективная реализация zТеперь основное время (81.49%) ушло на вычисления zЧтение из памяти стало coalesced и заняло всего 12.5 %

Эффективная работа с shared- памятью zДля повышения пропускной способности вся shared-память разбита на 16 банков zКаждый банк работает независимо от других zМожно одновременно выполнить до 16 обращений к shared-памяти zЕсли идет несколько обращений к одному банку, то они выполняются по очереди

Эффективная работа с shared- памятью zБанки строятся из 32-битовых слов zПодряд идущие 32-битовые слова попадают в подряд идущие банки zBank conflict – несколько нитей из одного half-warpа обращаются к одному и тому же банку zКонфликта не происходит если все 16 нитей обращаются к одному слову (broadcast)

Бесконфликтные паттерны доступа

Паттерны с конфликтами банков zСлева – конфликт второго порядка – вдвое меньшая скорость zСправа - несколько конфликтов, до 6-го порядка

Доступ к массиву элементов __shared__ float a [N]; float x = a [base + threadIdx.x]; __shared__ short a [N]; short x = a [base + threadIdx.x]; __shared__ char a [N]; char x = a [base + threadIdx.x]; Нет конфликтов Конфликты 2-го порядка Конфликты 4-го порядка

Пример – матрицы 16*16 zПеремножение A*B T двух матриц 16*16, расположенных в shared-памяти yДоступ к обеим матрицам идет по строкам yВсе элементы строки распределены равномерно по 16 банкам yКаждый столбец целиком попадает в один банк zПолучаем конфликт 16-го порядка при чтении матрицы B

Пример – матрицы 16*16

z Дополним каждую строку одним элементом zВсе элементы строки (кроме последнего) лежат в разных банках z Все элементы столбца также лежат в разных банках z Фактически за счет небольшого увеличения объема памяти полностью избавились от конфликтов

Константная память __constant__ float contsData [256]; float hostData [256]; cudaMemcpyToSymbol ( constData, hostData, sizeof ( data ), 0, cudaMemcpyHostToDevice ); //////////////////////////////////////////////////////// template cudaError_t cudaMemcpyToSymbol ( const T& symbol, const void * src, size_t count, size_t offset, enum cudaMemcpyKind kind ); template cudaError_t cudaMemcpyFromSymbol ( void * dst, const T& symbol, size_t count, size_t offset, enum cudaMemcpyKind kind );

Базовые алгоритмы на CUDA Основные алгоритмы: zReduce zScan (prefix sum) zHistogram zSort

Параллельная редукция (reduce) zДано: yМассив элементов yБинарная ассоциативная операция + zНеобходимо найти yА, zЛимитирующий фактор – доступ к памяти zВ качестве операции также может быть min/max

Реализация параллельной редукции zКаждому блоку сопоставляем часть массива zБлок zКопирует данные в shared-память zИерархически суммирует данные в shared-памяти zСохраняет результат в глобальной памяти

Иерархическое суммирование zПозволяет проводить суммирование параллельно, используя много нитей zТребует log(N) шагов

Редукция, вариант 1

__global__ void reduce1 ( int * inData, int * outData ) { __shared__ int data [BLOCK_SIZE]; int tid = threadIdx.x; int i = blockIdx.x * blockDim.x + threadIdx.x; data [tid] = inData [i]; // load into shared memory __syncthreads (); for ( int s = 1; s < blockDim.x; s *= 2 ) { if ( tid % (2*s) == 0 ) // heavy branching !!! data [tid] += data [tid + s]; __syncthreads (); } if ( tid == 0 ) // write result of block reduction outData[blockIdx.x] = data [0]; }

Редукция, вариант 2

__global__ void reduce2 ( int * inData, int * outData ) { __shared__ int data [BLOCK_SIZE]; int tid = threadIdx.x; int i = blockIdx.x * blockDim.x + threadIdx.x; data [tid] = inData [i]; // load into shared memory __syncthreads (); for ( int s = 1; s < blockDim.x; s if ( index < blockDim.x ) data [index] += data [index + s]; __syncthreads (); } if ( tid == 0 ) // write result of block reduction outData [blockIdx.x] = data [0]; }

Редукция, вариант 2 zПрактически полностью избавились от ветвления zОднако получили много конфликтов по банкам yДля каждого следующего шага цикла степень конфликта удваивается

Редукция, вариант 3 zИзменим порядок суммирования yРаньше суммирование начиналось с соседних элементов и расстояние увеличивалось вдвое yНачнем суммирование с наиболее удаленных (на dimBlock.x/2) и расстояние будем уменьшать вдвое на каждой итерации цикла

Редукция, вариант 3

__global__ void reduce3 ( int * inData, int * outData ) { __shared__ int data [BLOCK_SIZE]; int tid = threadIdx.x; int i = blockIdx.x * blockDim.x + threadIdx.x; data [tid] = inData [i]; __syncthreads (); for ( int s = blockDim.x / 2; s > 0; s >>= 1 ) { if ( tid < s ) data [tid] += data [tid + s]; __syncthreads (); } if ( tid == 0 ) outData [blockIdx.x] = data [0]; }

Редукция, вариант 3 zИзбавились от конфликтов по банкам zИзбавились от ветвления zНо, на первой итерации половина нитей простаивает yПросто сделаем первое суммирование при загрузке

Редукция, вариант 4 __global__ void reduce4 ( int * inData, int * outData ) { __shared__ int data [BLOCK_SIZE]; int tid = threadIdx.x; int i = 2 * blockIdx.x * blockDim.x + threadIdx.x; data [tid] = inData [i] + inData [i+blockDim.x]; // sum __syncthreads (); for ( int s = blockDim.x / 2; s > 0; s >>= 1 ) { if ( tid < s ) data [tid] += data [tid + s]; __syncthreads (); } if ( tid == 0 ) outData [blockIdx.x] = data [0]; }

Редукция, вариант 5 zПри s

Редукция, вариант 5... for ( int s = blockDim.x / 2; s > 32; s >>= 1 ) { if ( tid < s ) data [tid] += data [tid + s]; __syncthreads (); } if ( tid < 32 ) // unroll last iterations { data [tid] += data [tid + 32]; data [tid] += data [tid + 16]; data [tid] += data [tid + 8]; data [tid] += data [tid + 4]; data [tid] += data [tid + 2]; data [tid] += data [tid + 1]; }

Редукция, быстродействие Вариант алгоритмаВремя выполнения (миллисекунды) reduction reduction reduction reduction49.10 reduction58.67

Parallel Prefix Sum (Scan) Имеется входной массив и бинарная операция По нему строится массив следующего вида

Parallel Prefix Sum (Scan) zОчень легко делается последовательно zДля распараллеливания используем sum tree zВыполняется в два этапа zСтроим sum tree zПо sum tree получаем результат

Построение sum tree

1.Используем одну нить на 2 элемента массива 2.Загружаем данные 3.__syncthreads () 4.Выполняем log(n) проходов для построения дерева

Построение sum tree #define BLOCK_SIZE 256 __global__ void scan1 ( float * inData, float * outData, int n ) { __shared__ float temp [2*BLOCK_SIZE]; inttid = threadIdx.x; intoffset = 1; temp [tid] = inData [tid];// load into shared memory temp [tid+BLOCK_SIZE] = inData [tid+BLOCK_SIZE]; for ( int d = n >> 1; d > 0; d >>= 1 ){ __syncthreads (); if ( tid < d ) { int ai = offset * (2 * tid + 1) - 1; int bi = offset * (2 * tid + 2) - 1; temp [bi] += temp [ai]; } offset

Получение результата по sum tree

1.Одна нить на 2 элемента массива 2.Обнуляем последний элемент 3.Copy and increment 4.Выполняем log(n) проходов для получения результата

Получение результата по sum tree if ( tid == 0 ) temp [n-1] = 0;// clear the last element for ( int d = 1; d < n; d = 1; __syncthreads (); if ( tid < d ) { int ai = offset * (2 * tid + 1) - 1; int bi = offset * (2 * tid + 2) - 1; float t = temp [ai]; temp [ai] = temp [bi]; temp [bi] += t; } __syncthreads (); outData [2*tid] = temp [2*tid];// write results outData [2*tid+1] = temp [2*tid+1];

Scan - Оптимизация Возможные проблемы: yДоступ к глобальной памяти -> coalesced yBranching -> small yКонфликты банков -> конфликты до 16 порядка !

Scan - Оптимизация Добавим по одному «выравнивающему» элементу на каждые 16 элементов в shared-памяти К каждому индексу добавим соответствующее смещение #define LOG_NUM_BANKS 4 #define CONLICT_FREE_OFFS(i) ((i) >> LOG_NUM_BANKS)

Scan - Оптимизация __shared__ float temp [2*BLOCK_SIZE+CONFLICT_FREE_OFFS(2*BLOCK_SIZE)]; inttid = threadIdx.x; intoffset = 1; intai = tid intbi = tid + (n / 2); intoffsA = CONFLICT_FREE_OFFS(ai); intoffsB = CONFLICT_FREE_OFFS(bi); temp [ai + offsA] = inData [ai + 2*BLOCK_SIZE*blockIdx.x]; temp [bi + offsB] = inData [bi + 2*BLOCK_SIZE*blockIdx.x]; for ( int d = n>>1; d > 0; d >>= 1, offset

Scan - Оптимизация if ( tid == 0 ){ inti = n CONFLICT_FREE_OFFS(n-1); sums [blockIdx.x] = temp [i];// save the sum temp [i] = 0;// clear the last element } for ( int d = 1; d < n; d = 1; __syncthreads (); if ( tid < d ){ intai = offset * (2 * tid + 1) - 1; intbi = offset * (2 * tid + 2) - 1; floatt; ai += CONFLICT_FREE_OFFS(ai); bi += CONFLICT_FREE_OFFS(bi); t = temp [ai]; temp [ai] = temp [bi]; temp [bi] += t; } __syncthreads (); outData [ai + 2*BLOCK_SIZE*blockIdx.x] = temp [ai + offsA]; outData [bi + 2*BLOCK_SIZE*blockIdx.x] = temp [bi + offsB];

Scan больших массивов Рассмотренный код хорошо работает для небольших массивов, целиком, помещающихся в shared-память В общем случае: yВыполняем отдельный scan для каждого блока yДля каждого блока запоминаем сумму элементов (перед обнулением) yПрименяем scan к массиву сумм yК каждому элементу, кроме элементов 1-го блока добавляем значение, соответствующее данному блоку

Scan больших массивов voidscan ( float * inData, float * outData, int n ) { int numBlocks = n / (2*BLOCK_SIZE); float * sums, * sums2; if ( numBlocks < 1 ) numBlocks = 1; // allocate sums array cudaMalloc ( (void**)&sums, numBlocks * sizeof ( float ) ); cudaMalloc ( (void**)&sums2, numBlocks * sizeof ( float ) ); dim3 threads ( BLOCK_SIZE, 1, 1 ), blocks ( numBlocks, 1, 1 ); scan3 >> ( inData, outData, sums, 2*BLOCK_SIZE ); if ( n >= 2*BLOCK_SIZE ) scan ( sums, sums2, numBlocks ); else cudaMemcpy ( sums2, sums, numBlocks*sizeof(float), cudaMemcpyDeviceToDevice ); threads = dim3 ( 2*BLOCK_SIZE, 1, 1 ); blocks = dim3 ( numBlocks - 1, 1, 1 ); scanDistribute >> ( outData + 2*BLOCK_SIZE, sums2 + 1 ); cudaFree ( sums ); cudaFree ( sums2 ); }

Построение гистограммы yДан массив элементов и способ классификации элементов: каждому элементу сопоставляется один из k классов. yЗадача – по массиву получить для каждого класса число элементов, попадающих в него. yПолученная таблица частот для классов и является гистограммой

Построение гистограммы zОчень легко реализуется последовательным кодом zЕсли мы выделяем по одной нити на каждый входной элемент, то нужна операция atomicIncr zОчень частые обращения к счетчикам, лучше всего их разместить в shared- памяти zИдеальный случай – у каждой нити своя таблица счетчиков в shared-памяти

Построение гистограммы для 64 классов (bins) yНа каждый счетчик отводим 1 байт yВсего гистограмма – 64 байта (на одну нить) yВсего в разделяемой памяти SM можно разместить 256 таких гистограмм yРазмер блока – 64 нити, максимум 4 блока на SM yКаждая нить может обработать не более 255 байт

Построение гистограммы для 64 классов (bins) Посмотрим на конфликты банков: y64*tid+value ybank = ((64*tid+value)/4) & 0xF=(value>>2) & 0xF yНомер банка полностью определяется входными данными, если есть много повторений, то будет много конфликтов по банкам y64*value+tid ybank = ((64*value+tid)/4) & 0xF=(tid>>2) & 0xF yНомер банка определяется номером нити

Построение гистограммы для 64 классов (bins) yВ первом случае все определяется входными данными, очень высока вероятность конфликта банков вплоть до 16-го порядка. yВо втором случае номер банка определяется старшими битами номера нити и мы получаем постоянный конфликт четвертого порядка yНо зачем в качестве tid использовать именно номер нити в блоке – подойдет любое значение, получаемое из номера нити путем фиксированной перестановки битов

Построение гистограммы для 64 классов (bins) yНомер банка определяется битами 2..5 величины tid. yПостроим tid как следующую перестановку битов номера нити в блоке: ytid=(threadIdx.x>>4) | ((threadIdx.x & 0xF)

Построение гистограммы для 64 классов (bins) inline __device__ void addByte ( uchar * base, uint data ) { base[64*data]++; } inline __device__ void addWord ( uchar * base, uint data ) { addByte ( base, (data >> 2) & 0x3FU ); addByte ( base, (data >> 10) & 0x3FU ); addByte ( base, (data >> 18) & 0x3FU ); addByte ( base, (data >> 26) & 0x3FU ); } __global__ void histogram64Kernel ( uint * partialHist, uint * data, uint dataCount ) { __shared__ uchar hist [64*64]; int tid = (threadIdx.x >> 4) | ((threadIdx.x & 0x0F)

Построение гистограммы zБолее общий случай– zПросто не хватит shared-памяти давать каждой нити по своей гистограмме zДавайте выделять по своей таблице счетчиков на определенный набор нитей z(+) Уменьшаются затраты на shared-память z(-) Появляется проблема синхронизации с записью нитей этого набора

Построение гистограммы Когда проще всего обеспечивать атомарность записи: zКогда каждый такой набор нитей всегда лежит в пределах одного warpа zПо-прежнему сохраняется риск нескольких нитей, одновременно увеличивающих один и тот же элемент гистограммы, но с этим можно бороться zЕсли несколько нитей одновременно делают запись по одному адресу, то только одна из этих записей проходит

Построение гистограммы zПусть каждый warp нитей имеет свою таблицу счетчиков z192 нити в блоке дают 6 warpов, т.е. 6*256*4=6Кб shared-памяти на блок z5 старших битов каждого счетчика будут хранить номер нити (внутри warpа), сделавшей последнюю запись

Построение гистограммы __device__ inline void addData256 ( volatile unsigned * warpHist, unsigned data, unsigned threadTag ) { unsigned count; do { count = warpHist [data] & 0x07FFFFFFU; // mask thread tag bits count = threadTag | (count + 1); // increment count and tag it warpHist [data] = count; } while ( warpHist [data] != count ); // check whether we've modified value } zКаждая нить строит новое значение zУвеличить на единицу zВыставить старшие биты в номер нити в warpе zКак минимум одна запись пройдет и соответствующая нить выйдет из цикла

Построение гистограммы zКаждая нить меняет свой элемент таблицы zСразу же выходим, никаких расходов zДве нити пытаются увеличить один и тот же счетчик zОдной это получится (запишется ее значение) zДругой нет – ее значение будет отброшено zТа нить, которая записала выходит из цикла и оставшаяся нить со делает запись (со второй попытки)

Построение гистограммы #define WARP_LOG2SIZE 5// bits to identify warp #define WARP_N 6// warps per block __global__ void histogramKernel ( unsigned * result, unsigned * data, int n ) { int globalTid = blockIdx.x * blockDim.x + threadIdx.x; int numThreads = blockDim.x * gridDim.x; int warpBase = (threadIdx.x >> WARP_LOG2SIZE) * BIN_COUNT; unsigned threadTag = threadIdx.x > 0) & 0xFFU, threadTag ); addData256 ( hist + warpBase, (data4 >> 8) & 0xFFU, threadTag ); addData256 ( hist + warpBase, (data4 >> 16) & 0xFFU, threadTag ); addData256 ( hist + warpBase, (data4 >> 24) & 0xFFU, threadTag ); } __syncthreads(); for ( int i = threadIdx.x; i < BIN_COUNT; i += blockDim.x ){ unsigned sum = 0; for ( int base = 0; base < BLOCK_MEMORY; base += BIN_COUNT ) sum += hist [base + i] & 0x07FFFFFFU; result[blockIdx.x * BIN_COUNT + i] = sum; }

Сортировка. Битоническая сортировка Базовая операция – полуочиститель, упорядочивающий пары элементов на заданном расстоянии:

Битоническая сортировка zПоследовательность называется битонической, если она zСостоит из двух монотонных частей zПолучается из двух монотонных частей циклическим сдвигом zПримеры: z1,3,4,7,6,5,2 z5,7,6,4,2,1,3 (получена сдвигом 1,3,5,7,6,4,2)

Битоническая сортировка zЕсли к битонической последовательности из n элементов применить полуочиститель Bn, то в результате у полеченной последовательности zОбе половины будут битоническими zЛюбой элемент первой половины будет не больше любого элемента второй половины zХотя бы одна из половин будет монотонной

Битоническая сортировка Если к битонической последовательности длины n применить получистители B n,B n/2,…,B 8,B 4,B 2 то в результате мы получим отсортированную последовательность (битоническое слияние)! Если у нас произвольная последовательность: zПрименим В2 с чередующимся порядком, в результате каждые 4 подряд идущих элемента будут образовывать битоническую последовательность zПри помощи битонического слияния отсортируем каждую пару последовательностей из 4-элементов, но с чередующимся порядком упорядочивания. zПри помощи битонического слияния отсортируем каждую пару из 8 элементов

Битоническая сортировка Пусть есть произвольная последовательность длины n. Применим к каждой паре элементов полуочиститель B 2 с чередующимся порядком сортировки. Тогда каждая четверка элементов будет образовывать битоническую последовательность.

Битоническая сортировка Применим к каждой такой четверке элементов полуочиститель B 4 с чередующимся порядком сортировки. Тогда каждые восемь элементов будет образовывать битоническую последовательность. Применим к каждым 8 элементам полуочиститель B 4 с чередующимся порядком сортировки и так далее. Всего потребуется log(n)*log(n) проходов для полной сортировки массива

Битоническая сортировка zОчень хорошо работает для сортировки через шейдеры (steps3d.narod.ru/tutorials/gpu-sort- tutorial.html) zПлохо использует возможности CUDA, поэтому обычно не используется для сортировки больших массивов

Поразрядная сортировка (radix sort) Пусть задан массив из 32- битовых целых чисел: Отсортируем этот массив по старшему (31-му) биту, затем по 30-му биту и т.д. После того, как мы дойдем до 0-го бита и отсортируем по нему, последовательность будет отсортирована

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

Поразрядная сортировка Пусть нам надо отсортировать массив по k-му биту. Тогда рассмотрим массив, где из каждого элемента взят данный бит (b[i]=(a[i] >> k) & 1). Каждый элемент этого массива равен или нулю или единице. Применим к нему операцию scan, сохранив при этом общую сумму элементов

Поразрядная сортировка В результате мы получим сумму всех выбранных бит (т.е.число элементов исходного массива, где в рассматриваемой позиции стоит единичный бит) и массив частичных сумм битов s n Отсюда легко находится количество элементов исходного массива, где в рассматриваемой позиции стоит ноль (Nz). По этим данным легко посчитать новые позиции для элементов массива:

Поразрядная сортировка По этим данным легко посчитать новые позиции для элементов массива:

Поразрядная сортировка - float Поразрядная сортировка легко адаптируется для floating point-величин. Положительные значения можно непосредственно сортировать Отрицательные значения при поразрядной сортировке будут отсортированы в обратном порядке

Поразрядная сортировка - float uint flipFloat ( uint f ) { uint mask = -int(f >> 31) | 0x ; return f ^ mask; } uint unflipFloat ( uint f ) { uint mask = ((f >> 31) - 1) | 0x ; return f ^ mask; } Чтобы сортировать значения разных знаков достаточно произвести небольшое преобразование их тип uint, приводимое ниже

2011 Решение системы линейных алгебраических уравнений Традиционные методы ориентированы на последовательное вычисление элементов и нам не подходят Есть еще итеративные методы Ax=f, A – матрица размера N*N, f – вектор размера N

2011 Итеративные методы zЭффективны когда zМатрица А сильна разрежена zПараллельные вычисления zВ обоих случаях цена (по времени) одной итерации O(N)

2011 Сходимость zЕсли есть сходимость, то только к решению системы zЗаписав уравнения для погрешности получаем достаточное условие сходимости zЗа счет выбора достаточно малого значения параметра получаем сходимость

2011 Код на CUDA // // one iteration // __global__ void kernel ( float * a, float * f, float alpha, float * x0, float * x1, int n ) { int idx = blockIdx.x * blockDim.x + threadId.x; int ia = n * idx; float sum = 0.0f; for ( int i = 0; i < n; i++ ) sum += a [ia + i] * x0 [i]; x1 [idx] = x0 [idx] + alpha * (sum – f [idx] ); }

2011 Ресуры нашего курса CUDA.CS.MSU.SU – Место для вопросов и дискуссий – Место для материалов нашего курса – Место для ваших статей! Если вы нашли какой-то интересный подход! Или исследовали производительность разных подходов и знаете, какой из них самый быстрый! Или знаете способы сделать работу с CUDA проще!

2011 Вопросы