CUDA ile Matris Çarpımı – 1

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 R
C: 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:

A_{i,j} A matrisinin i. satırının j. sütunundaki eleman
C = AB
C_{i,k}=\sum^N_{j=1}A_{i,j} B_{j,k}
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. 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);

}
&#91;/code&#93;</pre>
<h2>Optimize Edilmemiş GPU Kerneli</h2>
<div>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.</div>
<div>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

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

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.
matMul1 ve matMul2 zamanlaması

Proje dosyasını açabilmek için farklı kaydederek uzantısını 7z yapınız.
Kaynakları indir: CUDA matris carpimi 1