Комментарии 14
А вот интересно: у меня сделаны самодельные нейронки и там в основе это самое матричное умножение (и для свёрток тоже) с использованием shared памяти на CUDA. Так вот, игрался я на 1070 ti и решил купить Tesla P100 потому как игры с памятью (банальная перестановка доступа к массиву) сильно влияют на быстродействие, а у Tesla P100 шина 4096 бит, и я полагал, что это сильно ускорит чтение/запись памяти (не shared). Купил. И ускорения практически не получил. И возникает вопрос, а где же эти 4096 бит тогда используют и где они дают существенный прирост производительности?
В моем понимании, для того чтобы воспользоваться преимуществом широкой шины надо что бы приложение могло за такт обеспечить требуемое количество запросов к памяти. Если доступ не сгруппированный то будет много отдельных транзакций которые забьют этот размер шины и пропускная способность упадёт. Или если SM недостаточно загружены, или много синхронизаций то тоже не получиться сформировать нужное количество запросов. Для более точного ответа нужно наверное начать с профилирования и определить узкие места.
Спасибо за статью, тема очень интересная.
Жалко, много фактических ошибок:
Если реализовывать это обычным способом для CPU, все выходит достаточно просто. Мы получаем три вложенных цикла, в которых перемножаем элементы и сохраняем результирующий элемент матрицы C:
void matmul(int M, int N, int K,
const float* A, const float* const B, float* C) {
for (int m = 0; m < M; m++)
for (int n = 0; n < N; n++)
for (int k = 0; k < K; k++)
C[m * M + n] += A[m * M + k] * B[k * K + n];
}Ладно матрица С не зануляется, так еще и индексация полностью неверна. Стоило запускать тесты на неквадратных матрицах, вероятно был бы сегфолт
Одному потоку еще доступна private memory или регистры в терминологии CUDA. Также они разделяют пространство памяти с кешем первого уровня: можно регулировать, сколько отойдет на регистры, сколько на кеш первого уровня.
Нет, регистровый файл стоит отдельно, его размер фиксирован. Настраивается только доля Shared Memory vs L1-cache.
Размер файла регистров ограничен для каждого стриминг-процессора. Например для A100 — 16384х32bit.
Это не для SM, а на partition (который вы называете streaming processor). Вот из спеки Ampere, отличие как раз в 4 раза:
The register file size is 64K 32-bit registers per SM.
И вот так выглядит код простого подхода для вычисления одного элемента матрицы C:
__kernel void naive_gemm(int m, int n, int k,
const __global float* a,
const __global float* b,
__global float* c) {
const int global_row = get_global_id(0);
const int global_col = get_global_id(1);
float acc = 0.0f;
for (int i = 0; i < k; ++i) {
acc += a[i * m + global_row] * b[global_col * i + k];
}
c[global_col * m + global_row] = acc;
}Опять беда с индексами, еще обязательно надо добавить гвард: if (global_row>=m || global_col>=n) return
По итогу дальнейшие сравнения бессмыслены -- умножается совсем не то, и паттерн доступа к памяти уже другой.
Хотелось бы еще увидеть методологию измерения перфа, на GPU это очень чувствительная тема. Без неё результаты остаются невоспроизводимыми. И результаты обычно приводят в GFLOP/s
Каждый поток во время выполнения будет дожидаться завершения всех операций c локальной памятью:
const int tiled_row = TS * t + row;
const int tiled_col = TS * t + col;
a_sub[col][row] = A[tiled_col * M + global_row];
b_sub[col][row] = B[global_col * K + tiled_row];
// синхронизируем что быть уверенным, что tile загружен
barrier(CLK_LOCAL_MEM_FENCE);for (int k=0; k < TS; k++) {
acc += a_sub[k][row] * b_sub[col][k];
}
// синхронизируем перед следующим tile
barrier(CLK_LOCAL_MEM_FENCE);Здесь сильно запутанно: почему и зачем загрузка в shared memory происходит с транпонированием? В каком вообще формате входные матрицы (col- или row- major)?
Например, для доступа к 32-байтовым элементам одна транзакция будет размером 128 байт
32-битным элементам
Важно помнить о нюансах точности таких специальных ядер: они могут выполнять вычисления в половинной точности (FP16), что подходит не для всех задач.
Стоило упомянуть, что TensorOps умеют умножать матрицы и в TF32, что ускоряет умножение с достаточной для многих задач точностью
Все предложенные статьи в конце действительно очень интересные, спасибо за список
Спасибо за развернутый комментарий и найденные неточности. Я действительно допустил ошибки в индексации для простых реализаций. В более сложных все правильно. Также я специально опустил проверку размерностей для упрощения и уменьшения количества кода, т.е. примеры работают только с квадратными матрицами в формате column-major. С размером регистрового файла и локальной памяти тоже было написано неверно. И под стриминг процессором я понимаю Streaming Processor(SP) а не Streaming Multiprocessor(SM), т.е. четвертую часть как вы и написали. Методика измерения производительности использовалась самая простая - просто N раз запустил и усреднил, конечно для серьезных измерений такое не подходит так есть много нюансов, включая плавающую частоту. Но думаю для демонстрации такой подход вполне приемлем учитывая разницу производительности в несколько порядков.
Еще раз спасибо что обратили внимание на ошибки, я внес изменения в статью.
Еще в статье не описаны многие переменные. Например, что такое TSK, K, RTSN, RTSM?
Хочу попробовать реализовать через регистры, и честно говоря, уже запутался в индексах и что куда суммируется. :)
Так же нужно бы описать сетку, на которой запускаются потоки. Вот при использовании разделяемой памяти сетка равна размер выходной матрицы деленный на размер тайла. А при регистрах? Делить на размер тайла умноженный на размер масштабного коэффициента?
Так же не понятно, вот эта игра с регистрами, она выполняется после того, как потоки положили данные тайла в разделяемую память? То есть, это просто замена цикла сложений результатов умножения по строке и столбцу тайла?
Посмотрите пожалуйста исходный код проекта https://gitverse.ru/kolkir/clgemm думаю вы там найдете ответы. Вот описания констант:
#define TS 32
#define TSM 128 // The tile-size in dimension M
#define TSN 128 // The tile-size in dimension N
#define TSK 16 // The tile-size in dimension K
#define WPTM 8 // The amount of work-per-thread in dimension M
#define WPTN 8 // The amount of work-per-thread in dimension N
#define RTSM (TSM / WPTM) // The reduced tile-size in dimension M (== number of threads)
#define RTSN (TSN / WPTN) // The reduced tile-size in dimension N (== number of threads)
...В файле optim_gemm.cc находится реализация ядра со всеми оптимизациями, в тут запуск ядер.
Спасибо, посмотрю.
Там же по сути всего лишь такая конструкция:
static const size_t m_size=512;
CTensor<type_t> cTensorA(1,m_size,m_size);
CTensor<type_t> cTensorB(1,m_size,m_size);
CTensor<type_t> cTensorC(1,m_size,m_size);
//задаём матрицы
for(size_t y=0;y<m_size;y++)
{
for(size_t x=0;x<m_size;x++)
{
cTensorA.SetElement(0,y,x,(x+y*m_size)%5);
cTensorB.SetElement(0,y,x,(x*m_size+y)%5);
}
}
//очищаем тензор результата
cTensorC.Zero();
//делаем умножение
for(size_t y=0;y<cTensorC.GetSizeY();y++)
{
for(size_t x=0;x<cTensorC.GetSizeX();x++)
{
type_t line[m_size];//массив размера cTensorA.GetSizeX()=cTensorB.GetSizeY()
//заполняем линию
for(size_t k=0;k<m_size;k++) line[k]=cTensorB.GetElement(0,x,k);
//берём элемент a
type_t a=cTensorA.GetElement(0,y,x);
//перемножаем с элементами b и складываем в выходном тензоре
for(size_t k=0;k<m_size;k++)
{
type_t v=cTensorC.GetElement(0,y,k);
v+=a*line[k];
cTensorC.SetElement(0,y,k,v);
}
}
}При этом x и y внешнего цикла это как раз координаты потока CUDA внутри выходного тензора.
Наврал. :)
Вот так должно быть.
static const size_t ma_size_x=512;
static const size_t ma_size_y=384;
static const size_t mb_size_x=128;
static const size_t mb_size_y=ma_size_x;
static const size_t mc_size_x=mb_size_x;
static const size_t mc_size_y=ma_size_y;
CTensor<type_t> cTensorA(1,ma_size_y,ma_size_x);
CTensor<type_t> cTensorB(1,mb_size_y,mb_size_x);
CTensor<type_t> cTensorC(1,mc_size_y,mc_size_x);
//заполняем матрицы
for(size_t y=0;y<ma_size_y;y++)
{
for(size_t x=0;x<ma_size_x;x++) cTensorA.SetElement(0,y,x,(x+y*ma_size_x)%5);
}
for(size_t y=0;y<mb_size_y;y++)
{
for(size_t x=0;x<mb_size_x;x++) cTensorB.SetElement(0,y,x,(x*mb_size_y+y)%5);
}
//очищаем тензор результата
cTensorC.Zero();
//делаем умножение
type_t line[mc_size_x];//размер выходного тензора по x
for(size_t ax=0;ax<ma_size_x;ax++)
{
for(size_t k=0;k<mc_size_x;k++) line[k]=cTensorB.GetElement(0,ax,k);
for(size_t ay=0;ay<ma_size_y;ay++)
{
//берём элемент a
type_t a=cTensorA.GetElement(0,ay,ax);
for(size_t k=0;k<mc_size_x;k++)
{
type_t v=cTensorC.GetElement(0,ay,k);
v+=a*line[k];
cTensorC.SetElement(0,ay,k,v);
}
}
}Очень странно. Вот если взять файл optim_gemm.cc и в файле запуска ядер глянуть параметры запуска
lrange = cl::NDRange{TSM / WPTM, TSN / WPTN};
grange = cl::NDRange{matrix_data.m / WPTM, matrix_data.n / WPTN};То получится что-то странное для grange. Потому что в ядре написано так:
const int tidm = get_local_id(0); // Local row ID (max: TSM/WPTM == RTSM)
const int tidn = get_local_id(1); // Local col ID (max: TSN/WPTN == RTSN)
const int offset_m = TSM * get_group_id(0); // Work-group offset
const int offset_n = TSN * get_group_id(1); // Work-group offsetТо есть, offset_m и offset_n дают смещение на сетке TSM и TSN. Вероятно, там всё-таки grange должно быть получено делением размеров матриц на эти самые TSM и TSN, соответственно. А так - WPTM=8. То есть, количество блоков при делении размера на 8 будет много больше, чем при делении на TSM=128. И быстродействие сильно просядет для матрицы 8192x8192.
А если адаптировать всё это для CUDA и попутно поменять все эти TSM и TSN в которых я путаюсь на удобные (мне) названия:
static const size_t TENSOR_OPERATION_TILE_SIZE_X=128;///<размер блока операций с тензорами X
static const size_t TENSOR_OPERATION_TILE_SIZE_Y=128;///<размер блока операций с тензорами Y
static const size_t TENSOR_OPERATION_TILE_SIZE_K=16;///<размер блока операций с тензорами K
static const size_t WORK_PER_TILE_SIZE_X=8;///<сколько элементов обрабатывает поток X
static const size_t WORK_PER_TILE_SIZE_Y=8;///<сколько элементов обрабатывает поток Y
static const size_t WORK_PER_TILE_X=TENSOR_OPERATION_TILE_SIZE_X/WORK_PER_TILE_SIZE_X;///<количество обрабатываемых блоков по WORK_PER_TILE_SIZE_X элементов
static const size_t WORK_PER_TILE_Y=TENSOR_OPERATION_TILE_SIZE_Y/WORK_PER_TILE_SIZE_Y;///<количество обрабатываемых блоков по WORK_PER_TILE_SIZE_Y элементов
//----------------------------------------------------------------------------------------------------
//функция CUDA для умножения тензоров
//----------------------------------------------------------------------------------------------------
template<class type_t,class kernel_output_t,class kernel_left_t,class kernel_right_t>
__global__ void CUDATensorMulTensorFunction(kernel_output_t tensor_output,kernel_left_t tensor_left,kernel_right_t tensor_right)
{
const size_t tidm=threadIdx.y;// Local row ID (max: TSM/WPTM == RTSM)
const size_t tidn=threadIdx.x;// Local col ID (max: TSN/WPTN == RTSN)
const size_t offset_m=TENSOR_OPERATION_TILE_SIZE_Y*blockIdx.y;// Work-group offset
const size_t offset_n=TENSOR_OPERATION_TILE_SIZE_X*blockIdx.x;// Work-group offset
const size_t z=blockIdx.z;
// Local memory to fit a tile of A and B
__shared__ type_t a_sub[TENSOR_OPERATION_TILE_SIZE_K][TENSOR_OPERATION_TILE_SIZE_Y];
__shared__ type_t b_sub[TENSOR_OPERATION_TILE_SIZE_X][TENSOR_OPERATION_TILE_SIZE_K+2];
// Allocate register space
type_t a_reg;
type_t b_reg[WORK_PER_TILE_SIZE_X];
type_t acc[WORK_PER_TILE_SIZE_Y][WORK_PER_TILE_SIZE_X];
// Initialise the accumulation registers
#pragma unroll
for(size_t wm=0;wm<WORK_PER_TILE_SIZE_Y;wm++)
{
#pragma unroll
for (size_t wn=0;wn<WORK_PER_TILE_SIZE_X;wn++) acc[wm][wn]=0;
}
// Loop over all tiles
size_t numTiles=tensor_left.Size_X/TENSOR_OPERATION_TILE_SIZE_K;
if (tensor_left.Size_X%(TENSOR_OPERATION_TILE_SIZE_K)) numTiles++;
size_t LPTA=(TENSOR_OPERATION_TILE_SIZE_K*WORK_PER_TILE_SIZE_Y*WORK_PER_TILE_SIZE_X)/(TENSOR_OPERATION_TILE_SIZE_X);// The amount of loads-per-thread for A
size_t LPTB=(TENSOR_OPERATION_TILE_SIZE_K*WORK_PER_TILE_SIZE_Y*WORK_PER_TILE_SIZE_X)/(TENSOR_OPERATION_TILE_SIZE_Y);// The amount of loads-per-thread for B
size_t t=0;
do
{
// Load one tile of A and B into local memory
#pragma unroll
for(size_t la=0;la<LPTA;la++)
{
size_t tid=tidn*WORK_PER_TILE_Y+tidm;
size_t id=la*WORK_PER_TILE_X*WORK_PER_TILE_Y+tid;
size_t col=id/TENSOR_OPERATION_TILE_SIZE_Y;
size_t row=id-col*TENSOR_OPERATION_TILE_SIZE_Y;
size_t tiled_index=TENSOR_OPERATION_TILE_SIZE_K*t+col;
a_sub[col][row]=tensor_left.GetElement(z,offset_m+row,tiled_index);//A[tiled_index*M+offset_m+row];
b_sub[row][col]=tensor_right.GetElement(z,offset_n+row,tiled_index);//B[tiled_index*N+offset_n+row];
}
// Synchronise to make sure the tile is loaded
__syncthreads();
// Loop over the values of a single tile
for(size_t k=0;k<TENSOR_OPERATION_TILE_SIZE_K;k++)
{
// Cache the values of b_sub in registers
#pragma unroll
for(size_t wn=0;wn<WORK_PER_TILE_SIZE_X;wn++)
{
size_t col=tidn+wn*WORK_PER_TILE_X;
b_reg[wn]=b_sub[col][k];
}
// Perform the computation
#pragma unroll
for(size_t wm=0;wm<WORK_PER_TILE_SIZE_Y;wm++)
{
size_t row=tidm+wm*WORK_PER_TILE_Y;
a_reg=a_sub[k][row];
#pragma unroll
for(size_t wn=0;wn<WORK_PER_TILE_SIZE_X;wn++) acc[wm][wn]+=a_reg*b_reg[wn];
}
}
// Synchronise before loading the next tile
__syncthreads();
// Next tile
t++;
}
while(t<numTiles);
// Store the final results in C
#pragma unroll
for(size_t wm=0;wm<WORK_PER_TILE_SIZE_Y;wm++)
{
size_t global_row=offset_m+tidm+wm*WORK_PER_TILE_Y;
#pragma unroll
for(size_t wn=0;wn<WORK_PER_TILE_SIZE_X;wn++)
{
size_t global_col=offset_n+tidn+wn*WORK_PER_TILE_X;
tensor_output.SetElement(z,global_row,global_col,acc[wm][wn]);
}
}
}
//----------------------------------------------------------------------------------------------------
//умножить тензоры
//----------------------------------------------------------------------------------------------------
template<class type_t,class kernel_output_t,class kernel_left_t,class kernel_right_t>
__host__ void MulAbstract(CTensor<type_t> &cTensor_Output,kernel_output_t &sTensorKernel_Output,const CTensor<type_t> &cTensor_Left,kernel_left_t &sTensorKernel_Left,const CTensor<type_t> &cTensor_Right,kernel_right_t &sTensorKernel_Right)
{
if (sTensorKernel_Left.Size_X!=sTensorKernel_Right.Size_Y || sTensorKernel_Left.Size_Z!=sTensorKernel_Right.Size_Z ||
sTensorKernel_Output.Size_Y!=sTensorKernel_Left.Size_Y || sTensorKernel_Output.Size_X!=sTensorKernel_Right.Size_X ||
sTensorKernel_Output.Size_Z!=sTensorKernel_Right.Size_Z)
{
throw "CTensor::MulAbstract: Размерности тензоров не совпадают!";
}
//копируем данные с процессора
cTensor_Left.CopyToDevice();
cTensor_Right.CopyToDevice();
dim3 thread(WORK_PER_TILE_X,WORK_PER_TILE_Y);
size_t block_x=sTensorKernel_Right.Size_X/(TENSOR_OPERATION_TILE_SIZE_X);
size_t block_y=sTensorKernel_Left.Size_Y/(TENSOR_OPERATION_TILE_SIZE_Y);
size_t block_z=sTensorKernel_Output.Size_Z;
dim3 blocks(block_x,block_y,block_z);
if (blocks.x==0) blocks.x=1;
if (blocks.y==0) blocks.y=1;
if (blocks.z==0) blocks.z=1;
CUDATensorMulTensorFunction<type_t,kernel_output_t,kernel_left_t,kernel_right_t><<<blocks,thread>>>(sTensorKernel_Output,sTensorKernel_Left,sTensorKernel_Right);
HANDLE_ERROR(cudaGetLastError());
HANDLE_ERROR(cudaDeviceSynchronize());
cTensor_Output.SetDeviceOnChange();
}Ядро тензора для CUDA:
//****************************************************************************************************
///!структура ядра тензора
//****************************************************************************************************
template<class type_t>
struct STensorKernel
{
size_t Size_X;///<размер по x
size_t Size_Y;///<размер по y
size_t Size_Z;///<размер по z
size_t StrideX;///<строка по X
size_t StrideZ;///<размер блока по Z
type_t *TensorData_Ptr;///<указатель на данные тензора на стороне GPU
__host__ __device__ STensorKernel(void)///<конструктор
{
}
__host__ __device__ STensorKernel(const CTensor<type_t> &cTensor)///<конструктор
{
Set(cTensor);
}
__forceinline__ __host__ __device__ type_t* GetTensorDataPtr(size_t z)///<получить указатель на элементы с глубиной z
{
return(&TensorData_Ptr[z*StrideZ]);
}
__forceinline__ __host__ __device__ type_t GetElement(size_t z,size_t y,size_t x)
{
if (x>=Size_X || y>=Size_Y || z>=Size_Z) return(0);
return(TensorData_Ptr[z*StrideZ+y*StrideX+x]);
}
__forceinline__ __host__ __device__ void SetElement(size_t z,size_t y,size_t x,type_t value)
{
if (x>=Size_X || y>=Size_Y || z>=Size_Z) return;
TensorData_Ptr[z*StrideZ+y*StrideX+x]=value;
}
__forceinline__ __host__ __device__ size_t GetSizeX(void) const
{
return(Size_X);
}
__forceinline__ __host__ __device__ size_t GetSizeY(void) const
{
return(Size_Y);
}
__forceinline__ __host__ __device__ size_t GetSizeZ(void) const
{
return(Size_Z);
}
__host__ __device__ void Reset(void)
{
Size_X=0;
Size_Y=0;
Size_Z=0;
StrideX=0;
StrideZ=0;
TensorData_Ptr=NULL;
}
__host__ __device__ void Set(const CTensor<type_t> &cTensor)
{
Size_Z=cTensor.Size_Z;
Size_Y=cTensor.Size_Y;
Size_X=cTensor.Size_X;
StrideX=cTensor.Size_X;
StrideZ=cTensor.Size_X*cTensor.Size_Y;
TensorData_Ptr=cTensor.DeviceItem.get();
}
};И мой текущий вариант умножения:
static const size_t TENSOR_OPERATION_BLOCK_SIZE_SCALE=2;///<размер множителя блока операций с тензорами
static const size_t TENSOR_OPERATION_BLOCK_SIZE=16;///<размер блока операций с тензорами
static const size_t TENSOR_OPERATION_TILE_SIZE=TENSOR_OPERATION_BLOCK_SIZE*TENSOR_OPERATION_BLOCK_SIZE_SCALE;
//----------------------------------------------------------------------------------------------------
//функция CUDA для умножения тензоров
//----------------------------------------------------------------------------------------------------
template<class type_t,class kernel_output_t,class kernel_left_t,class kernel_right_t>
__global__ void CUDATensorMulTensorFunction(kernel_output_t tensor_output,kernel_left_t tensor_left,kernel_right_t tensor_right)
{
//блок TENSOR_OPERATION_BLOCK_SIZE x TENSOR_OPERATION_BLOCK_SIZE в выходном тензоре
size_t block_x=blockIdx.x;
size_t block_y=blockIdx.y;
//координаты элементов блока в выходном тензоре
size_t out_in_block_x=threadIdx.x*TENSOR_OPERATION_BLOCK_SIZE_SCALE;
size_t out_in_block_y=threadIdx.y*TENSOR_OPERATION_BLOCK_SIZE_SCALE;
//координаты блока в выходном тензоре
size_t out_block_x=TENSOR_OPERATION_TILE_SIZE*block_x;
size_t out_block_y=TENSOR_OPERATION_TILE_SIZE*block_y;
//глобальные координаты в выходном тензоре
size_t out_x=out_in_block_x+out_block_x;
size_t out_y=out_in_block_y+out_block_y;
size_t out_z=blockIdx.z;
//получаем подматрицу выходной матрицы
type_t Cvalue[TENSOR_OPERATION_BLOCK_SIZE_SCALE][TENSOR_OPERATION_BLOCK_SIZE_SCALE];
for(size_t ky=0;ky<TENSOR_OPERATION_BLOCK_SIZE_SCALE;ky++)
{
for(size_t kx=0;kx<TENSOR_OPERATION_BLOCK_SIZE_SCALE;kx++)
{
Cvalue[ky][kx]=0;
}
}
//считаем, сколькно нужно проходов блоком по X
size_t m_max=tensor_left.Size_X/(TENSOR_OPERATION_TILE_SIZE);
if (tensor_left.Size_X%(TENSOR_OPERATION_TILE_SIZE)) m_max++;
__shared__ type_t As[TENSOR_OPERATION_TILE_SIZE][TENSOR_OPERATION_TILE_SIZE];
__shared__ type_t Bs[TENSOR_OPERATION_TILE_SIZE][TENSOR_OPERATION_TILE_SIZE];
for(size_t m=0;m<m_max;m++)
{
size_t offset=m*TENSOR_OPERATION_TILE_SIZE;
size_t py_left=out_in_block_y+out_block_y;
size_t py_right=out_in_block_y+offset;
size_t py=out_in_block_y;
for(size_t ky=0;ky<TENSOR_OPERATION_BLOCK_SIZE_SCALE;ky++,py_left++,py_right++,py++)
{
size_t px_left=out_in_block_x+offset;
size_t px_right=out_in_block_x+out_block_x;
size_t px=out_in_block_x;
for(size_t kx=0;kx<TENSOR_OPERATION_BLOCK_SIZE_SCALE;kx++,px_left++,px_right++,px++)
{
As[py][px]=tensor_left.GetElement(out_z,py_left,px_left);
Bs[py][px]=tensor_right.GetElement(out_z,py_right,px_right);
}
}
__syncthreads();
//выполняем умножение только для тех элементов, которые не выходят за размер выходного тензора
if ((out_y<tensor_output.Size_Y) &&
(out_x<tensor_output.Size_X))
{
for(size_t ky=0;ky<TENSOR_OPERATION_BLOCK_SIZE_SCALE;ky++)
{
for(size_t kx=0;kx<TENSOR_OPERATION_BLOCK_SIZE_SCALE;kx++)
{
type_t &cv=Cvalue[ky][kx];
size_t ay=out_in_block_y+ky;
size_t bx=out_in_block_x+kx;
for(size_t e=0;e<TENSOR_OPERATION_TILE_SIZE;e++) cv+=As[ay][e]*Bs[e][bx];
}
}
}
__syncthreads();
}
size_t py=out_y;
for(size_t ky=0;ky<TENSOR_OPERATION_BLOCK_SIZE_SCALE;ky++,py++)
{
size_t px=out_x;
for(size_t kx=0;kx<TENSOR_OPERATION_BLOCK_SIZE_SCALE;kx++,px++)
{
tensor_output.SetElement(out_z,py,px,Cvalue[ky][kx]);
}
}
__syncthreads();
}
//----------------------------------------------------------------------------------------------------
//умножить тензоры
//----------------------------------------------------------------------------------------------------
template<class type_t,class kernel_output_t,class kernel_left_t,class kernel_right_t>
__host__ void MulAbstract(CTensor<type_t> &cTensor_Output,kernel_output_t &sTensorKernel_Output,const CTensor<type_t> &cTensor_Left,kernel_left_t &sTensorKernel_Left,const CTensor<type_t> &cTensor_Right,kernel_right_t &sTensorKernel_Right)
{
if (sTensorKernel_Left.Size_X!=sTensorKernel_Right.Size_Y || sTensorKernel_Left.Size_Z!=sTensorKernel_Right.Size_Z ||
sTensorKernel_Output.Size_Y!=sTensorKernel_Left.Size_Y || sTensorKernel_Output.Size_X!=sTensorKernel_Right.Size_X ||
sTensorKernel_Output.Size_Z!=sTensorKernel_Right.Size_Z)
{
throw "CTensor::MulAbstract: Размерности тензоров не совпадают!";
}
//копируем данные с процессора
cTensor_Left.CopyToDevice();
cTensor_Right.CopyToDevice();
dim3 thread(TENSOR_OPERATION_TILE_SIZE,TENSOR_OPERATION_TILE_SIZE);
size_t block_x=sTensorKernel_Output.Size_X/thread.x;
if (sTensorKernel_Output.Size_X%thread.x) block_x++;
size_t block_y=sTensorKernel_Output.Size_Y/thread.y;
if (sTensorKernel_Output.Size_Y%thread.y) block_y++;
size_t block_z=sTensorKernel_Output.Size_Z;
dim3 blocks(block_x,block_y,block_z);
if (blocks.x==0) blocks.x=1;
if (blocks.y==0) blocks.y=1;
if (blocks.z==0) blocks.z=1;
dim3 thread_basic(TENSOR_OPERATION_BLOCK_SIZE,TENSOR_OPERATION_BLOCK_SIZE);
CUDATensorMulTensorFunction<type_t,kernel_output_t,kernel_left_t,kernel_right_t><<<blocks,thread_basic>>>(sTensorKernel_Output,sTensorKernel_Left,sTensorKernel_Right);
HANDLE_ERROR(cudaGetLastError());
HANDLE_ERROR(cudaDeviceSynchronize());
cTensor_Output.SetDeviceOnChange();
}Так вот, такая вот штука для матрицы 8192x8192 работает на I5-3570 так:
Тип 2 - это мой вариант. Тип 3 - это умножение на регистрах.
Вроде как 500 мс выиграны. Но я не уверен в правильности работы умножений на регистрах во всех размерностях матриц (неквадратные я не проверял ещё).
Завтра проверю на p106-100 (аналог 1060 с урезанной pci-e x1). Там вроде как результаты были противоположные. Вариант с регистрами дал 3700 мс, а мой вариант дал 933 мс. Проверю завтра. Может, что-то не так запускал.
Нет, всё-таки, похоже
grange = cl::NDRange{matrix_data.m / WPTM, matrix_data.n / WPTN};
Но тогда быстродействие просто в ноль уходит. То, что мой вариант делает за 0.11 мс, этот делает аж за 3.82. И это всего лишь небольшая матрица. С большой я ждал результата около минуты-двух. :O
А ещё вот тут b_sub[row][col]=tensor_right.GetElement(z,offset_n+row,tiled_index);//B[tiled_index*N+offset_n+row];
должно быть
b_sub[row][col]=tensor_right.GetElement(z,tiled_index,offset_n+row);//B[tiled_index*N+offset_n+row];
Смущает только tiled_index*N у B. Ведь матрица B задана как N по X и K по Y. Следовательно, точка (x,y) имеет смещение x*K+y (тут хранение массива как [y][x]); А в индексах x умножается на N, а не на K.
Это всё вообще работало у автора репозитория?
Вы в предыдущем посте правильно обратили внимание на момент с расчетом элементов в глобальном NDRange, дело в том, что в OpenCL global_range должен задавать общее количество элементов сетки, и local_range будет использоваться для определения количества групп. Т.е. не так как в CUDA когда задается количество блоков. Я действительно не упомянул про это в статье, сейчас добавил. Спасибо что обратили внимание. И у меня используются матрицы в column-major формате это немного запутывает индексацию.
То есть, число блоков будет отношением global к local?
Да, 
, вот ссылка на спецификацию.
Учимся разрабатывать для GPU на примере операции GEMM