CUDA高性能計(jì)算經(jīng)典問題：歸約

本系列為CUDA進(jìn)階，通過具體的經(jīng)典問題，講述高性能編程的一些基本原則以及方法。建議讀者先閱讀 NVIDIA 官方的編程指南完成 CUDA 入門，基礎(chǔ)比較少的同學(xué)也建議閱讀【工程實(shí)踐 | CUDA優(yōu)化之LayerNorm性能優(yōu)化實(shí)踐】，【CUDA高性能計(jì)算經(jīng)典問題：前綴和】，【GPU架構(gòu)介紹】— https://mp.weixin.qq.com/s/Jq9CFbCgi5cxVlE6Xvu9uw。本文如有不對(duì)的地方歡迎指正。

首先我們不嚴(yán)謹(jǐn)?shù)囟x一下Reduction（歸約）， 給N個(gè)數(shù)值，求出其總和/最大值/最小值/均值這一類的操作，稱為Reduction。如果是使用CPU，我們可以很簡單的寫一個(gè)循環(huán)遍歷一遍即可完成。而在GPU上，我們?nèi)绾尾⑿欣脦浊€(gè)線程去做Reduction？

本文選取求總和為例子編寫代碼，相比內(nèi)存訪問，由于數(shù)值加法并不是很重的計(jì)算，所以這個(gè)問題中，主要注意的是如何利用好各級(jí)Memory的帶寬。

Serial

我們先做一個(gè)Baseline，使用GPU上一個(gè)線程去遍歷得到結(jié)果，如下：

__global__?void?SerialKernel(const?float*?input,?float*?output,?size_t?n)?{??
??float?sum?=?0.0f;??
??for?(size_t?i?=?0;?i?????sum?+=?input[i];??
??}??
??*output?=?sum;??
}??
??
void?ReduceBySerial(const?float*?input,?float*?output,?size_t?n)?{??
??SerialKernel<<<1,?1>>>(input,?output,?n);??
}??

其中n的值為4 * 1024 * 1024，也即輸入的物理大小為4MByte。在作者的環(huán)境里，這段代碼耗時(shí)為 100307us，我們后續(xù)的算法可以與這個(gè)作為對(duì)比。

TwoPass

對(duì)于并發(fā)reduce這個(gè)問題，可以很容易想到一個(gè)樸素解法，以n=8為例，如下

我們把n均分為m個(gè)part，第一步啟動(dòng)m個(gè)block計(jì)算每個(gè)part的reduce結(jié)果，第二步啟動(dòng)一個(gè)單獨(dú)的block匯總每個(gè)part的結(jié)果得到最終結(jié)果。其中每個(gè)block內(nèi)部再把其負(fù)責(zé)的部分均分到每個(gè)線程，這樣就可以得到一個(gè)樸素的代碼如下

__global__?void?TwoPassSimpleKernel(const?float*?input,?float*?part_sum,??
????????????????????????????????????size_t?n)?{??
??//?n?is?divided?to?gridDim.x?part??
??//?this?block?process?input[blk_begin:blk_end]??
??//?store?result?to?part_sum[blockIdx.x]??
??size_t?blk_begin?=?n?/?gridDim.x?*?blockIdx.x;??
??size_t?blk_end?=?n?/?gridDim.x?*?(blockIdx.x?+?1);??
??//?after?follow?step,?this?block?process?input[0:n],?store?result?to?part_sum??
??n?=?blk_end?-?blk_begin;??
??input?+=?blk_begin;??
??part_sum?+=?blockIdx.x;??
??//?n?is?divided?to?blockDim.x?part??
??//?this?thread?process?input[thr_begin:thr_end]??
??size_t?thr_begin?=?n?/?blockDim.x?*?threadIdx.x;??
??size_t?thr_end?=?n?/?blockDim.x?*?(threadIdx.x?+?1);??
??float?thr_sum?=?0.0f;??
??for?(size_t?i?=?thr_begin;?i?????thr_sum?+=?input[i];??
??}??
??//?store?thr_sum?to?shared?memory??
??extern?__shared__?float?shm[];??
??shm[threadIdx.x]?=?thr_sum;??
??__syncthreads();??
??//?reduce?shm?to?part_sum??
??if?(threadIdx.x?==?0)?{??
????float?sum?=?0.0f;??
????for?(size_t?i?=?0;?i???????sum?+=?shm[i];??
????}??
????*part_sum?=?sum;??
??}??
}??
void?ReduceByTwoPass(const?float*?input,?float*?part_sum,?float*?sum,??
?????????????????????size_t?n)?{??
??const?int32_t?thread_num_per_block?=?1024;??//?tuned??
??const?int32_t?block_num?=?1024;?????????????//?tuned??
??//?the?first?pass?reduce?input[0:n]?to?part[0:block_num]??
??//?part_sum[i]?stands?for?the?result?of?i-th?block??
??size_t?shm_size?=?thread_num_per_block?*?sizeof(float);??//?float?per?thread??
??TwoPassSimpleKernel<<>>(input,??
?????????????????????????????????????????????????????????????????????part,?n);??
??//?the?second?pass?reduce?part[0:block_num]?to?output??
??TwoPassSimpleKernel<<<1,?thread_num_per_block,?shm_size>>>(part,?output,??
?????????????????????????????????????????????????????????????block_num);??
}??

