問題描述：求這樣一個step*16矩陣中每列的最小值,並存到第一行返回輸出。

具體情況：初始想法希望通過歸約演算法實現，但是在cuda的核函式中出現瞭如下問題。

錯誤程式碼：

呼叫該核函式，在核函式內自己迴圈4次，則計算結果出錯（i的步長是2，限制條件是<9)。

__global__ void qumin_error (float *idata_dev)
{
	int blockId = blockIdx.x + blockIdx.y * gridDim.x;
	int threadId = blockId * (blockDim.x * blockDim.y)+ (threadIdx.y * blockDim.x) + threadIdx.x;
	unsigned int index = threadId;
    int step=3360*64;
	for (int i = 1; i <9; i =i* 2)//16組資料，進行4次迴圈歸約
	{
		if ((index % (2 * i * step))<(step* i))
		{
			if (idata_dev[index] > idata_dev[index + (i)* step])
			{
				idata_dev[index] = idata_dev[index + (i)* step];
 			}
		}
	}
}
int main(void)
{
		qumin_error << <dimgrid1, dimblock >> >(idata_dev_min);//錯誤呼叫方法
}

為了解決這個問題，本文給出了三種優化方法得到正確答案，三種方法的實現逐步加深，層層遞進，可依次檢視。

優化1：迴圈多次呼叫核函式

核函式內只做一次歸約運算，主函式多次單獨呼叫下面的函式qumin_simple函式，也就是把原來核函式的迴圈拿了出來，結果正確。

__global__ void qumin_proj1(float *idata_dev,int num)
{
	int blockId = blockIdx.x + blockIdx.y * gridDim.x;
	int threadId = blockId * (blockDim.x * blockDim.y)+ (threadIdx.y * blockDim.x) + threadIdx.x;
	unsigned int index = threadId;
	for (int i = num; i <num +1; i = i * 2)
	{	
		if ((index % (2 * i * Step_Sum))<(Step_Sum * i))
		{
			if (idata_dev[index] > idata_dev[index + (i)* Step_Sum])
			{
				idata_dev[index] = idata_dev[index + (i)* Step_Sum];
 			}
		}
	}
}
int main(void)
{
	for (int i = 1; i < 9; i=i*2)
	{
		qumin_proj1<< <dimgrid1, dimblock >> >(idata_dev_min,i);//正確呼叫方法
	}
}

下面記錄一下正確結果的前20個值：

​

錯誤結果的前20個值：

​

其中：呼叫核函式語句沒有任何改變，沒有具體考慮執行緒數與資料量的對應關係。

dimgrid1=(2,16*105) ; dimblock=(32,32)

由於核函式的功能是歸約求最小值，我們可以看出正確結果的前20個值比錯誤結果的更小，可以猜想錯誤函式可能並沒有遍歷到所有待比較資料。

優化2：改變呼叫架構

為了處理此問題，我們首先選擇了一種優化演算法，它放棄了歸約求和的邏輯架構。

這是一個thread執行緒求16個數的最小值的運算。由於Step_Sum=64*3360,一個block無法處理這麼多資料，我們設定一個block中有1024個thread，一個grid中，有210個block(1024*210=step),block0可以計算出前1024列資料的最小值，以此類推。

並且我們用tmp作為中間變數去儲存計算過程中的最小值，最後再將最值資料賦給idata_dev，這樣得到的訪存效率是最高的。這種計算架構雖然迴圈次數增加了，但是對硬體的計算效能發揮到了最大，比原始歸約演算法更高效。

__global__ void qumin_proj2(float *idata_dev)//最優解法
{
	int blockId = blockIdx.x;//210
	int threadId = blockId * (blockDim.x * blockDim.y) + threadIdx.x;// (0~240)*1024+1024	
	unsigned int index = threadId;
	float  tmp;
	tmp = idata_dev[index];
	for (int i = 1; i <16; i++ )
	{
	       if (tmp > idata_dev[index + (i)* Step_Sum])
		    { 
		     tmp = idata_dev[index + (i)* Step_Sum];				
		    }			
	}
	idata_dev[index] = tmp;
	return;
}
int main(void)
{
qumin_proj2<<<210,1024>>>(idata_dev_min);
}

反思：如果一定要用歸約，如何實現？

應該是一個數組中資料量太大，在裝置端上各執行緒的計算速度不同，在執行中，速度快的執行緒可能會讀取速度慢的執行緒還未寫好的那塊記憶體空間。這是很危險的操作。此時，我們可以選擇用共享記憶體，利用同一執行緒塊的執行緒的同步來抑制這種情況。

優化3： 正確的歸約運算實現

正確實現歸約演算法的核心思想是：再一次迴圈內，要讓一個block中的執行緒去做這一次資料歸約，不同block中的thread絕對不能讀寫重複地址的資料，從而避免髒讀髒寫。

