很簡單,別人不會的你會。
比如,Bitonic Sort(雙調排序,時間復雜度 O(n(logn)^2) )是一種並列排序演算法,特別適合在GPU上實作。我自己用OpenCL實作了一個並進行了調優,演算法是我用紙筆推導的,分為5個kernel。在GPU上編程不同於CPU,沒有偵錯程式可用,也沒有printf,陣列越界電腦直接黑屏。
我就直接用C++實作了一個在CPU上執行的OpenCL kernel,實作了__local限定符 get_global_id get_local_id barrier等函數,用它實作的我的排序演算法。這樣就可以很方便的進行偵錯了,驗證完成後移植在GPU上。
本人看過多個Bitonic Sort的實作, 在github搜尋的基本上就是一堆玩具,沒有我實作的好,Microsoft DirectX SDK中的也是玩具。DirectX-Graphics-Samples的 使用了Indirect模式,但是其numthreads是寫死了的。 CUDA SDK Samples的 4個Kernel ,還算可以。
相比於其他的實作,我的這個使用#define設定不同的local size,在不同的硬件上實作更好的效能。
這就是Bitonic Sort的kernel程式碼(完全是自己實作的):
#if VEC_STEP == 4
#define VEC_SEQ (uint4)(0,1,2,3)
#define SCALAR_TYPE int
#define VEC_TYPE int4
#define VEC_TYPE_AS as_int4
#define LOGICAL_TYPE uint4
#define LOGICAL_TYPE_AS as_uint4
#define INDEX_TYPE uint4
#define INDEX_TYPE_AS as_uint4
#define VLOADN vload4
#define VSTOREN vstore4
#define INPUT_TO_PRIVATE(A,B,IA,IB,INPUT) \
A.s0=INPUT[IA.s0];\
A.s1=INPUT[IA.s1];\
A.s2=INPUT[IA.s2];\
A.s3=INPUT[IA.s3];\
B.s0=INPUT[IB.s0];\
B.s1=INPUT[IB.s1];\
B.s2=INPUT[IB.s2];\
B.s3=INPUT[IB.s3];
#define PRIVATE_TO_INPUT(A,B,IA,IB,INPUT) \
INPUT[IA.s0] = A.s0;\
INPUT[IA.s1] = A.s1;\
INPUT[IA.s2] = A.s2;\
INPUT[IA.s3] = A.s3;\
INPUT[IB.s0] = B.s0;\
INPUT[IB.s1] = B.s1;\
INPUT[IB.s2] = B.s2;\
INPUT[IB.s3] = B.s3;
#else
#error Need VEC_SEQ
#endif
#ifndef LOCAL_VEC_SIZE
#error Need LOCAL_VEC_SIZE
#endif
#ifndef WG_SIZE
#error Need WG_SIZE
#endif
//#define LOCAL_VEC_SIZE 4096
//設 global_size 為 2^m local_size為 2^n 那麽
void
BitnoicSortLocalStage1_1
(
__local
SCALAR_TYPE
*
input
,
const
unsigned
int
num_local
)
{
size_t
gid
=
get_group_id
(
0
);
size_t
num_groups
=
get_num_groups
(
0
);
size_t
lid
=
get_local_id
(
0
);
size_t
lsize
=
get_local_size
(
0
);
#pragma unroll
for
(
unsigned
int
step1
=
1
;
step1
<
LOCAL_VEC_SIZE
;
step1
<<=
1
)
{
unsigned
int
sh
=
0
;
for
(
unsigned
int
step2
=
step1
;
step2
>=
1
;
step2
>>=
1
)
{
for
(
unsigned
int
loop
=
lid
;
loop
<
LOCAL_VEC_SIZE
/
2
;
loop
+=
lsize
*
VEC_STEP
)
{
INDEX_TYPE
iloop
=
VEC_SEQ
*
(
uint
)
lsize
+
loop
;
INDEX_TYPE
gloop
=
iloop
+
(
LOCAL_VEC_SIZE
/
2
)
*
num_local
;
INDEX_TYPE
t1
=
gloop
&
step1
;
INDEX_TYPE
t2
=
t1
^
step1
;
INDEX_TYPE
u1
=
(
iloop
&
(
~
(
step2
-
1
)))
<<
1
;
INDEX_TYPE
u2
=
(
step2
-
1
)
&
iloop
;
INDEX_TYPE
a
=
u1
+
u2
+
(
t1
>>
sh
);
INDEX_TYPE
b
=
u1
+
u2
+
(
t2
>>
sh
);
VEC_TYPE
ia
;
// = input[a];
VEC_TYPE
ib
;
// = input[b];
INPUT_TO_PRIVATE
(
ia
,
ib
,
a
,
b
,
input
);
LOGICAL_TYPE
icomp
=
LOGICAL_TYPE_AS
(
ia
<
ib
);
LOGICAL_TYPE
ncomp
=
LOGICAL_TYPE_AS
(
!
icomp
);
VEC_TYPE
oa
=
select
(
ia
,
ib
,
icomp
);
VEC_TYPE
ob
=
select
(
ia
,
ib
,
ncomp
);
PRIVATE_TO_INPUT
(
oa
,
ob
,
a
,
b
,
input
);
}
barrier
(
CLK_LOCAL_MEM_FENCE
);
sh
++
;
}
}
}
void
BitnoicSortLocalStage1_2
(
__local
SCALAR_TYPE
*
input
,
const
unsigned
int
step1_begin
,
const
unsigned
int
shift
,
const
unsigned
int
num_local
)
{
size_t
gid
=
get_group_id
(
0
);
size_t
num_groups
=
get_num_groups
(
0
);
size_t
lid
=
get_local_id
(
0
);
size_t
lsize
=
get_local_size
(
0
);
unsigned
int
sh
=
shift
;
unsigned
int
step1
=
step1_begin
;
#pragma unroll
for
(
unsigned
int
step2
=
(
LOCAL_VEC_SIZE
/
2
);
step2
>=
1
;
step2
>>=
1
)
{
for
(
unsigned
int
loop
=
lid
;
loop
<
LOCAL_VEC_SIZE
/
2
;
loop
+=
lsize
*
VEC_STEP
)
{
INDEX_TYPE
iloop
=
VEC_SEQ
*
(
uint
)
lsize
+
loop
;
INDEX_TYPE
gloop
=
iloop
+
(
LOCAL_VEC_SIZE
/
2
)
*
num_local
;
INDEX_TYPE
t1
=
step1
&
gloop
;
INDEX_TYPE
t2
=
t1
^
step1
;
INDEX_TYPE
u1
=
(
iloop
&
(
~
(
step2
-
1
)))
<<
1
;
INDEX_TYPE
u2
=
(
int
)(
step2
-
1
)
&
iloop
;
INDEX_TYPE
a
=
u1
+
u2
+
(
t1
>>
sh
);
INDEX_TYPE
b
=
u1
+
u2
+
(
t2
>>
sh
);
VEC_TYPE
ia
;
// = input[a];
VEC_TYPE
ib
;
// = input[b];
INPUT_TO_PRIVATE
(
ia
,
ib
,
a
,
b
,
input
);
LOGICAL_TYPE
icomp
=
LOGICAL_TYPE_AS
(
ia
<
ib
);
LOGICAL_TYPE
ncomp
=
LOGICAL_TYPE_AS
(
!
icomp
);
VEC_TYPE
oa
=
select
(
ia
,
ib
,
icomp
);
VEC_TYPE
ob
=
select
(
ia
,
ib
,
ncomp
);
PRIVATE_TO_INPUT
(
oa
,
ob
,
a
,
b
,
input
);
}
barrier
(
CLK_LOCAL_MEM_FENCE
);
sh
++
;
}
}
//執行1次
__kernel
__attribute__
((
reqd_work_group_size
(
WG_SIZE
,
1
,
1
)))
void
BitonicSortLocalStage1_1Kernel
(
__global
SCALAR_TYPE
*
input
,
const
unsigned
int
step1
,
const
unsigned
int
step2
,
const
unsigned
int
len
)
{
size_t
gid
=
get_group_id
(
0
);
size_t
num_groups
=
get_num_groups
(
0
);
size_t
lid
=
get_local_id
(
0
);
size_t
lsize
=
get_local_size
(
0
);
unsigned
int
loop_count
=
len
/
LOCAL_VEC_SIZE
;
__local
SCALAR_TYPE
value
[
LOCAL_VEC_SIZE
];
for
(
unsigned
int
ng
=
gid
;
ng
<
loop_count
;
ng
+=
num_groups
)
{
async_work_group_copy
(
value
,
input
+
(
ng
*
LOCAL_VEC_SIZE
)
,
LOCAL_VEC_SIZE
,(
event_t
)
0
);
barrier
(
CLK_LOCAL_MEM_FENCE
);
//wait_group_events(1,&e1);
BitnoicSortLocalStage1_1
(
value
,
ng
);
async_work_group_copy
(
input
+
(
ng
*
LOCAL_VEC_SIZE
),
value
,
LOCAL_VEC_SIZE
,(
event_t
)
0
);
barrier
(
CLK_LOCAL_MEM_FENCE
);
//wait_group_events(1,&e2);
}
}
//執行m-n次
__kernel
__attribute__
((
reqd_work_group_size
(
WG_SIZE
,
1
,
1
)))
void
BitonicSortLocalStage1_2Kernel
(
__global
SCALAR_TYPE
*
input
,
const
unsigned
int
step1
,
const
unsigned
int
step2
,
const
unsigned
int
len
,
const
unsigned
int
sh
)
{
size_t
gid
=
get_group_id
(
0
);
size_t
num_groups
=
get_num_groups
(
0
);
size_t
lid
=
get_local_id
(
0
);
size_t
lsize
=
get_local_size
(
0
);
unsigned
int
loop_count
=
len
/
LOCAL_VEC_SIZE
;
__local
SCALAR_TYPE
value
[
LOCAL_VEC_SIZE
];
for
(
unsigned
int
ng
=
gid
;
ng
<
loop_count
;
ng
+=
num_groups
)
{
async_work_group_copy
(
value
,
input
+
(
ng
*
LOCAL_VEC_SIZE
)
,
LOCAL_VEC_SIZE
,(
event_t
)
0
);
barrier
(
CLK_LOCAL_MEM_FENCE
);
BitnoicSortLocalStage1_2
(
value
,
step1
,
sh
,
ng
);
async_work_group_copy
(
input
+
(
ng
*
LOCAL_VEC_SIZE
),
value
,
LOCAL_VEC_SIZE
,(
event_t
)
0
);
barrier
(
CLK_LOCAL_MEM_FENCE
);
}
}
void
BitonicSortLocalStage2
(
__local
SCALAR_TYPE
*
input
)
{
size_t
gid
=
get_group_id
(
0
);
size_t
num_groups
=
get_num_groups
(
0
);
size_t
lid
=
get_local_id
(
0
);
size_t
lsize
=
get_local_size
(
0
);
size_t
gsize
=
get_global_size
(
0
);
#pragma unroll
for
(
unsigned
int
step1
=
LOCAL_VEC_SIZE
/
2
;
step1
>=
1
;
step1
>>=
1
)
{
for
(
int
loop
=
lid
;
loop
<
LOCAL_VEC_SIZE
/
2
;
loop
+=
lsize
*
VEC_STEP
)
{
INDEX_TYPE
iloop
=
VEC_SEQ
*
(
uint
)
lsize
+
loop
;
INDEX_TYPE
u1
=
(
~
(
step1
-
1
)
&
iloop
)
<<
1
;
INDEX_TYPE
u2
=
((
step1
-
1
)
&
iloop
);
INDEX_TYPE
a
=
u1
+
u2
;
INDEX_TYPE
b
=
a
+
step1
;
VEC_TYPE
ia
;
// = input[a];
VEC_TYPE
ib
;
// = input[b];
INPUT_TO_PRIVATE
(
ia
,
ib
,
a
,
b
,
input
);
LOGICAL_TYPE
icomp
=
LOGICAL_TYPE_AS
(
ia
<
ib
);
LOGICAL_TYPE
ncomp
=
LOGICAL_TYPE_AS
(
!
icomp
);
VEC_TYPE
icomp_
=
VEC_TYPE_AS
(
icomp
&
1
);
VEC_TYPE
ncomp_
=
VEC_TYPE_AS
(
ncomp
&
1
);
VEC_TYPE
oa
=
ia
*
icomp_
+
ib
*
ncomp_
;
VEC_TYPE
ob
=
ia
*
ncomp_
+
ib
*
icomp_
;
PRIVATE_TO_INPUT
(
oa
,
ob
,
a
,
b
,
input
);
}
barrier
(
CLK_LOCAL_MEM_FENCE
);
}
}
//執行1次
__kernel
__attribute__
((
reqd_work_group_size
(
WG_SIZE
,
1
,
1
)))
void
BitonicSortLocalStage2Kernel
(
__global
SCALAR_TYPE
*
input
,
const
unsigned
int
step1
,
const
unsigned
int
step2
,
const
unsigned
int
len
)
{
size_t
gid
=
get_group_id
(
0
);
size_t
num_groups
=
get_num_groups
(
0
);
size_t
lid
=
get_local_id
(
0
);
size_t
lsize
=
get_local_size
(
0
);
unsigned
int
loop_count
=
len
/
LOCAL_VEC_SIZE
;
__local
SCALAR_TYPE
value
[
LOCAL_VEC_SIZE
];
for
(
unsigned
int
ng
=
gid
;
ng
<
loop_count
;
ng
+=
num_groups
)
{
async_work_group_copy
(
value
,
input
+
(
ng
*
LOCAL_VEC_SIZE
)
,
LOCAL_VEC_SIZE
,(
event_t
)
0
);
barrier
(
CLK_LOCAL_MEM_FENCE
);
BitonicSortLocalStage2
(
value
);
async_work_group_copy
(
input
+
(
ng
*
LOCAL_VEC_SIZE
),
value
,
LOCAL_VEC_SIZE
,(
event_t
)
0
);
barrier
(
CLK_LOCAL_MEM_FENCE
);
}
}
//執行(m-n) + (m-n)(m-n-1)/2 次
__kernel
__attribute__
((
reqd_work_group_size
(
WG_SIZE
,
1
,
1
)))
void
BitonicSortGlobalStage1Kernel
(
__global
SCALAR_TYPE
*
input
,
const
unsigned
int
step1
,
const
unsigned
int
step2
,
const
unsigned
int
len
,
const
unsigned
int
sh
)
{
size_t
gid
=
get_global_id
(
0
);
size_t
num_groups
=
get_num_groups
(
0
);
size_t
lid
=
get_local_id
(
0
);
size_t
lsize
=
get_local_size
(
0
);
size_t
gsize
=
get_global_size
(
0
);
for
(
unsigned
int
loop
=
gid
*
VEC_STEP
;
loop
<
len
/
2
;
loop
+=
gsize
*
VEC_STEP
)
{
unsigned
int
t1
=
step1
&
loop
;
unsigned
int
t2
=
t1
^
step1
;
unsigned
int
u1
=
((
~
(
step2
-
1
))
&
loop
)
<<
1
;
unsigned
int
u2
=
(
step2
-
1
)
&
loop
;
unsigned
int
a
=
u1
+
u2
+
(
t1
>>
sh
);
unsigned
int
b
=
u1
+
u2
+
(
t2
>>
sh
);
VEC_TYPE
ia
;
VEC_TYPE
ib
;
ia
=
VLOADN
(
a
/
VEC_STEP
,
input
);
//input[a];
ib
=
VLOADN
(
b
/
VEC_STEP
,
input
);
//input[b];
LOGICAL_TYPE
icomp
=
LOGICAL_TYPE_AS
(
ia
<
ib
);
LOGICAL_TYPE
ncomp
=
LOGICAL_TYPE_AS
(
!
icomp
);
VEC_TYPE
oa
=
select
(
ia
,
ib
,
icomp
);
VEC_TYPE
ob
=
select
(
ia
,
ib
,
ncomp
);
VSTOREN
(
oa
,
a
/
VEC_STEP
,
input
);
//input[a];
VSTOREN
(
ob
,
b
/
VEC_STEP
,
input
);
//input[b];
}
}
//執行m-n次
__kernel
__attribute__
((
reqd_work_group_size
(
WG_SIZE
,
1
,
1
)))
void
BitonicSortGlobalStage2Kernel
(
__global
SCALAR_TYPE
*
input
,
const
unsigned
int
step1
,
const
unsigned
int
step2
,
const
unsigned
int
len
)
{
size_t
gid
=
get_global_id
(
0
);
size_t
num_groups
=
get_num_groups
(
0
);
size_t
lid
=
get_local_id
(
0
);
size_t
lsize
=
get_local_size
(
0
);
size_t
gsize
=
get_global_size
(
0
);
for
(
unsigned
int
loop
=
gid
*
VEC_STEP
;
loop
<
len
/
2
;
loop
+=
gsize
*
VEC_STEP
)
{
unsigned
int
u1
=
(
~
(
step1
-
1
)
&
loop
)
<<
1
;
unsigned
int
u2
=
((
step1
-
1
)
&
loop
);
unsigned
int
a
=
u1
+
u2
;
unsigned
int
b
=
a
+
step1
;
VEC_TYPE
ia
;
VEC_TYPE
ib
;
ia
=
VLOADN
(
a
/
VEC_STEP
,
input
);
//input[a];
ib
=
VLOADN
(
b
/
VEC_STEP
,
input
);
//input[b];
LOGICAL_TYPE
icomp
=
LOGICAL_TYPE_AS
(
ia
<
ib
);
LOGICAL_TYPE
ncomp
=
LOGICAL_TYPE_AS
(
!
icomp
);
VEC_TYPE
icomp_
=
VEC_TYPE_AS
(
icomp
&
1
);
VEC_TYPE
ncomp_
=
VEC_TYPE_AS
(
ncomp
&
1
);
VEC_TYPE
oa
=
ia
*
icomp_
+
ib
*
ncomp_
;
VEC_TYPE
ob
=
ia
*
ncomp_
+
ib
*
icomp_
;
VSTOREN
(
oa
,
a
/
VEC_STEP
,
input
);
//input[a];
VSTOREN
(
ob
,
b
/
VEC_STEP
,
input
);
//input[b];
}
}
後來找到工作後,將其改成DirectX Compute Shader實作的。完全獨自一人完成,公司除本人外沒有一個人會。這個是為了做粒子效果用的,而粒子效果需要進行排序。
以後打算使用HLSL Shader Model 6 的新指令(vote, ballot, shuffle)進行效能最佳化,效能可以更快。
在GPU上還有一種排序 Radix Sort(基數排序),Bullet物理引擎用的就是Radix Sort。與Bitonic Sort相比 速度會快一些,但是會占額外的記憶體空間(Bitonic Sort不需要,可以就地進行排序,Radix Sort需要2個Buffer,一個讀一個寫),浮點數比較需要處理(Bullet引擎的FlipFloat)。本人也用OpenCL實作了一個並進行了測試,包含7個Kernel。其中有ExclusiveScan,Histogram等平行計算常用演算法。
=============
2019年1月20日更新:
二叉樹,估計大家都知道是什麽。
但在GPU上建造二叉樹各位聽說過嗎?
Bullet物理引擎就用到了,用於射線相交檢測的BVH,需要將AABB包圍盒轉換成Morton Code後排序,然後根據排序後的Morton Code建造二叉樹。
具體原理請看
和論文 Maximizing Parallelism in the Construction of BVHs, Octrees,and k-d Trees (該連結裏就有)
借用論文中的一張圖:
如圖
紅色代表中間節點 綠色代表葉子節點。一共有8個葉子節點,7個中間節點。
水平的箭頭代表二進制公共字首的高度,如3-4之間就為0,0-1之間為3
讓0為最高的節點,開始尋找最高的公共字首高度,找到3-4之間最高,於是0節點指向3,4節點
3節點 2-3公共字首高度為4,3-4公共字首高度為0,2-3比3-4低所以向左尋找。尋找比3-4高度低的最高的節點。
找到1-2之間公共字首高度為2,最高。3節點指向1,2節點。
4節點 3-4公共字首高度為0,4-5公共字首高度為1,3-4比4-5高所以向右尋找。尋找比3-4高度低的最高的節點。
找到4-5之間公共字首高度為1,最高。4節點指向4,5節點。由於4節點指向自身,因此指向葉子節點。
葉子節點間,中間節點間完全無依賴,可以完全並列的建造。
Bullet物理引擎使用的方法與上述有些不同。
Bullet物理引擎的做法是:
每個葉子節點,如果左邊的公共字首高度比右邊大,那麽就指向左邊的中間節點,如果右邊大就指向右邊的。最左的節點 公共字首高度最小,最右的節點公共字首高度最大。
每個中間節點,分別從左右方向尋找離該節點最近的公共字首高度比該節點高的節點,如果找不到, 公共字首高度為最小。
如果找到的左邊的公共字首高度比右邊大,那麽就指向左邊的中間節點,如果右邊大就指向右邊的。
當時為了搞懂演算法原理,我把執行時產生的數據輸出成文本,才找出了其中的規律。
這就是Bullet物理引擎,一個內建的Demo,看看左下角有多少個立方體。
這個Demo中的拾取演算法(射線-三角形相交)就用到了BVH,每一幀的拾取都會建造BVH。碰撞檢測也可以用BVH,但這個Demo用的是SAP。
寫本文時入職才2個月,第一份工作。
可以看看我的另一個回答:
//=====================
2019年1月27日更新:
看到評論中有些人用過CUDA,我就說一說CUDA 和DirectX12 做計算的區別,難度不是一個等級的:
CUDA分配Buffer 直接用 cudaMalloc 。DirectX12用CreateCommittedResource,再根據資源的使用方式建立view。
CUDA上傳用cudaMemcpy,而DirectX12 有2種方法:第一種直接建立一個 upload buffer ,然後memcpy。第二種建立一個 default buffer,再建立一個upload buffer,先memcpy到upload buffer,再根據資源類別 CopyBufferRegion CopyTextureRegion 到default buffer。
Root Signature :shader的資源繫結布局,相當於函數參數的定義,這個CUDA沒有。
Descriptor Heap: CUDA沒有,存放 view的地方。相當於預先定義的一個陣列,存放Buffer 的指標。
Resource Barrier:資源狀態轉換,CUDA沒有。上傳資源時upload buffer狀態轉換為 COPY_SOURCE ,default buffer狀態轉換為COPY_DEST。當參數使用時default buffer狀態轉換為SHADER_RESOURCE。
Command Queue Command List Command Allocator :所有的渲染命令如 Draw Dispatch Copy Clear 都要記錄到 Command List中,而 Command List需要指定Command Allocator 。完成後Close,Command Queue 再ExecuteCommandList執行記錄的命令。
Fence Signal :ExecuteCommandList 是非阻塞的,執行後立即返回。需要等待的話建立一個Fence, 然後Command Queue再Signal建立的Fence,Fence再Wait。