這種分為兩步的方法就稱為Two-Pass，這個(gè)方法的時(shí)間為92us，相比之前的100307us確實(shí)快了許多。但是這個(gè)方法仍然比較樸素，沒有利用好GPU特性。

首先讀取Global Memory計(jì)算單線程的結(jié)果時(shí)，由于給單個(gè)線程劃分了一塊連續(xù)地址進(jìn)行局部reduce，導(dǎo)致了同一個(gè)warp內(nèi)的不同線程任意時(shí)刻讀取的地址非連續(xù)。

打個(gè)比方，假設(shè)我們有9個(gè)數(shù)，使用3個(gè)線程去做局部reduce, 0號(hào)線程處理第0,1,2的數(shù)，1號(hào)線程處理第3,4,5的數(shù)，而第三個(gè)線程處理第6,7,8的數(shù)。于是有如下表：

由于從Global Memory到SM的數(shù)據(jù)傳輸也是類似CPU Cache Line的方式，比如一次讀取32*4=128字節(jié)，如果同一Warp內(nèi)所有線程同時(shí)訪問Global Memory連續(xù)的且對(duì)齊的128字節(jié)，那么這次讀取就可以合并為一個(gè)Cache Line的讀取。而在上面的情況下，由于每個(gè)線程讀取的地址都不連續(xù)，意味著每次要觸發(fā)多個(gè)CacheLine的讀取。但是從Global Memory到SM的帶寬并不是無限的，SM內(nèi)以及L2能緩存下的數(shù)據(jù)也不是無限的，這意味著很可能會(huì)導(dǎo)致實(shí)際讀了多次Global Memory，同時(shí)帶寬也被擠壓導(dǎo)致延時(shí)上升。

需要注意，事實(shí)上GPU也支持32/64字節(jié)大小的CacheLine，詳情參考：https://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#device-memory-accesses_

為了應(yīng)對(duì)上文這種現(xiàn)象，我們應(yīng)努力使得Warp內(nèi)不同線程訪問的地址連續(xù)，這會(huì)導(dǎo)致單個(gè)線程處理的地址不連續(xù)，在直覺上與CPU性能優(yōu)化相反。仍然使用9數(shù)3線程的例子，列表如下：

我們可以寫出優(yōu)化后的kernel：

__global__?void?TwoPassInterleavedKernel(const?float*?input,?float*?part_sum,??
?????????????????????????????????????????size_t?n)?{??
??int32_t?gtid?=?blockIdx.x?*?blockDim.x?+?threadIdx.x;??//?global?thread?index??
??int32_t?total_thread_num?=?gridDim.x?*?blockDim.x;??
??//?reduce??
??//???input[gtid?+?total_thread_num?*?0]??
??//???input[gtid?+?total_thread_num?*?1]??
??//???input[gtid?+?total_thread_num?*?2]??
??//???input[gtid?+?total_thread_num?*?...]??
??float?sum?=?0.0f;??
??for?(int32_t?i?=?gtid;?i?????sum?+=?input[i];??
??}??
??//?store?sum?to?shared?memory??
??extern?__shared__?float?shm[];??
??shm[threadIdx.x]?=?sum;??
??__syncthreads();??
??//?reduce?shm?to?part_sum??
??if?(threadIdx.x?==?0)?{??
????float?sum?=?0.0f;??
????for?(size_t?i?=?0;?i???????sum?+=?shm[i];??
????}??
????part_sum[blockIdx.x]?=?sum;??
??}??
}??

