上星期課程的作業是  Performance of sqrt in CUDA,老師希望我們可以實作完成並做出類似的比較圖,我用了CUDA完成後也做了一份在Xeon Phi上面運作,最後附上Code與比較圖做點筆記,這個程式碼主要實作:有N個點,計算自己與所有點的距離之和,也就是說需要透過兩點平方和開根號求距離的公式,在把所有答案加在一起,就完成了。

這個工作其實不難,在CPU上面應該這樣運行:

void getDistCPU(float *h_result, float *h_X, float *h_Y, int nElems)
{	
	int i,j;
	for ( i=0 ; i < nElems; i++)
	{
		
		
		float localX = h_X[i];
		float localY = h_Y[i];

		float accumDist = 0;
		for ( j=0; j < nElems; j++)
		{

			float curX = h_X[j]-localX;
			float curY = h_Y[j]-localY;

			accumDist += sqrt((curX * curX) + (curY*curY));
		}
		h_result[i] = accumDist;
	}
}

程式碼與參考網址一模一樣,做的事情也如果所說的XD,但若要移植到CUDA上面呢?

這個Function擁有兩層Loop,移植到CUDA上面最簡單的方式就是讓外層Loop交給CUDA Threads去做,也就是外層Loop要做幾次,就要幾個Threads去做,Kernel function如下:

__global__ void CalculateDist(float *pX, float *pY, float *pResult, int totalElements)
{
	int index = blockIdx.x * blockDim.x + threadIdx.x;

	
	float localX = pX[index];
	float localY = pY[index];

	float accumulatedDist = 0;
	
	for (int i=0; i < totalElements; i++)
	{
		float curX = pX[i];
		float curY = pY[i];

		curX = curX - localX;	//curX is now distance
		curY = curY - localY;

		accumulatedDist += sqrt((curX * curX) + (curY*curY));
	}
	pResult[index] = accumul

我對CUDA其實不太熟悉(練習太少有點慚愧Orz),這個簡單的範例還算讓我更加熟悉了一點CUDA的平行方式與用法,Kernal function寫完之後,重點在於要去呼叫它時的Thread數量與Block數量:

        cudaMalloc((void **)&d_result, sizeof(float) * nElems );
	cudaMalloc((void **)&d_x, sizeof(float) * nElems );
	cudaMalloc((void **)&d_y, sizeof(float) * nElems );

	
        int blockNum = nElems / 1024;
    
		
	cudaEvent_t start,stop;
	cudaEventCreate(&start);
	cudaEventCreate(&stop);	
	
	cudaEventRecord(start, 0);
	
	cudaMemcpy(d_x, x, sizeof(float) * nElems , cudaMemcpyHostToDevice);
	cudaMemcpy(d_y, y, sizeof(float) * nElems , cudaMemcpyHostToDevice);
        CalculateDist<<<blockNum,1024>>>(d_x,d_y,d_result,nElems );
	cudaMemcpy(result, d_result, sizeof(float) * nElems , cudaMemcpyDeviceToHost);
	
	cudaEventRecord(stop, 0);
        cudaEventSynchronize(stop);
	
	
	float elapsedTime = 0;
	cudaEventElapsedTime(&elapsedTime, start, stop);

上面的Code包含了計時以及Memory allocate、copy in、copy out,較重要的Call kernel function的部份,Block數量我是用總Element數量去除以最大Threads數1024(不過先決條件就是要整除拉!不然就會出錯的),由於N的大小有我決定所以測試過程中是不會出錯的。
無標題

這是秒數比較圖,單位其實不是很清楚XD,不過可以看到的就是,CPU的運算時間是成對數成長的,因為O(N*N),所以運行時間會因為N的增加而影響更大,CUDA則否,在超過最大的Threads為止(問了學長差不多是65535*1024這麼多),都只要進行O(N),時間當然就差很多拉!

最後補上運行在Phi的程式碼,實作方式為Cilk offload:

//首先在全域變數宣告儲存x,y與result的空間

_Cilk_shared float x[N];
_Cilk_shared float y[N];
_Cilk_shared float result[N];

接下來把CPU的Function改成:

_Cilk_shared void getDistPhi() 
{
	
#ifdef __MIC__
	printf("In MICn");
	int i,j;
	
#pragma omp parallel for
	for ( i=0 ; i < N; i++)
	{
		
		
		float localX = x[i];
		float localY = y[i];

		float accumDist = 0;
#pragma omp simd
		for ( j=0; j < N; j++)
		{

			float curX = x[j]-localX;
			float curY = y[j]-localY;

			accumDist += sqrt((curX * curX) + (curY*curY));
		}
		result[i] = accumDist;
	}
#else
printf("Offload to coprocessor failed!n"); 
#endif
}

如此一來可以避免無法Offload的機器運行這個程式碼而出錯,重點差別在於在外層Loop加上了 #pragma omp parallel for,以及在內層迴圈加上了#pragma omp simd指令,並且透過icc編譯器來編譯則可透過openmp以及cilk offload的方式來進行Phi的平行化運算。

clock_t t1, t2;
printf("Start offload. n");		
double startTime = dsecnd();
_Cilk_offload  getDistPhi();
printf("Offload done. n");
double endTime = dsecnd();
runtime = endTime-startTime;

最後透過這樣的程式碼來呼叫Offload function,這樣的方式又稱為MYO – (Virtual-Shared) Memory Model ,這個計時的方法是mkl的,所以compile的時候記得加上-mkl的選項。

 

就這樣!收工!

 

Reference:

Performance of sqrt in CUDA

Parallel Programming and Optimization with Intel Xeon Phi Coprocessors, Colfax 2013 http://www.colfax-intl.com/nd/xeonphi/book.aspx .