Не квадратное матричное умножение в CUDA

Question

Код, который я использую для умножения матриц в CUDA, позволяет мне умножать как квадратные, так и неквадратные матрицы, однако ширина и высота ДОЛЖНЫ быть кратны размеру блока.

Так, например, я могу умножить [3] [6] * [6] [3] (используя blockize = 3), но я не могу умножить [3] [2] * [2][3].

Кто-нибудь знает способ сделать это?Это мое ядро:

#include <stdio.h>

#include <limits.h>

#include <stdlib.h>
#define blocksize 3
#define HM (1*blocksize) 
#define WM (2*blocksize) 
#define WN (1*blocksize)
#define HN WM 
#define WP WN   
#define HP HM  
#define PTH WM
#define PTW HM

__global__ void nonsquare(float*M, float*N, float*P, int uWM,int uWN)

{
__shared__ float MS[blocksize][blocksize];
__shared__ float NS[blocksize][blocksize];


int tx=threadIdx.x, ty=threadIdx.y, bx=blockIdx.x, by=blockIdx.y;
int rowM=ty+by*blocksize;
int colN=tx+bx*blocksize;
float Pvalue=0;


for(int m=0; m< uWM/blocksize;++m){
    MS[ty][tx]=M[rowM*uWM+(m*blocksize+tx)];
    NS[ty][tx]=M[colN + uWN*(m*blocksize+ty)];
    __syncthreads();

    for(int k=0;k<blocksize;k++)
        Pvalue+=MS[ty][k]*NS[k][tx];
    __syncthreads();
    P[rowM*WP+colN]=Pvalue;
     }
    }

Заранее спасибо!

Answer 1

Я думаю, что проще всего было бы просто заполнить блоки на конце нулями:

for(int m=0; m< uWM/blocksize;++m){
    colM = m*blocksize+tx;
    rowN = m*blocksize+ty;
    if (rowM > uWN || rowN > uWM || colM > uWM || colN > uWN) {
        MS[ty][tx]=0.;
        NS[ty][tx]=0.;
    } else {
        MS[ty][tx]=M[rowM*uWM+colM];
        NS[ty][tx]=N[colN + uWN*rowN];
    }

плюс или минус.(Эта строка NS должна ссылаться на N, а не на M, верно?)

Но, поскольку я, похоже, здесь единственный, кто выступает за использование существующих настроенных библиотек, когда это возможно - почему бы не использовать CUBLAS или МАГМА вместо того, чтобы кататься самостоятельно?Они быстрые и проверены сотнями пользователей.

Answer 2

Основное требование к производительности здесь заключается в том, чтобы первое или второе измерение «плитки» совместно используемой памяти было круглым кратным 16 - исторически это то, что необходимо для достижения оптимальной глобальной пропускной способности памяти (т. Е. Слитых транзакций с половиной деформации),Должно ли это быть первое или второе измерение мозаичного элемента, определяется тем, хранятся ли матрицы в главном порядке столбцов или строк.Ничто не говорит о том, что мозаика совместно используемой памяти должна быть квадратной, только то, что начальный размер хранилища (LDA в обозначении BLAS) должен быть кратен 16.

Вы можете легко шаблонизировать ядро ​​с размерами мозаики.в качестве аргументов шаблона и создания нескольких версий в зависимости от размеров матрицы.Для данной архитектуры, вероятно, существует оптимальный размер элемента мозаичного изображения, который уравновешивает занятость и параллелизм на уровне команд.«Умный» способ решить эту проблему, вероятно, состоит в том, чтобы разложить умножение матриц на две операции: первая выполняет большую часть работы при оптимальном размере тайла, а вторая - при другом размере остальных столбцов.Если результат сразу возвращается к памяти хоста после завершения продукта, лучше всего выполнить вторую операцию на хосте, используя оптимизированный BLAS, перекрывающийся с ядром графического процессора.Именно такой подход используют многие подпрограммы в библиотеке UTK Magma.