這個(gè)優(yōu)化把時(shí)間從92us降低到了78us。

緊接著，更進(jìn)一步，我們可以看到之前的代碼，在把shared memory歸約到最終值時(shí)采取了單線程遍歷的簡單方法，接下來我們優(yōu)化這個(gè)步驟。

這里回顧一下shared memory的特性，其存在于SM上，意味著極快地訪問延時(shí)與帶寬，但其被分成32個(gè)Bank，與Warp的32線程對(duì)應(yīng)。如果一個(gè)Warp內(nèi)的32線程同時(shí)訪問了32個(gè)不同的bank，也即沒有任意兩個(gè)線程訪問同一bank，這時(shí)達(dá)到了最快的訪存速度。否則，如果有兩個(gè)線程同時(shí)訪問了同一個(gè)bank，那么就會(huì)發(fā)生bank conflict，對(duì)這個(gè)bank的訪存無法并發(fā)，形成順序執(zhí)行，也就意味著降低了訪存速度。

如果訪問了互斥的bank，那一定不會(huì)有bank conflict。如果訪問了相同的bank，但是訪問的是bank內(nèi)的同一連續(xù)地址空間，也不會(huì)有bank conflict，這種情況下，如果都為讀操作，則會(huì)廣播給訪問線程們，如果是寫，則只有一個(gè)線程的寫會(huì)成功，具體是哪個(gè)線程是未定義行為。

Shared Memory有4字節(jié)模式和8字節(jié)模式：

• 4字節(jié)模式：其中屬于Bank 0的地址有[0, 4), [128, 132), [256, 260)...，而屬于Bank 1的地址有[4, 8), [132, 136), [260, 264) ...，依次類推每個(gè)bank的地址。

• 8字節(jié)模式：其中屬于Bank 0的地址有[0, 8), [256, 264), [512, 520)...，而屬于Bank 1的地址有[8, 16), [264, 272), [520, 528) ...，依次類推每個(gè)bank的地址。

現(xiàn)在回到我們的問題，將存在于shared memory上的數(shù)據(jù)reduce，同時(shí)盡可能避免bank conflict，一般來說我們假設(shè)單個(gè)block的線程數(shù)是32的倍數(shù)，當(dāng)然我們代碼里實(shí)際也是如此，我們可以每次把數(shù)據(jù)的后半部分加到前半部分上，每次讀寫都沒有bank conflict，代碼如下：

__global__?void?TwoPassSharedOptimizedKernel(const?float*?input,??
?????????????????????????????????????????????float*?part_sum,?size_t?n)?{??
??int32_t?gtid?=?blockIdx.x?*?blockDim.x?+?threadIdx.x;??//?global?thread?index??
??int32_t?total_thread_num?=?gridDim.x?*?blockDim.x;??
??//?reduce??
??//???input[gtid?+?total_thread_num?*?0]??
??//???input[gtid?+?total_thread_num?*?1]??
??//???input[gtid?+?total_thread_num?*?2]??
??//???input[gtid?+?total_thread_num?*?...]??
??float?sum?=?0.0f;??
??for?(int32_t?i?=?gtid;?i?????sum?+=?input[i];??
??}??
??//?store?sum?to?shared?memory??
??extern?__shared__?float?shm[];??
??shm[threadIdx.x]?=?sum;??
??__syncthreads();??
??//?reduce?shm?to?part_sum??
??for?(int32_t?active_thread_num?=?blockDim.x?/?2;?active_thread_num?>=?1;??
???????active_thread_num?/=?2)?{??
????if?(threadIdx.x???????shm[threadIdx.x]?+=?shm[threadIdx.x?+?active_thread_num];??
????}??
????__syncthreads();??
??}??
??if?(threadIdx.x?==?0)?{??
????part_sum[blockIdx.x]?=?shm[0];??
??}??
}??

