cuda_course_task_matmul (Задания)

Лихогруд Николайn.lihogrud@gmail.comЗаданиеПостановкаm×n × n×p = m×p = ==1< , >Постановка Нить (, )вычисляет элементмассива ,Грид Размер блока – 32x32 Понадобитсявертикали,3232блоков попо горизонтали Если m и p не кратны 32 – будутпростаивающие нити, это нестрашноЯдро__global__ void matmul(A,B,C) {< Вычислить глобальные индексы нити (idx, idy) в матрице C >;if (idx < p && idy < m) {< вычислить скалярное произведение строки idyматрицы A на столбец idx матрицы B >;< записать результат скалярного произведения в элемент[idy, idx] матрицы C >;}}Глобальные индексы Каковы глобальные индексы нити (idx, idy), если индексы нити в блоке - (threadIdx.x = 22, threadIdx.y = 3), индексы блока в гриде - (blockIdx.x = 3, blockIdx.y = 1), размер каждого блока - (blockDim.x = blockDim.y = 32)???ПростаиваютИндексы CUDA и матричные Очень вероятно, что вы будете путаться между (i,j) и (x,y): i – номер строки, соответствует координате по y в терминологииCUDA j – номер столбца, соответствует координате по x Путаница будет оттого что элементу [i,j] в матричной нотациисоответсвует элемент [y, x] в индексах CUDA (пространственных), но вCUDA первым пишется индекс x Чтобы не путаться, можно глобальные индексы сразу назвать [i,j], ане [idx, idy]Структура cudaPitchedPtr Определена в CUDA Toolkit:struct cudaPitchedPtr {size_t pitch; // число байтов между началами строкvoid *ptr;// указатель на памятьsize_t xsize;// логическая ширина матрицы в элементахsize_t ysize;// логическая высота матрицы в элементах} cudaPitchedPtr make_cudaPitchedPtr(void *d, size_t p,size_t xsz, size_t ysz) создать такую структуруОбертка матриц матриц Вместе с матрицей часто передают её размеры Следует обернуть все матрицы в cudaPitchedPtr С обернутыми матрицами легче будет сделать следующую часть Матрицы плотные, расположены по строкам Число байтов между началами строк (pitch) =<число столбцов> * <размер элемента> Для матрицы с элементами типа T из m строк и n столбцов:cudaPitchedPtr C =make_cudaPitchedPtr(malloc(n * m * sizeof(T)),n * sizeof(T), n, m);get_elem Обращение к элементам матрицы, заданной при помощи cudaPitchedPtr:Type a = ((Type *)((char*)matrix.ptr + Row * matrix.pitch)) [Column]; Удобно определить макрос:#define get_elem(array, Row, Column) \(((Type*)((char*)array.ptr + (Row) * array.pitch))[(Column)])Type a = get_elem(array, 2, 10);Требование к программе Матрицы прямоугольные, размеры передаются через параметрызапуска Тип элементов float, инициализация матриц случайными числами поформуле rand() / RAND_MAX Обернуть ВСЕ матрицы в cudaPitchedPtr, обращаться черезget_elem Реализовать возможность проверки правильности результата (понеобязательному последнему флагу в параметрах запуска) Вычисление эталонного произведения матриц на хосте исравнение с результатом CUDAДобавляем padding (набивку) На хосте заменяем cudaMalloc на cudaMallocPitch cudaMemcpy на cudaMemcpy2D На устройстве обращаемся к элементам матриц с учетом ихpitch Именно для этого обернули матрицы в cudaPitchedPtr –чтобы не плодить лишних переменных и хранить pitchвместе с указателямиcudaMallocPitch & cudaPitchedPtr Выделение памяти под матрицу с элементами типа Type илогическими размерами (в элементах) width x height:cudaPitchedPtr matrix =make_cudaPitchedPtr(0,0, width, height);cudaMallocPitch(&matrix.ptr, &matrix.pitch,width * sizeof( Type ), height); Если Type==int32_t, a width==1000, то будет выделено 4096* heightбайтов и в matrix.pitch запишется 4096cudaMemcpy2D & cudaPitchedPtr Пример пересылки c хоста на GPU:cudaPitchedPtr aHost, aDev;…// выделение памяти под матрицу на GPU и на хостеcudaMemcpy2D( aDev.ptr, aDev.pitch, aHost.ptr,aHost.pitch, aHost.pitch, aHost.ysize,cudaMemcpyHostToDevice); aHost – матрица на хосте, сплошная (без набивки), значит еёфактическая ширина(pitch) равна обычной ширине в байтах aDev – матрица на устройстве, память под неё выделена на при помощиcudaMallocPitch, как было описано вышеПравильная работа с памятью Каждая нить рассчитывает один элемент матрицы-результата Желательно, чтобы обращения в память нитей одного варпа попадали впромежуток памяти размером 128 байт адрес начала этого промежутка делится на 128 ( выровнен по128) В этом случае обращения варпа в память будут выполнены заодну транзакцию к кешу L1Правильная работа с памятью При вычислении скалярного произведения нить пробегаетстолбец матрицы B и строку матрицы AНужно ли транспонировать матрицу B?Нужно ли транспонировать матрицу A?Нужно ли транспонировать матрицу B?Нужно ли транспонировать матрицу A?Пусть ширина каждого блока кратна размеру варпа• т.е.

все нити одного варпа лежат в одной строке блокаи вычисляют соседние элементы в строке результатаМатрицы 10*10, варп 8 нитейa.ptr[idy * a.xsize + k]k = 0Матрицы 10*10, варп 8 нитейb.ptr[k * b.xsize + idx]k = 0Матрицы 10*10, варп 8 нитейa.ptr[idy * a.xsize + k]k = 1Матрицы 10*10, варп 8 нитейb.ptr[k * b.xsize + idx]k = 1Матрицы 10*10, варп 8 нитейa.ptr[idy * a.xsize + k]k = 2Матрицы 10*10, варп 8 нитейb.ptr[k * b.xsize + idx]k = 2Матрицы 10*10, варп 8 нитейc.ptr[idy * c.xsize + idx] = tmpНужно ли транспонировать матрицу B?Нужно ли транспонировать матрицу A? Если ширина блока кратна размеру варпа, то все нити варпаобращаются к одному элементу матрицы А и соседним элементам встроке матрицы B Транспонировать не нужноНужно ли транспонировать матрицу B?Нужно ли транспонировать матрицу A? Что если ширина блока не кратна размеру варпа? т.е.

нити одного варпа лежат не в одной строке блокаМатрицы 10*10, варп 8 нитей, блок 6x2a.ptr[idy * a.xsize + k]k = 0Матрицы 10*10, варп 8 нитей, блок 6x2b.ptr[k * b.xsize + idx]k = 0Нужно ли транспонировать матрицу B?Нужно ли транспонировать матрицу A? Если ширина блока не кратна размеру варпа, то все нитиварпа обращаются к соседним элементам столбцаматрицы А и соседним элементам в строке матрицы B Нужно транспонировать матрицу АТребования к программе Программа должна запускать ядра на устройствах сcudaDeviceProp::major == 2 Память под все матрицы выделять черезcudaMallocPitch Выводить время работы ядер Время работы считать через события Реализовать возможность задания высоты блока гридаТребования к программе Таким образом, параметры программы: M, N, K - размеры матриц Высота блоков 16 или 32? Проверять результат? На выходе: Время работы ядра без набивки Время работы ядра с набивкойРезультаты тестовИдея Пусть нужно перемножить две «полоски»шириной 32 элемента Пусть размер блока 32х32 Делаем два массива в общей памяти сразмерами 32х32 (наш управляемый кеш) В первый копируем блок матрицы А, вовторой – В (закешировали) Вычисляем блок матрицы C используяданные из общей памятиИдея Вычисляем блок матрицы C используяданные из общей памяти В итоге за время работы блока к каждомуэлементу в глобальной памяти мыобратимся всего один раз Без общей памяти – 32 разаИдея• Теперь пусть нужно перемножить матрицыпроизвольного размера Разрежем их на «полоски» Матрица C есть сумма матриц, получаемыхпри перемножении соответствующихполосокИдея• Матрица C есть сумма матриц, получаемыхпри перемножении соответствующихполосок Так же делаем 2 массива по 32х32 вобщей памяти, кешируем блоки полосок Число этапов – (A.width - 1) / 32 + 1 Накапливаем результатМатрицы произвольного размера Чтобы не возиться с выходом за границы матриц, следует явнодополнить матрицы нулями до размеров, кратных 32 new_width = ((width - 1)/32 + 1) new_height = ((height - 1)/32 + 1)Оптимизация В интернете легко можно найти примеры реализации перемноженияматриц с общей памятью, поэтому проведем некоторые оптимизации: Высота блока – 16 нитей Каждая нить вычисляет два элемента: a[i,j] и a[i+16, j] Получается ускорение 1.5x, по сравнению с обычным вариантом По аналогии можно сделать ядра, где каждая нить считает > 2элементов, которые будут еще эффективнееТребования к программе Общую память выделять динамически Дополнить матрицы нулями до размеров, кратных 32 Для наглядности ускорения добавить в графики из прошлой части: График скорости работы варианта с общей памятью при блоке32x32, нить вычисляет один элемент График скорости работы варианта с общей памятью при блоке32x16, нить вычисляет два элементаend.