Apr 102011
Matris Çarpımı CUDA
Matrislere iki boyutlu vektörler dersek vektörleri tek boyutlu değişkenler gibi görmek, değişkenlere de “boyutsuz” benzetmesini yapmak yerinde olacaktır. Matris nedir sorusunun cevabını aramayacağımız bu yazımızda matris çarpımını C dilinde CPU üzerinde ve GPU üzerinde en basit haliyle nasıl yapabileceğimizi göreceğiz.
Örnek kodlarda kodun basitliği açısından kare matrisler kullanılmıştır.
A: M x N ( M satır, N sütun),B: N x RC: M x R boyutlarında seçilmiştir.Örnek kodlarda matrislerin boyutları #define ile tanımlanacaktır fakat bu yazımızdaki örneklerde
M = N = R olarak seçilmelidir.Tanımlar:
Bu yazımızdaki örnekler:
- C ile yazılmış naif matris çarpımı
- Optimize edilmemiş GPU kerneli
- 1. Optimizasyon
Optimizasyon çok basit olarak yapılmıştır öyleki optimizasyon demeye dilim varmıyor. İleriki yazılarda kernel adım adım optimize edilerek GPU’nun suyu sıkılacaktır.
Naif CPU Matris Çarpımı
C matrisinin her elemanı için A matrisinin N sütunu ve B matrisinin N satırı çarpılarak toplanır.
Bu algoritma, matris çarpımının “text-book” çözümüdür.
M x N x R kere loop döner ( O(n^3), n kare matrisimizin bir boyutu )Bu algoritma, matris çarpımının “text-book” çözümüdür.
inline
float* h_loc(float* h_B, int width, int i, int j)
{
return h_B + i*width + j;
}
void hostMatrixMultiply(
float* h_C,
float* h_A,
float* h_B,
int m,
int n,
int r)
{
clock_t start = clock(), diff;
for(int i=0; i<m; i++)
{
for(int k=0; k<r; k++)
{
float sum = 0;
for(int j=0; j<n; j++)
{
sum += *h_loc(h_A, n, i, j) * *h_loc(h_B, r, j, k);
}
*h_loc(h_C, r, i, k) = sum;
}
}
diff = clock() - start;
int msec = diff * 1000 / CLOCKS_PER_SEC;
printf("hostMatrixMultiply: %i ms\n", msec);
}
Optimize Edilmemiş GPU Kerneli
deviceMultiply1 fonksiyonunun çağırdığı matMul1 kerneli, performans hakkında en ufak tedirginlik hissedilmeden yazılmış bir matris çarpım kernelidir. deviceMultiply1 girizgah fonksiyonunda GPU üzerinde hafıza ayrılmakta, matrisler GPU hafızasına kopyalanmakta ve hesaplanan sonuç CPU hafızasına alındıktan sonra GPU üzerinde ayrılan hafızalar bırakılmaktadır. Sonucun hesaplanması matMul1 kernelinde gerçekleşir.
deviceMultiply1 fonksiyonunda matMul1 kerneli çağrılmadan önce grid ve block boyutları hesaplanır. GPU üzerinde çalışacak thread sayısının C matrisindeki eleman sayısı kadar olmasını sağlayacak değerler hesaplanır. Daha sonra
matMul1<<< dimGrid, dimBlock >>>( d_C, d_A, d_B, m, n, r);
satırıyla kernel çağrılır.
Memory access sayısının düşürülmesi, memory coalescing, tiling gibi hiçbir teknik kullanılmadan yazılmış olan naif CUDA matris çarpım kernelini birlikte inceleyelim.
__global__
void matMul1(
float* d_C,
float* d_A,
float* d_B,
int m,
int n,
int r)
{
int i = blockIdx.x * BLOCK_SIZE + threadIdx.x; // GPU threadinin hesaplayacağı C matrisinin satırı.
int k = blockIdx.y * BLOCK_SIZE + threadIdx.y; // GPU threadinin hesaplayacağı C matrisinin sütunu
int cIdx = i*m + k; // Hesaplanacak elemanın C indisi.
d_C[ cIdx ] = 0;
// C matrisinin her bir hücresi için
for(int j=0; j<n; j++)
{
// her bir eleman için C matrisinin bir elemanı okunur ve yazılır.
d_C[ cIdx ] += *loc( d_A, n, i, j ) * *loc( d_B, r, j, k );
}
}
1. Optimizasyon (Hafıza Erişimi)
matMul1 kernelinde yapılacak çok ufak bir değişiklik ile performans arttırılabilir. Bu değişikliğin getirdiği hız kazancı, GPU’da hafıza erişim işlemlerinin CPU hafızasına erişimden çok daha yavaş olması ile gün ışığına çıkmaktadır.
matMul1 kernelinde d_C, C matrisinin GPU hafızasındaki yerine işaretçidir(pointer). Bu işaretçinin kullanıldığı her yer GPU hafızasına erişmek demektir. İşaretçi başta sıfırlanmıştır ve toplama işlemi işaretçi üzerinde yapılmıştır. += operatörü, d_C[ Idx ] adresinin önce okunması, toplam işleminden sonra da tekrar yazılması demektir. Yani d_A ve d_B’yi saymazsak for döngüsünde toplam 2*n kere hafıza erişimi yapılmaktadır. Halbuki C matrisinin değerini hesaplayıp tek seferde hafızaya yazabiliriz.
matMul2 kernelinde yaratılan değişken ile global GPU RAM erişimi tek thread başına 1 adet yazmaya (d_A ve d_B hariç) düşmüştür.
__global__
void matMul2(
float* d_C,
float* d_A,
float* d_B,
int m,
int n,
int r)
{
int i = blockIdx.x * BLOCK_SIZE + threadIdx.x;
int k = blockIdx.y * BLOCK_SIZE + threadIdx.y;
int cIdx = i*m + k;
float val = 0; // ara toplam değişkeni
// C matrisinin her bir hücresi için
for(int j=0; j<n; j++)
{
val += *loc( d_A, n, i, j ) * *loc( d_B, r, j, k );
}
// C matrisine bir kere yazılır.
d_C[ cIdx ] = val;
}
Zamanlama
Yapılan bu ufak değişiklikle kazanılan 128 x 128 boyutlarında iki matrisin çarpım sürelerini karşılaştırabiliriz.
Proje dosyasını açabilmek için farklı kaydederek uzantısını 7z yapınız.
Kaynakları indir: CUDA matris carpimi 1
Merhaba. Güzel çalışma. Merak ettiğim matMul1<<>>( d_C, d_A, d_B, m, n, r);
satırıyla belirtilen dimgrid ve dimblock değerlerini neye göre belirliyorsunuz. Bunu merak ediyorum. Yardımcı olabilirseniz sevinirim.
int i = blockIdx.x * BLOCK_SIZE + threadIdx.x;
bu satırı kullanmakla tam olarak ne demiş oluyoruz. Burada bu thread indexi i değişkenine atamamızın sebebi nedir?
@ceylan:
dimGrid ve dimBlock değerleri, hedef GPU, problemin yapısal modeli, doluluk oranı(occupancy) gibi parametrelere göre belirlenebilir. Doğrusu, bazı uygun olabilecek değerleri çalışmaya başlamadan test edip dinamik olarak bu değerleri kullanmaktır.
@mustafa:
blocklar BLOCK_SIZE x BLOCK_SIZE boyutlarında matris olsunlar. blockIdx.x, bu matriste yatayda hangi blokta olduğumuzu tanımlar. thread‘ler de blokların içerisindeki lokasyonlar olarak düşünülebilir.
Bu durumda i değeri, yataydaki blokların içindeki thread’e ait olan hücreyi verir.