這一步優(yōu)化從78us降低到了46us. 接下來更進(jìn)一步的，如果活躍線程數(shù)少于等于32，也即只剩最后一個(gè)warp時(shí)，我們是不需要block級(jí)別的同步的，因?yàn)閣arp內(nèi)必然同步（注意這里是無分支情況，不需要syncwarp，感謝評(píng)論區(qū)指正，之前的說法有誤導(dǎo)嫌疑），改寫為：

__global__?void?TwoPassWarpSyncKernel(const?float*?input,?float*?part_sum,??
??????????????????????????????????????size_t?n)?{??
??int32_t?gtid?=?blockIdx.x?*?blockDim.x?+?threadIdx.x;??//?global?thread?index??
??int32_t?total_thread_num?=?gridDim.x?*?blockDim.x;??
??//?reduce??
??//???input[gtid?+?total_thread_num?*?0]??
??//???input[gtid?+?total_thread_num?*?1]??
??//???input[gtid?+?total_thread_num?*?2]??
??//???input[gtid?+?total_thread_num?*?...]??
??float?sum?=?0.0f;??
??for?(int32_t?i?=?gtid;?i?????sum?+=?input[i];??
??}??
??//?store?sum?to?shared?memory??
??extern?__shared__?float?shm[];??
??shm[threadIdx.x]?=?sum;??
??__syncthreads();??
??//?reduce?shm??
??for?(int32_t?active_thread_num?=?blockDim.x?/?2;?active_thread_num?>?32;??
???????active_thread_num?/=?2)?{??
????if?(threadIdx.x???????shm[threadIdx.x]?+=?shm[threadIdx.x?+?active_thread_num];??
????}??
????__syncthreads();??
??}??
??//?the?final?warp??
??if?(threadIdx.x?32)?{??
????volatile?float*?vshm?=?shm;??
????if?(blockDim.x?>=?64)?{??
??????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?32];??
????}??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?16];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?8];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?4];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?2];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?1];??
????if?(threadIdx.x?==?0)?{??
??????part_sum[blockIdx.x]?=?vshm[0];??
????}??
??}??
}??

需要注意，使用warp隱式同步時(shí)使用shared memory需要配合volatile關(guān)鍵字。這個(gè)優(yōu)化從46us降低到了40us。

更進(jìn)一步的，我們可以把第二個(gè)for循環(huán)展開，同時(shí)我們要求在編譯時(shí)就知道blockDim.x，而這可以通過template做到，這里我們?nèi)匀患僭O(shè)blockDim.x是32的倍數(shù)，此外限制其最大是1024，代碼如下：

template???
__global__?void?TwoPassUnrollKernel(const?float*?input,?float*?part_sum,??
????????????????????????????????????size_t?n)?{??
??int32_t?gtid?=?blockIdx.x?*?blockDim.x?+?threadIdx.x;??//?global?thread?index??
??int32_t?total_thread_num?=?gridDim.x?*?blockDim.x;??
??//?reduce??
??//???input[gtid?+?total_thread_num?*?0]??
??//???input[gtid?+?total_thread_num?*?1]??
??//???input[gtid?+?total_thread_num?*?2]??
??//???input[gtid?+?total_thread_num?*?...]??
??float?sum?=?0.0f;??
??for?(int32_t?i?=?gtid;?i?????sum?+=?input[i];??
??}??
??//?store?sum?to?shared?memory??
??extern?__shared__?float?shm[];??
??shm[threadIdx.x]?=?sum;??
??__syncthreads();??
??//?reduce?shm??
??if?(block_thread_num?>=?1024)?{??
????if?(threadIdx.x?512)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?512];??
????}??
????__syncthreads();??
??}??
??if?(block_thread_num?>=?512)?{??
????if?(threadIdx.x?256)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?256];??
????}??
????__syncthreads();??
??}??
??if?(block_thread_num?>=?256)?{??
????if?(threadIdx.x?128)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?128];??
????}??
????__syncthreads();??
??}??
??if?(block_thread_num?>=?128)?{??
????if?(threadIdx.x?64)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?64];??
????}??
????__syncthreads();??
??}??
??//?the?final?warp??
??if?(threadIdx.x?32)?{??
????volatile?float*?vshm?=?shm;??
????if?(blockDim.x?>=?64)?{??
??????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?32];??
????}??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?16];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?8];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?4];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?2];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?1];??
????if?(threadIdx.x?==?0)?{??
??????part_sum[blockIdx.x]?=?vshm[0];??
????}??
??}??
}??

