Ядро 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 раз быстрее, проанализировав код Кристиана Сарофина:

  1. это быстрее, если вы напишите сумму в целом вместо использования + = и — =.
  2. объявив выход, особенно если вы разделите реальный & мнимые расчеты на выходе увеличивают скорость.

Ядро ниже представляет собой оптимизированное ядро, рассчитывающее 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");
}

0

Решение

Я бы назначил поток для каждого умножения матриц и выгрузил бы все в регистры. Не беспокойтесь о занятости, так как я подозреваю, что это будет связано с памятью. Попав в регистры, каждый поток будет выполнять матричное умножение и сохранять его в 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);
}
1

Другие решения

Я настроил свое ядро, используя код Кристиана Сарофина. Его ядру нужно 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");
}
-2