Ядро CUDA 2×2 ZgemmBatched: в 5 раз быстрее, чем CuBLAS. Это может идти быстрее?
Я видел низкую производительность API CuBLAS при использовании небольших матриц, особенно пакетного общего умножения матриц. Моя цель — написать быстрое ядро для вычисления пакетных сложных двойных матричных умножений для размера матрицы 2×2. Ядро должно быть в состоянии извлечь эти матрицы из больших квадратных матриц, то есть умножения блочных матриц.
Ядро, которое я написал, рассчитывает 8 двойных комплексных матричных умножений 2×2 на блок (32 потока). Я могу достичь примерно 1,65 двойных GFlops (анализ NSight) на GTX 970 (вычислительная возможность 5.2), который должен выполнить 109 двойных GFlops. Требуется около 300 мс, чтобы рассчитать 9,1 млн. Умножений матриц.
*a, *b а также *c являются указателями на массив из N раз nxn матриц. Эти матрицы должны быть в квадрате.
lda, ldb, ldc являются ведущими размерами всех входных матриц. Используя эти параметры, любой квадратный матричный массив можно использовать для извлечения матриц 2х2.
Вопрос: Как мне повысить производительность этого ядра?
#define THREADMULTIPLIER 8
__global__ void bulkMatrixMul(const cuDoubleComplex *a, int lda, const cuDoubleComplex *b, int ldb, cuDoubleComplex *c, int ldc){
int pos = (THREADMULTIPLIER * blockIdx.x + threadIdx.z);
int ia = pos * lda * lda + threadIdx.x;
int ib = (pos * ldb + threadIdx.y) * ldb;
c[(pos * ldc + threadIdx.y) * ldc + threadIdx.x] = cuCfma(a[ia], b[ib], cuCmul(a[ia + LD], b[ib + 1]));
};void bulkBatchedMatrixMul(complex<double> *a, int lda, complex<double> *b, int ldb, complex<double> *c, int ldc, int batchSize, cudaStream_t *stream){
dim3 threads(2, 2, THREADMULTIPLIER);
bulkMatrixMul << <batchSize / THREADMULTIPLIER, threads, 0, *stream >> >((cuDoubleComplex*)a, lda, (cuDoubleComplex*)b, ldb, (cuDoubleComplex*)c , ldc);
if (cudaSuccess != cudaGetLastError())
printf("Error!\n");
}
Я сделал ядро в 10 раз быстрее, проанализировав код Кристиана Сарофина:
- это быстрее, если вы напишите сумму в целом вместо использования + = и — =.
- объявив выход, особенно если вы разделите реальный & мнимые расчеты на выходе увеличивают скорость.
Ядро ниже представляет собой оптимизированное ядро, рассчитывающее 16 умножений матриц на блок. Он может вычислить 9,1 миллиона умножений матриц за 31 мс: я использую batchSize 10000 и вызываю это ядро 70 раз на поток в 13 потоках. Согласно NSight он достигает около 15 двойных GFlops. GTX970 был использован.
__global__ void bulkMatrixMul(const cuDoubleComplex *a, int lda, const cuDoubleComplex *b, int ldb, cuDoubleComplex *c, int ldc){
int pos = (blockDim.z * blockIdx.x + threadIdx.z);
int ia = pos * lda * lda + threadIdx.x;
int ib = (pos * ldb + threadIdx.y) * ldb;
cuDoubleComplex cR;
cR.x = a[ia].x * b[ib].x - a[ia].y * b[ib].y - a[ia + lda].y * b[ib + 1].y + a[ia + lda].x * b[ib + 1].x;
cR.y = a[ia].x * b[ib].y + a[ia].y * b[ib].x + a[ia + lda].x * b[ib + 1].y + a[ia + lda].y * b[ib + 1].x;
c[(pos * ldc + threadIdx.y) * ldc + threadIdx.x] = cR;
};void bulkBatchedMatrixMul(complex<double> *a, int lda, complex<double> *b, int ldb, complex<double> *c, int ldc, int batchSize, cudaStream_t *stream){
dim3 threads(2, 2, 16);
bulkMatrixMul <<< batchSize / 16, threads, 0, *stream >>>((cuDoubleComplex*)a, lda, (cuDoubleComplex*)b, ldb, (cuDoubleComplex*)c , ldc);
if (cudaSuccess != cudaGetLastError())
printf("Error!\n");
}
Решение
Я бы назначил поток для каждого умножения матриц и выгрузил бы все в регистры. Не беспокойтесь о занятости, так как я подозреваю, что это будет связано с памятью. Попав в регистры, каждый поток будет выполнять матричное умножение и сохранять его в c. Ниже приведена демонстрация этого.
#include <stdio.h>
#include <iostream>
#include <cuComplex.h>
#include <cuda.h>
#include <cuda_runtime.h>
#define N_MATRIX 10000
using namespace std;
__device__ void D_C_Mult_Add(cuDoubleComplex &a, cuDoubleComplex &b, cuDoubleComplex &c){
c.x+=a.x*b.x;
c.y+=a.x*b.y;
c.x-=a.y*b.y;
c.y+=a.y*b.x;
}
__global__ void multMatrix(cuDoubleComplex* a, cuDoubleComplex* b, cuDoubleComplex* c) {
int tid = blockIdx.x*blockDim.x + threadIdx.x;
int stride=gridDim.x*blockDim.x;
if(tid>=N_MATRIX) return;
cuDoubleComplex a_sub[4], b_sub[4], c_sub[4];
for(int i=0; i<4; i++){
a_sub[i]=a[tid+stride*i];
b_sub[i]=a[tid+stride*i];
c_sub[i].x=0.0;
c_sub[i].y=0.0;
}
D_C_Mult_Add(a_sub[0], b_sub[0], c_sub[0]);
D_C_Mult_Add(a_sub[0], b_sub[1], c_sub[1]);
D_C_Mult_Add(a_sub[2], b_sub[0], c_sub[2]);
D_C_Mult_Add(a_sub[2], b_sub[1], c_sub[3]);
D_C_Mult_Add(a_sub[1], b_sub[2], c_sub[0]);
D_C_Mult_Add(a_sub[1], b_sub[3], c_sub[1]);
D_C_Mult_Add(a_sub[3], b_sub[2], c_sub[2]);
D_C_Mult_Add(a_sub[3], b_sub[3], c_sub[3]);
for(int i=0; i<4; i++)
c[tid+stride*i]=c_sub[i];
}
int main(int argc, char **argv) {cuDoubleComplex *a, *b, *c;
cudaMalloc(&a, sizeof(cuDoubleComplex)*N_MATRIX*4);
cudaMalloc(&b, sizeof(cuDoubleComplex)*N_MATRIX*4);
cudaMalloc(&c, sizeof(cuDoubleComplex)*N_MATRIX*4);
multMatrix<<<N_MATRIX/128+1, 128>>>(a, b, c);
}
Другие решения
Я настроил свое ядро, используя код Кристиана Сарофина. Его ядру нужно 60 мс, чтобы рассчитать 9,1 млн. Умножений матриц на GTX 970.
Мое настроенное ядро может сделать это за 31 мс. Спасибо, Кристиан!
__global__ void bulkMatrixMul(const cuDoubleComplex *a, int lda, const cuDoubleComplex *b, int ldb, cuDoubleComplex *c, int ldc){
int pos = (blockDim.z * blockIdx.x + threadIdx.z);
int ia = pos * lda * lda + threadIdx.x;
int ib = (pos * ldb + threadIdx.y) * ldb;
cuDoubleComplex cR;
cR.x = a[ia].x * b[ib].x - a[ia].y * b[ib].y - a[ia + lda].y * b[ib + 1].y + a[ia + lda].x * b[ib + 1].x;
cR.y = a[ia].x * b[ib].y + a[ia].y * b[ib].x + a[ia + lda].x * b[ib + 1].y + a[ia + lda].y * b[ib + 1].x;
c[(pos * ldc + threadIdx.y) * ldc + threadIdx.x] = cR;
};void bulkBatchedMatrixMul(complex<double> *a, int lda, complex<double> *b, int ldb, complex<double> *c, int ldc, int batchSize, cudaStream_t *stream){
dim3 threads(2, 2, 16);
bulkMatrixMul << <batchSize / 16, threads, 0, *stream >> >((cuDoubleComplex*)a, lda, (cuDoubleComplex*)b, ldb, (cuDoubleComplex*)c , ldc);
if (cudaSuccess != cudaGetLastError())
printf("Error!\n");
}