從40us降低到了38us。為了后文方便，我們把這個(gè)對(duì)shared memory reduce的過程提取一個(gè)函數(shù)稱為ReduceSharedMemory，如下：

template???
__device__?void?ReduceSharedMemory(float*?shm,?float*?result)?{??
??if?(block_thread_num?>=?1024)?{??
????if?(threadIdx.x?512)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?512];??
????}??
????__syncthreads();??
??}??
??if?(block_thread_num?>=?512)?{??
????if?(threadIdx.x?256)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?256];??
????}??
????__syncthreads();??
??}??
??if?(block_thread_num?>=?256)?{??
????if?(threadIdx.x?128)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?128];??
????}??
????__syncthreads();??
??}??
??if?(block_thread_num?>=?128)?{??
????if?(threadIdx.x?64)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?64];??
????}??
????__syncthreads();??
??}??
??//?the?final?warp??
??if?(threadIdx.x?32)?{??
????volatile?float*?vshm?=?shm;??
????if?(blockDim.x?>=?64)?{??
??????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?32];??
????}??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?16];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?8];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?4];??
????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?2];??????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?1];??
????if?(threadIdx.x?==?0)?{??
??????*result?=?vshm[0];??
????}??
??}??
}??

Single Pass

上一節(jié)中分為兩次kernel，經(jīng)過多次優(yōu)化，降低到了38us。現(xiàn)在我們想辦法，一次kernel完成我們的目標(biāo)。

首先一個(gè)簡單的辦法，就是把之前的兩步合并到一個(gè)kernel中，為了實(shí)現(xiàn)這個(gè)目標(biāo)，我們需要一個(gè)計(jì)數(shù)器，每個(gè)block完成part sum計(jì)算時(shí)把計(jì)數(shù)器加1，當(dāng)發(fā)現(xiàn)自己加完就是最后一個(gè)block，由這最后一個(gè)block去做part sum的歸約，代碼如下

__device__?int32_t?done_block_count?=?0;

template?
__global__?void?SinglePassMergedKernel(const?float*?input,?float*?part_sum,
???????????????????????????????????????float*?output,?size_t?n)?{
??int32_t?gtid?=?blockIdx.x?*?blockDim.x?+?threadIdx.x;??//?global?thread?index
??int32_t?total_thread_num?=?gridDim.x?*?blockDim.x;
??//?reduce
??//???input[gtid?+?total_thread_num?*?0]
??//???input[gtid?+?total_thread_num?*?1]
??//???input[gtid?+?total_thread_num?*?2]
??//???input[gtid?+?total_thread_num?*?...]
??float?sum?=?0.0f;
??for?(int32_t?i?=?gtid;?i?????sum?+=?input[i];
??}
??//?store?sum?to?shared?memory
??extern?__shared__?float?shm[];
??shm[threadIdx.x]?=?sum;
??__syncthreads();
??//?reduce?shared?memory?to?part_sum
??ReduceSharedMemory(shm,?part_sum?+?blockIdx.x);
??//?make?sure?when?a?block?get?is_last_block?is?true,
??//?all?the?other?part_sums?is?ready
??__threadfence();
??//?check?if?this?block?is?the?last
??__shared__?bool?is_last_block;
??if?(threadIdx.x?==?0)?{
????is_last_block?=?atomicAdd(&done_block_count,?1)?==?gridDim.x?-?1;
??}
??__syncthreads();
??//?reduce?part_sum?to?output
??if?(is_last_block)?{
????sum?=?0.0f;
????for?(int32_t?i?=?threadIdx.x;?i???????sum?+=?part_sum[i];
????}
????shm[threadIdx.x]?=?sum;
????__syncthreads();
????ReduceSharedMemory(shm,?output);
????done_block_count?=?0;
??}
}

void?ReduceBySinglePass(const?float*?input,?float*?part,?float*?output,
????????????????????????size_t?n)?{
??const?int32_t?thread_num_per_block?=?1024;
??const?int32_t?block_num?=?1024;
??size_t?shm_size?=?thread_num_per_block?*?sizeof(float);
??SinglePassMergedKernel
??????<<>>(input,?part,?output,?n);
}