共享記憶體的適用情況：對於本機，sharedmemoryperblock=48k，而每個float為32bit=4byte，所以我們可以初始化的shared陣列的最大為48*1024/4=12288.。只要不超出這個數就都可以，在這裡為了方便計算我選擇sharedmemory大小為1024。

目前有兩種設定block數和thread數的方案：

1.blockDim.x=1;threadDim.x=1024

2.blockDim.x=3360;threadDim.x=1024

3.blockDim.x=3360*16；threadDim.x=64

1.方案1：每一個block中的資料個數剛好為1024個，他每次只處理16組數中每64個數據的最小值。迴圈3360次，實現功能。

如下圖所示：每一列表示一組資料，每一次迴圈都是將16組資料的最小值放到黃色的第一列中。

2.方案2是方案1的優化，將方案1的3360次迴圈改成了呼叫3360個block，每3360/16=210個block負責一組資料，每一個block中的資料個數剛好為1024個，所以每一個thread負責一個數據，它會將這個資料和後面15組中同位置的資料進行歸約比較，將最小值放到第一個位置處。

經過試驗：方案1、2成功。

程式碼：

方案1：

__global__ void qumin_project3_1(float *idata_dev)
{
__shared__ float cache[1024];
int blockId = blockIdx.x;
int threadId = blockId * blockDim.x + threadIdx.x;
int tid = threadIdx.x; //(0,63)
int shared_step = 64;
int shang = tid / shared_step;//一共1024個執行緒，對64取商相當於計算該執行緒應該負責第幾組資料
int yushu = tid % shared_step;//對64取餘相當於計算該執行緒應該負責這組資料中64個數中的哪一個
int cache_idx = tid + shang * 64;
for (int j = 0; j<3360; j++)
{
	//if (64 * j < yushu<64 * (j + 1))
		cache[tid] = idata_dev[shang * Step_Sum + yushu+j*64];//store the data into cache [64*16]
	for (int i = 1; i <9; i*=2)
	{
		if (tid % (2 * i*shared_step) < (shared_step*i))
		{
			if (cache[tid] > cache[tid + (i)* shared_step])
			{
				cache[tid] = cache[tid + (i)* shared_step];//compare and change the num
			}
		}
			__syncthreads();
			if (tid<64)
				idata_dev[tid + j*64] = cache[tid];//put the min data back into idata_dev
	}
}
}
void main()
{
qumin_type2 << <1, 1024 >> >(idata_dev_min);
}

​

測試成功後，我嘗試了一下去掉__syncthreads()這條執行緒同步命令，輸出結果就和先前的錯誤結果相同，所以我們可以斷定：前面歸約演算法的問題出在執行緒同步上。

方案2：優化實現很簡單，只是將方案1中的外層迴圈j，換成blockIdx.x就可以了

__global__ void qumin_project3_2(float *idata_dev)
{
	__shared__ float cache[1024];
	int blockId = blockIdx.x;
	int threadId = blockId * blockDim.x + threadIdx.x;
	int tid = threadIdx.x; //(0,63)
	int index = threadId;//(0,63)(64,127)...(3359*64,3360*64-1)
	int shared_step = 64;
	int yushu = tid % shared_step;
	int shang = tid / shared_step;
	int cache_idx = tid + shang * 64;                            
		
    cache[tid] = idata_dev[shang * Step_Sum + yushu + blockIdx.x*shared_step];//store the data into cache [64*16]
		for (int i = 1; i <9; i *= 2)
		{
			if (tid % (2 * i*shared_step) < (shared_step*i))
			{
				if (cache[tid] > cache[tid + (i)* shared_step])
				{
					cache[tid] = cache[tid + (i)* shared_step];//compare and change the num
				}
			}
			__syncthreads();
			if (tid<64)
				idata_dev[tid + blockIdx.x * 64] = cache[tid];//put the min data back into idata_dev
		}
}
void main()
{
qumin_type2 << <3360, 1024 >> >(idata_dev_min);
}

初始bug出現原因思考：

具體衝突表現描述為：折半歸約的過程中，黃色區域與紅色區域比較大小，將小值寫入黃色區域。在歸約的下一次迴圈中，黃色區域要去和綠色區域比較，但是由於硬體限制，綠色區域的某些執行緒仍然沒有完成上一輪次迴圈的資料寫入，就被讀取走了資料，導致發生錯誤。這樣得到的最小值是有問題的。

在原始程式中，我們雖然加上了__syncthreads()命令，但是該命令同步的是同一個執行緒塊內的資料，原情況中是由不同的執行緒塊之間計算最值，所以__syncthreads()無法同步這些計算執行緒。這也解釋了為什麼優化1中，只是確保了上一次寫完了在進行下一次讀取和計算就能得到正確結果。