這個(gè)版本的時(shí)間為39us，相比Two-Pass的時(shí)間略有增加。這些時(shí)間的對(duì)比結(jié)論并不是一定的，需要取決于數(shù)據(jù)的規(guī)模，以及gridDim和blockDim的選擇，另外和GPU本身的型號(hào)也有關(guān)系，甚至和CPU的繁忙程度也有關(guān)系，所以數(shù)據(jù)對(duì)比作為參考即可。

另一個(gè)方法就是直接做Atomic：

__global__?void?SinglePassAtomicKernel(const?float*?input,?float*?output,
???????????????????????????????????????size_t?n)?{
??int32_t?gtid?=?blockIdx.x?*?blockDim.x?+?threadIdx.x;??//?global?thread?index
??int32_t?total_thread_num?=?gridDim.x?*?blockDim.x;
??//?reduce
??//???input[gtid?+?total_thread_num?*?0]
??//???input[gtid?+?total_thread_num?*?1]
??//???input[gtid?+?total_thread_num?*?2]
??//???input[gtid?+?total_thread_num?*?...]
??float?sum?=?0.0f;
??for?(int32_t?i?=?gtid;?i?????sum?+=?input[i];
??}
??atomicAdd(output,?sum);
}

這個(gè)方法的耗時(shí)是2553us，作為參考即可。

對(duì)于一個(gè)warp內(nèi)的reduce，我們還可以使用warp級(jí)別的指令直接做這件事情，比如WarpReduce，由于這個(gè)指令不支持float，所以我們退而求其次，使用WarpShuffle也可以做到類似的事情，代碼如下：

template???
__device__?void?ReduceSharedMemoryByShuffle(float*?shm,?float*?result)?{??
??if?(block_thread_num?>=?1024)?{??
????if?(threadIdx.x?512)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?512];??
????}??
????__syncthreads();??
??}??
??if?(block_thread_num?>=?512)?{??
????if?(threadIdx.x?256)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?256];??
????}??
????__syncthreads();??
??}??
??if?(block_thread_num?>=?256)?{??
????if?(threadIdx.x?128)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?128];??
????}??
????__syncthreads();??
??}??
??if?(block_thread_num?>=?128)?{??
????if?(threadIdx.x?64)?{??
??????shm[threadIdx.x]?+=?shm[threadIdx.x?+?64];??
????}??
????__syncthreads();??
??}??
??//?the?final?warp??
??if?(threadIdx.x?32)?{??
????volatile?float*?vshm?=?shm;??
????if?(blockDim.x?>=?64)?{??
??????vshm[threadIdx.x]?+=?vshm[threadIdx.x?+?32];??
????}??
????float?val?=?vshm[threadIdx.x];??
????val?+=?__shfl_xor_sync(0xffffffff,?val,?16);??
????val?+=?__shfl_xor_sync(0xffffffff,?val,?8);??
????val?+=?__shfl_xor_sync(0xffffffff,?val,?4);??
????val?+=?__shfl_xor_sync(0xffffffff,?val,?2);??
????val?+=?__shfl_xor_sync(0xffffffff,?val,?1);??
????if?(threadIdx.x?==?0)?{??
??????*result?=?val;??
????}??
??}??
}??

使用這個(gè)指令的好處就是不需要shared memory，不過我們這個(gè)例子里體現(xiàn)不出優(yōu)勢(shì)來，耗時(shí)沒有什么變化。

最后我們的最優(yōu)耗時(shí)是38us，在本人機(jī)器下測(cè)試了cub的結(jié)果是49us，當(dāng)然這也僅作為一個(gè)參考，具體的影響因素和具體的測(cè)試數(shù)據(jù)規(guī)模也有關(guān)，和GPU型號(hào)等也有關(guān)。不過通過這個(gè)例子，還是把很多知識(shí)點(diǎn)串起來了。

本文如有問題，歡迎指正，歡迎多交流（線上線下均可），共同學(xué)習(xí)共同進(jìn)步，下一篇文章見。

原文鏈接：https://zhuanlan.zhihu.com/p/416959273_

題圖源自Pixabay