|
|
楼主 |
发表于 2008-9-13 16:24:00
|
显示全部楼层
Re:请问GPU是如何进行大型矩阵乘法的？
http://www.video.com.cn/club/thread-7636-1-1.html
第二个 CUDA程序
前面介绍的计算平方和的程序，似乎没有什么实用价值。所以我们的第二个 CUDA 程序，要做一个确实有（某些）实用价值的程序，也就是进行矩阵乘法。而且，这次我们会使用浮点数。
, {& B! z+ ^4 o- e4 M6 ?* {' m# w3 q# j+ o/ b' d5 ^
虽然矩阵乘法有点老套，不过因为它相当简单，而且也可以用来介绍一些有关 CUDA 的有趣性质。
4 n5 {: D& c) D) s矩阵乘法
为了单纯起见，我们这里以方形的矩阵为例子。基本上，假设有两个矩阵 A 和 B，则计算 AB = C 的方法如下：4 U5 z& i& ^9 f+ _& d6 I0 v0 O' `
for(i = 0; i < n; i++) {! r* I, b# G* H, j) F) B
for(j = 0; j < n; j++) {, Y# |, j/ l4 q
C[j] = 0;
. Q5 ~- C# e. ~ for(k = 0; k < n; k++) {
& H, I3 _0 G ] C[j] += A[k] * B[k][j];
0 R+ R7 @/ b3 N! A' T9 j$ u }
. q5 w g+ C% Q- Z% I }
+ V0 D. Q I2 ~( w }
/ l$ o/ H) y5 j8 v1 P) z$ s) M% k1 A3 y2 ~5 h5 N+ f [+ C
一开始，我们先准备好产生数据、设定 CUDA 等等的工作：7 E8 Y) X! O. U9 D
int main()
8 C) ~7 M" m# J6 z, h( P/ t& Z# s {" c) b8 n( f- {6 _' Q
float *a, *b, *c, *d;
0 s! t0 v* B8 r g$ a$ p int n = 1000;- \- }0 {- [8 _/ F" n8 f- h
* x- R) F ^) W+ |. p if(!InitCUDA()) return 0;
) p( Q& k! ^8 n0 J T7 m" l; G8 U% V! D
a = (float*) malloc(sizeof(float) * n * n);
) B4 D0 `, @- Y4 Q" _: S b = (float*) malloc(sizeof(float) * n * n);
3 `; y- H; M- v" }1 N6 a- P* B c = (float*) malloc(sizeof(float) * n * n);. H! O1 W+ t0 E% J R7 u' @$ m. M5 U8 S
d = (float*) malloc(sizeof(float) * n * n);* M, h5 E& G- Q# a4 z
, @% F b7 v& [: V/ h srand(0);
" M5 f+ c0 N9 q. r' s# i
$ p9 [4 y% c2 a5 P6 @4 Z& [: q matgen(a, n, n);8 d. o$ x1 O. @) Z; {
matgen(b, n, n);
/ a( ~3 h+ I1 S/ N' m2 v; G& N8 Z/ h9 N6 m
clock_t time = matmultCUDA(a, n, b, n, c, n, n);& e. Y! m: y3 P1 W8 q9 e0 R# K' m
# M$ |; D. L8 m7 S) e6 R
matmult(a, n, b, n, d, n, n);+ r* j: A+ Y; P6 E
compare_mat(c, n, d, n, n);% H3 ]. I: Z# e ~. X1 {
. o9 Q2 ]( V0 C6 i+ Z* i( N+ Y
double sec = (double) time / CLOCKS_PER_SEC;
" h! v0 w; Y/ @ Z- ?+ z! l printf("Time used: %.2f (%.2lf GFLOPS)\n", sec,) i7 z L& Z0 f- u& U
2.0 * n * n * n / (sec * 1E9)); b5 s/ I5 x9 K5 P
: B7 |: M2 P: G- m( e- L2 ~$ G- W5 M return 0;/ F" J& m' R2 o" I0 a8 v9 G
}
2 Z8 M' ]6 v* K" J: GInitCUDA 函式和第一个 CUDA 程序一样，可以直接参考前面的文章。以下是上面用到的一些其它的函式：6 F9 @! X0 |' a& q: H& H
产生矩阵：
/ o! s" g% `5 P$ i. [ void matgen(float* a, int lda, int n)
8 [9 @0 ^* W8 C) w8 D) M {
, N |: l/ B# e int i, j;
) G u. n K. N I) O( x' m1 ^' q5 B$ e
for(i = 0; i < n; i++) {
# Y/ p0 k5 U( y( | for(j = 0; j < n; j++) {* n( [% d$ t( A6 g; o
a[i * lda + j] = (float) rand() / RAND_MAX + & A" t- T v' z! {0 t6 O
(float) rand() / (RAND_MAX * RAND_MAX);8 P% x1 A+ l5 y: o' w+ C/ s& G7 c( E& v
}
$ ^2 v% f {2 I& F3 L( `6 g7 a }/ G3 f+ t; O+ J9 @* Z
}
* E0 R0 Q: z- i2 N: K- A这个函式只是利用随机数生成器把矩阵填满 0 ~ 1 之间的数字。特别注意到因为 C 语言中无法声明变动大小的二维矩阵，所以我们使用 i * lda + j 的方式。$ z \; W5 Z/ r* b8 V" [7 Z
进行矩阵乘法：5 ]6 [2 Q( `' [
void matmult(const float* a, int lda, const float* b, int ldb,+ Q+ v0 ]1 c& b) U' l+ ^: E$ A& g
float* c, int ldc, int n)
8 U' I+ [& D5 [# `" J% L2 Y$ y {* g$ n& C* c5 b+ L" ^+ A
int i, j, k;
7 d F$ h& I% ^, k& d
% l/ O, \- G& H for(i = 0; i < n; i++) {$ ~5 d! B5 I( ]
for(j = 0; j < n; j++) {
& e! A+ y5 w6 U' D# Z double t = 0;, b4 k! S. q; D! ]
for(k = 0; k < n; k++) {
4 M' {" @7 W& _# y5 z3 M t += a[i * lda + k] * b[k * ldb + j];6 l3 i2 x4 D3 c5 |2 f( `0 |
}4 F0 C8 j8 j* d7 K1 J; [( X: t
c[i * ldc + j] = t;
4 O5 Y& j( [5 F+ k. J: H }. k8 W4 n2 c$ I( A1 F/ V- E' g$ e. U
}
7 U1 d4 c8 j* f7 q: L }, o. ~: A; s* T$ X; F6 l; S
这是以 CPU 进行矩阵乘法、用来进行验证答案正确与否的程序。特别注意到它用 double 来储存暂时的计算结果，以提高精确度。
- r) f ?0 l& H8 ?4 O- g/ W0 }2 _( d' A. h0 E `# [3 C
验证结果：
4 p! f v2 ~' X, \& ?4 U void compare_mat(const float* a, int lda,! o% Y* j& u$ r( w; W
const float* b, int ldb, int n)
0 e; ^3 u: l( Z9 A8 v {7 ]! U3 ` @' ]9 a) W8 J/ g. |
float max_err = 0;
, A3 H3 Q. Q7 f) ]. g" } float average_err = 0;& p8 {% w- o, ?4 P$ @+ q
int i, j;. _+ I% Q3 @2 l! x: Q$ A
* `! a" h" | D; |$ s
for(i = 0; i < n; i++) {
' f' u9 j4 ~' I9 J for(j = 0; j < n; j++) {- ~6 {# U! Q3 @/ P1 P+ T3 ]( {
if(b[i * ldb + j] != 0) { l% z) ?" O4 a3 s& Z
float err = fabs((a[i * lda + j] - ( _$ U7 C7 C# O. j
b[i * ldb + j]) / b[i * ldb + j]);
8 k3 e& g8 C5 g if(max_err < err) max_err = err;
, Y7 o0 d C3 x" z! D% k0 l- X: u average_err += err;
9 C9 M8 ^* q8 U$ O4 ] }& O3 g+ E# _ h
}2 \4 |1 Q. m: y8 [, _( W+ U
}+ w" u: J) c- m5 d
& b ~0 X2 |6 O: \* K1 p printf("Max error: %g Average error: %g\n",: g8 Z; I$ ^! T+ c& n4 F8 B
max_err, average_err / (n * n));1 n- x2 U9 x7 e* P
}0 D+ w! A3 _# m! Q
这个函式计算两个矩阵的最大相对误差和平均相对误差，并把结果印出来。" q8 F" T* S! T2 x" h; p8 T
最后是 CUDA 的矩阵乘法的部份：3 c# }; {0 V! f$ M% }4 N; w
#define NUM_THREADS 2565 F' E5 w+ Q& h- M! ~" H0 _% [
( b0 q9 \' l5 h clock_t matmultCUDA(const float* a, int lda,4 k' H) c# I) d# U b' v% x" E, o+ u
const float* b, int ldb, float* c, int ldc, int n)
* {( K- V% W2 c {: @. \1 [" e0 E, }2 \( ?& o7 r
float *ac, *bc, *cc;
, Y" K6 g6 E5 L1 N clock_t start, end;) R) `# g2 C; u. Q3 z/ z2 _2 V/ |
% W f* ]+ D, q5 x U( S+ |$ ^0 D
start = clock();) i* R: j( d; c
cudaMalloc((void**) &ac, sizeof(float) * n * n);
. p! e; a; K- _0 ?- z8 O8 R0 V cudaMalloc((void**) &bc, sizeof(float) * n * n);
) s1 _& G2 L3 b6 y6 d9 L cudaMalloc((void**) &cc, sizeof(float) * n * n);# X( I3 r& Z5 R0 v# r% H
- I: ^$ }9 @; x1 K cudaMemcpy2D(ac, sizeof(float) * n, a, sizeof(float) * lda,
, e7 P7 Z- G& m. ~0 x# P$ K( T- T sizeof(float) * n, n, cudaMemcpyHostToDevice);
! R' Y: J( c6 ]9 N' o5 s- O# _ cudaMemcpy2D(bc, sizeof(float) * n, b, sizeof(float) * ldb,
- s; n" w6 I8 g9 a, d, i" s. y sizeof(float) * n, n, cudaMemcpyHostToDevice);. Q X7 P2 b* s" O; C) ?: v
7 Z" x* Q4 e% N1 w. Y' O7 I int blocks = (n + NUM_THREADS - 1) / NUM_THREADS;
! l0 {* V8 N& q% e. o( M6 y& g matMultCUDA<<<blocks * n, NUM_THREADS>>>8 L @* Z |0 G4 t
(ac, n, bc, n, cc, n, n);: U4 _4 e0 k4 k& r: C: T* e, \/ T
2 |0 [7 \. C( j; C, p
cudaMemcpy2D(c, sizeof(float) * ldc, cc, sizeof(float) * n,+ }6 c @ i$ `' f8 n+ l
sizeof(float) * n, n, cudaMemcpyDeviceToHost); ^6 ?) e3 r6 N/ m3 C& N" f& {9 x4 P
/ h( T# ~$ ~2 t; E! u2 ]9 x/ t
cudaFree(ac);" B1 V+ ^1 p: c# N( A* e2 ?
cudaFree(bc);
9 q7 l0 ^8 K" w) [/ x# R+ h3 s cudaFree(cc);3 K5 b# x0 c" h& u# S- V8 Q
) i* L2 X+ g. a1 X r# s
end = clock();
, @ S0 W) o$ F$ a4 I0 c0 j" C/ `+ p8 \) u O. K/ }
return end - start;7 c3 W3 ~7 o) S( a3 z+ o
}' K1 R$ I& |# d" q% I
这个函式相当单纯，就是在显卡内存中配置存放矩阵的内存，然后把主内存中的矩阵数据复制到显卡内存上。不过，因为我们的矩阵乘法函式可以指定 pitch（即 lda、ldb、和 ldc），所以如果用一般的 cudaMemcpy 函式来复制内存的话，会需要每个 row 都分开复制，那会需要呼叫很多次 cudaMemcpy 函式，会使效率变得很差。因此，在这里我们用了一个新的 cudaMemcpy2D 函式，它是用来复制二维数组，可以指定数组的 pitch。这样就可以透过一次函数调用就可以了。; T+ K0 l8 ~- b. d8 z
进行计算的 kernel 如下：
' W3 t7 C0 I' t6 [( b: l( I __global__ static void matMultCUDA(const float* a, size_t lda,$ b% S1 Z) S; H1 {6 A0 S7 K
const float* b, size_t ldb, float* c, size_t ldc, int n). t9 \" q# o+ L$ x6 r
{
3 k5 h" J1 z7 _+ Q" L1 H const int tid = threadIdx.x;1 ~; |* f( E4 L# O- h8 s% |$ [, e
const int bid = blockIdx.x;4 {, T( F1 |8 s9 s* E
const int idx = bid * blockDim.x + tid;
0 [) e6 d- I; Z( i8 O& B const int row = idx / n;
& }# q _* P0 S( a/ m0 D const int column = idx % n;& x$ @2 u! }, N8 S! Q" j; ^0 q
int i;% D5 |8 y0 C3 V+ t
( \/ g6 [3 B! ^
if(row < n && column < n) {) ^; ?4 o2 R9 j* w
float t = 0;8 G, V2 a% ]0 Y. C. U( c
for(i = 0; i < n; i++) {
" k" T: v5 p+ a: N( H0 Z: z t += a[row * lda + i] * b[i * ldb + column];! g8 i' _% D1 _# h- `
}
4 X) m, N; J; ^0 Z) B, \8 g c[row * ldc + column] = t;
8 ?' u: V _+ P v* ?" O- n }
( A, ?2 I3 j6 V! _0 r }6 t' s4 }9 M$ M& |
这个函式一开始先从 bid 和 tid 计算出这个 thread 应该计算的 row 和 column，在判断 row 和 column 在范围内之后，就直接进行计算，并把结果写到 c 矩阵中，是非常单纯的函式。
' \3 s# h9 R1 Z2 ~6 K V: m! d在 GeForce 8800GT 上实际执行的结果如下：; b9 b1 \5 e4 l0 E7 }
Max error: 2.01484e-006 Average error: 3.36637e-0078 `1 j4 ^% g+ f' C9 R* o
Time used: 1.1560 (1.73 GFLOPS)2 c5 U7 Q5 U/ q5 Q0 T4 \7 Q9 B
可以看到两个问题：
! B7 B9 c6 r$ T) E- f. I) B" a
很明显的，执行效率相当低落。
最大相对误差偏高。理想上应该要低于 1e-6。
计算结果的误差偏高的原因是，在 CPU 上进行计算时，我们使用 double（即 64 bits 浮点数）来累进计算过程，而在 GPU 上则只能用 float（32 bits 浮点数）。在累加大量数字的时候，由于累加结果很快会变大，因此后面的数字很容易被舍去过多的位数。1 d, A5 M5 ?* t* Y
由于 CUDA 的浮点数运算，在进行加、减、乘法时是符合 IEEE 754 规定的精确度的，因此，我们可以利用 Kahan's Summation Formula 来提高精确度。把程序改成：" P% t- x7 c6 M( j' [1 m
if(row < n && column < n) {0 Q7 x0 i" w* F- W% Y8 n
float t = 0;8 a7 G$ q% t0 h( ^
float y = 0;
9 b# C7 Q( U4 {8 X3 q( H( ^' F for(i = 0; i < n; i++) {
6 f# [3 [, }( W9 V8 Z float r;
0 a) y& p' J2 B. L) h ~3 X! l" m y -= a[row * lda + i] * b[i * ldb + column];
% f+ M& `; c( n' H8 o; M r = t - y;2 \9 b ~4 J( |- S
y = (r - t) + y;
! P% z+ e+ K0 g; t' R: K t = r;6 O) K) l0 S/ U7 \$ @; W' {
}
, m: O3 a9 d0 \+ ~ }
2 W& Y' v; w1 d* m修改后的程序的执行结果是：% Z$ f+ W% M& J- w' _0 Q/ X; i
Max error: 1.19209e-007 Average error: 4.22751e-008
3 h% c' G% l2 C1 m+ ?- j- e+ b Time used: 1.1560 (1.73 GFLOPS)8 a! p h. j0 A+ @, Y" E. J/ v
可以看到相对误差有很大的改善，效率则没什么变化。0 p! _" n5 u o9 i
由于 Kahan's Summation Formula 需要的运算量提高，但是效率却没有什么改变，可以看出这个 kernel 主要的瓶颈应该是在内存的存取动作上。这是因为有大量的内存读取是重复的。例如，矩阵 a 的一个 row 在每次进行计算时都被重复读入，但这是相当浪费的。这样的计算方式，总共需要读取 2*n3 次内存。如果让一个 row 只需要读入一次的话，就可以减到为 n3+n2 次。# J# P) f( b* p' ]" c0 F0 o
第一个改良
和我们的第一个 CUDA 程序一样，我们可以利用 shared memory 来储存每个 row 的数据。不过，因为只有同一个 block 的 threads 可以共享 shared memory，因此现在一个 row 只能由同一个 block 的 threads 来进行计算。另外我们也需要能存放一整个 row 的 shared memory。因此，把先把呼叫 kernel 的部份改成：5 d3 M0 R$ @+ t6 {* a$ n
matMultCUDA<<<n, NUM_THREADS, sizeof(float) * n>>>
. g1 i' j+ ?" ^5 o (ac, n, bc, n, cc, n, n);% j5 C& `8 s! Z, R/ a. I5 ^
kernel 的部份则改成：
6 I4 d3 P$ V5 o( I __global__ static void matMultCUDA(const float* a, size_t lda,! T9 s+ |: y2 ^4 J9 R1 Q
const float* b, size_t ldb, float* c, size_t ldc, int n)
7 J- [, s1 t; D: W& m {
* b4 J, c. i7 i: H, u4 u extern __shared__ float data[];# B G0 p- b' `
const int tid = threadIdx.x;
: l: J; q( n. N) T7 i! D: z const int row = blockIdx.x;+ o/ x/ B. b6 ]8 G% _
int i, j;1 b- y* y; D c0 S) z1 j0 q! {9 p
# u- c# \9 V0 X) r7 A
for(i = tid; i < n; i += blockDim.x) {% V) ?8 R8 S4 A \+ E! A0 N
data = a[row * lda + i];
% m9 K7 F4 o1 P2 D' |+ o }8 E6 m( I& H' X, A6 [& f) z ]
6 r' [8 a% v1 G. V# s% M __syncthreads();
! L4 h! y( Q: z, P0 j# ^
! U6 N, J: h/ O8 B4 X7 Y2 l2 K for(j = tid; j < n; j += blockDim.x) {
* D. Q" ?" j/ T6 ` float t = 0;
; s# ?' B/ F- G float y = 0;
" I/ u9 `9 P+ R( I for(i = 0; i < n; i++) {
; y3 L0 Z, h' l. `5 ~: x- f5 ?2 r float r;
& }" s, e B/ A( ^& ] }5 F y -= data * b[i * ldb + j];
( [* B3 G4 f4 V! ?, a7 G4 B r = t - y;
~2 G$ s- D' ]: A$ J3 i y = (r - t) + y;
( S. R( q" X. Z# B" k t = r;
: f, \/ l, O" p# l* W }
{5 [) B/ Y7 G) N5 J c[row * ldc + j] = t;4 p: w: R" }3 d0 n
}0 \: ^2 X3 x1 x+ e% i; Z% d) e2 c
}# ~) _7 z3 C' Y! `. p) t: c* b
@0 Q b, Y7 L) L7 }2 j# b$ V$ o
第一个部份先把整个 row 读到 shared memory 中，而第二个部份则进行计算，并没有太大的变化。主要的差别是现在一个 row 只由一个 block 进行计算。
1 j. s' e. `9 t0 Z- }1 C在 GeForce 8800GT 上，执行的结果是：5 b: t; H# X3 G* o; D$ P
Max error: 1.19209e-007 Average error: 4.22751e-0087 v! X* Y; n- {+ F
Time used: 0.4220 (4.74 GFLOPS)8 E' E2 G# o; C# N1 v8 V. g+ E
很明显的，计算的结果并没有改变，不过速度则提高了超过一倍。虽然如此，但是这样的效率仍不尽理想，因为理论上 GeForce 8800GT 有超过 300GFLOPS 的运算性能。即使是把 Kahan's Summation Formula 所需要的额外运算考虑进去，这样的效率仍然连理论最大值的十分之一都不到。2 [ L2 y" s: e8 A* k6 y8 b
会有这样的结果，原因其实还是同样的：对内存的存取次数太多了。虽然现在 A 矩阵的 row 的数据已经不再需要重复读取，但是 B 矩阵的 column 的数据仍然一直被重复读取。5 F* |% u* o7 i0 \* ~: q/ L
另一个问题比较不是那么明显：对 B 矩阵的读取，虽然看起来不连续，但实际上它是连续的。这是因为不同的 thread 会读取不同的 column，因此同时间每个 thread 读取的各个 column 加起来，就是一个连续的内存区块。那么，为什么效率还是不佳呢？这是因为，GPU 上的内存控制器，从某个固定的倍数地址开始读取，才会有最高的效率（例如 16 bytes 的倍数）。由于矩阵大小并不是 16 的倍数（这里使用的是 1000x1000 的矩阵），所以造成效率不佳的情形。$ O8 y1 o2 H' @( g0 m
要解决这个问题，我们可以在 cudaMalloc 的时候稍微修改一下，让宽度变成 适当的倍数就可以了。但是，适当的倍数是多少呢？幸运的是，我们并不需要知道这些细节。CUDA 提供了一个 cudaMallocPitch 的函式，可以自动以最佳的倍数来配置内存。因此，我们可以把 cudaMalloc 的部份改成：
7 n) D7 F1 d* o+ Q7 D0 s: c size_t pitch_a, pitch_b, pitch_c;; J" y) W1 f( P1 P4 H& cudaMallocPitch((void**) &ac, &pitch_a, sizeof(float) * n, n);
& j7 g7 x m% ^5 q cudaMallocPitch((void**) &bc, &pitch_b, sizeof(float) * n, n);
0 g* U: A; ?+ G y cudaMallocPitch((void**) &cc, &pitch_c, sizeof(float) * n, n);3 w1 _% f. C2 R* S! A4 z+ |1 D
cudaMallocPitch 函式会以适当的倍数配置内存，并把配置的宽度传回。因此，在把矩阵复制到显卡内存上时，要使用它传回的宽度：
) {1 y- q5 g% F' l6 h0 p; b1 U cudaMemcpy2D(ac, pitch_a, a, sizeof(float) * lda,
& z8 y! m4 U( M( P sizeof(float) * n, n, cudaMemcpyHostToDevice);
) L' @7 J( ^5 O \; W5 p cudaMemcpy2D(bc, pitch_b, b, sizeof(float) * ldb,
8 a0 k" u) q2 ~) r$ J sizeof(float) * n, n, cudaMemcpyHostToDevice);* m0 a4 O0 J5 l B5 x6 v
呼叫 kernel 的部份也需要修改：
- y& d# r2 {2 k matMultCUDA<<<n, NUM_THREADS, sizeof(float) * n>>>5 E/ } L& } g6 v9 O
(ac, pitch_a / sizeof(float), bc, pitch_b / sizeof(float),
% t( O( O# b5 `9 p0 H4 F5 Q2 q cc, pitch_c / sizeof(float), n);1 k% k9 T0 n4 f5 [" g Z) L
同样的，把计算结果复制回到主内存时，也要使用传回的宽度值：3 e+ ?. \9 Z' _5 {. O) P
cudaMemcpy2D(c, sizeof(float) * ldc, cc, pitch_c,
9 @8 A) n7 b" N; I( \8 i; ^6 q2 i% ` sizeof(float) * n, n, cudaMemcpyDeviceToHost);% I! A, M F# `$ j, j5 g! X( f) U
这样就完成了。Kernel 部份则不需要修改。
; a! G- p; {$ t1 ], Y4 w H% T这样的修改有多大的效果呢？在 GeForce 8800GT 上执行，结果如下：
5 M3 Y4 f$ I4 w& J" [0 A Max error: 1.19209e-007 Average error: 4.22751e-008
. y( t/ B7 I2 r6 ]6 ]- k0 I Time used: 0.1250 (16.00 GFLOPS): A# m. c% q9 [ e2 C
可以看到，执行速度又再大幅提高了三倍多！而这只是把内存的配置方式稍微修改一下而已。2 \8 o' V1 V' T
虽然执行速度提高了很多，但是，和前面提到的理论值相比，其实还是有相当的差距。这是因为，前面也提到过，这样的做法需要 n3+n2 次的内存读取，和 n2 次的内存写入动作。由于 n = 1000，每个数字的大小是 32 bits，所以总共的内存存取数据量约为 4GB。除以实际执行的时间 0.125 秒，得到的带宽数值是约 32GB/s，这已经接近 GeForce 8800GT 显卡内存的带宽了。由于我们计算时间的时候，把配置内存、以及数据的复制动作也计算进去，因此实际上花费在 kernel 的时间是更短的（约 0.09 秒）。因此，可以很明显的看出，这个程序的效率，是受限于内存带宽的。2 k$ }) p+ A" [. 进一步的改良
上一节的结论显示出，矩阵乘法的程序，效率是受限于内存带宽的。那有没有办法降低内存的存取次数呢？答案当然是有的，不然就不会有这一节了  ; U6 ~1 f3 Q) E# `; N7 Q. p& J
要进一步降低内存带宽的使用，可以注意到，在上一节的方法中，虽然 A 矩阵的存取次数被减至最低，但是 B 矩阵的存取次数并没有减少。这是因为我们只将 A 矩阵的 row 加载到 shared memory 中，但是 B 矩阵的 column 也是有被重复使用的。理想上应该也可以避免重复加载才对。不过，由于 B 矩阵的 column 使用的时机，和 A 矩阵的 row 是不同的，所以并不能直接这样做。! h5 `9 {. |1 P' Q
解决方法是 "blocking"。也就是把整个矩阵乘法的动作，切割成很多小矩阵的乘法。例如，要计算 C 矩阵的 (0, 0) ~ (15, 15) 的值，可以把它想成：
! P7 W Z. @: S) y1 | A(0~15, 0~15) * B(0~15, 0~15) + A(0~15,16~31) * B(16~31, 0~15)* a& I( Z z2 s8 g. _4 H s
+ A(0~15, 32~47) * B(32~47, 0~15) + ...! n1 A( K- W0 l9 M" Q4 \7 这样一来，我们就可以把两个小矩阵加载到 shared memory，则小矩阵本身的乘法就不需要再存取任何外部的内存了！这样一来，假设小矩阵的大小是 k，则实际上需要的内存存取次数就会变成约 2k2(n/k)3 = 2n3/k。 g( i& p+ u v9 E+ n
由于目前 CUDA 每个 block 的 thread 数目最多是 512，因此 k = 16 似乎是一个相当理想的数字（共 256 个 threads）。因此，对于一个 n = 1000 的矩阵来说，我们可以把内存存取的量减少到约 500MB，也就是上一节的存取量的 1/8。理论上，这样应该可以让效率提高八倍（假设没有遇到别的瓶颈）。! V6 ]& V2 l/ l: Q" x
为了方便进行区块的计算，我们让每个 block 有 16x16 个 threads，再建立 (n/16)x(n/16) 个 blocks。把呼叫 kernel 的地方改成：: A5 @* s2 E3 D- k) ?) G5 [1 T
int bx = (n + BLOCK_SIZE - 1) / BLOCK_SIZE;3 i X8 q0 V& g
dim3 blocks(bx, bx);* N) j) n, U+ Z- g& O
dim3 threads(BLOCK_SIZE, BLOCK_SIZE);% [0 O* x+ }; ^; W) x/ h
matMultCUDA<<<blocks, threads>>>(ac, pitch_a / sizeof(float),$ O4 K% c+ D" w
bc, pitch_b / sizeof(float), cc, pitch_c / sizeof(float), n);
$ b" n4 C$ o( g7 D" n3 t$ sBLOCK_SIZE 则是定义成 16。dim3 是 CUDA 的一种数据型态，表示一个 3D 的向量。在这里，我们透过 dim3 来建立 16x16 个 threads 的 block，和 (n/16)x(n/16) 个 blocks。
`1 q5 @ m: q( S# rKernel 程序的部份，则改成：
! q1 Y4 j9 \% ^2 ^' Q, ] N __global__ static void matMultCUDA(const float* a, size_t lda,, k8 R6 G0 D; W4 F) H" H0 t
const float* b, size_t ldb, float* c, size_t ldc, int n)
( E( q: A5 H4 C7 p/ |5 T* g {1 K% }1 R- v# w0 v4 x, s& r+ ]
__shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];
' w0 L9 o% g( b6 l9 J- i+ S __shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];" j% M9 Z% S: Y2 s0 P+ _
const int tidc = threadIdx.x;
+ {0 y+ q6 R: Y9 j I3 n const int tidr = threadIdx.y;# K7 R$ d0 V: W/ x( E. i; Q2 ]
const int bidc = blockIdx.x * BLOCK_SIZE;
+ a- r" S+ G( V6 Q2 ]) K$ \2 u const int bidr = blockIdx.y * BLOCK_SIZE;
4 o( Z5 T$ J, N& q0 L( ^5 d1 y int i, j;6 ? L9 `3 }1 O0 n/ e1 C. Z2 V5 c7 u/ W
4 Y% S; W) P* E$ U' r float results = 0;# g! X4 s$ b3 b# \/ T/ f
float comp = 0;
- V0 ?; l3 }1 K( x
. U: r7 E7 G: x5 F; u for(j = 0; j < n; j += BLOCK_SIZE) {
; |$ z! J$ F$ B3 k if(tidr + bidr < n && tidc + j < n) {$ q2 j* Y" f7 j1 `( L6 G. t+ w
matA[tidr][tidc] = a[(tidr + bidr) * lda + tidc + j];5 h# q6 z, Y1 N# P
}
$ D+ q, M8 B I5 Q else {
c2 ~0 }) M: t) B8 `% ~ matA[tidr][tidc] = 0;
$ E* P7 ?# D- g2 U! d }
% s5 M+ L5 i7 t2 z8 I$ X
/ t6 N9 U& c) D4 G' y2 { if(tidr + j < n && tidc + bidc < n) {- D) Q# w7 [; s2 ?: U
matB[tidr][tidc] = b[(tidr + j) * ldb + tidc + bidc]; O" ?# P7 E8 M0 y9 Y3 s, `
}
" U# [/ O2 W& A" r2 s+ o5 C5 v else {8 V4 ^' D( B+ {3 h
matB[tidr][tidc] = 0;4 I% q8 Y8 A! k
}
& p( r3 m' u4 ~$ n6 M: x
) e) _5 v! i$ V6 \% X/ e __syncthreads();3 O9 R+ ^/ f" G p! J. {7 [
; T' v% J3 V, @: M1 V$ D: H5 ? }
for(i = 0; i < BLOCK_SIZE; i++) {
5 Z4 E: f2 d& Z float t;5 g2 M. ]) ]4 E: @
comp -= matA[tidr] * matB[tidc];
* B6 g& h; K: x& h3 D1 O& i& n t = results - comp;
, A8 D0 H6 L4 F; R0 E4 M9 i: S comp = (t - results) + comp;* C* _4 T: b2 ]/ E
results = t;2 c$ I5 t6 `2 Y) F
}
4 N& C; a( k/ g' l9 C% j
3 M( Q4 E* ?( D __syncthreads();
! s! p8 U2 r' V f }
6 `8 d2 X4 a# o( j
) F4 G; C( ?; K* k5 r8 ` if(tidr + bidr < n && tidc + bidc < n) {
. I- b, E9 g/ L( l/ F6 ^' T; B% [3 h c[(tidr + bidr) * ldc + tidc + bidc] = results;" L9 [7 N% G+ @1 } g0 L
}
5 z) D3 V! H& W" \6 T; P7 T; w }3 T; \# N2 w% m" U# ?# o$ Y
注意到因为我们现在使用 16x16 的 threads，因此 threadIdx 变量可以取得 threadIdx.x 和 threadIdx.y，范围分别是 0 ~ 15。blockIdx.x 和 blockIdx.y 变量也是同样的情形，范围分别是 0 ~ n/16。6 T$ L6 L' A7 \+ L) U
在程序中，因为矩阵的大小不一定会是 16 的倍数，因此需要使用 if 判断式检查是否超出矩阵范围。) Z" O h% K( ^0 ~9 q2 _5 p
这个版本在 GeForce 8800GT 上的执行结果如下：8 a; w* {) U0 L- S$ f7 d1 ^
Max error: 1.19209e-007 Average error: 4.22751e-008" b! K i( \: _ }
Time used: 0.0780 (25.64 GFLOPS)1 h4 Y0 p3 ^* i" w: X& b
速度虽然提高了，但是似乎并没有达到预期中的八倍。当然，前面提到过，我们在计算时间时，把一些复制内存、配置内存的动作也计算在内，这些动作的时间并不会缩短。实际上 kernel 的运行时间，大约是 0.053 秒左右（约略相当于 38GFLOPS），比上一节的版本快了将近一倍。! {0 {. V# n2 }- D |, Z( v- J% ?
如果这一版程序已经不再限于内存带宽，那为什么没有达到预期的效率呢？这是因为这一版程序已经是限于指令周期了。除了使用 Kahan's Summation Formula 会需要更多的运算之外，程序中也有大量计算矩阵地址的乘法等等，这都会需要花费运算资源。另外，那些用来判断超出矩阵范围的 if 判断式，也会有一定的影响。7 p7 o% n8 d/ Q2 k
要把那些 if 判断式去掉，有一个方法是，在配置内存时，就配置成 16 的倍数，并在复制矩阵到显卡内存之前，先将它清为 0。如下所示：0 P3 J+ s' [& Y9 R
int newn = ((n + BLOCK_SIZE - 1) / BLOCK_SIZE) * BLOCK_SIZE;
' H* \4 [+ T6 d' p5 m# ~, |4 a% Z8 W2 F4 x; l. |, m% s' j7 a& H
cudaMallocPitch((void**) &ac, &pitch_a,
% w( Q( t8 a& V; A1 B8 r sizeof(float) * newn, newn);) [& L! \3 E) Z, H
cudaMallocPitch((void**) &bc, &pitch_b,0 R- z& \4 ]& N* e. `$ Q: Y
sizeof(float) * newn, newn);, P4 @3 v \. y! y$ l, g
cudaMallocPitch((void**) &cc, &pitch_c,
8 ^/ I+ \! z) V: X- C- [$ J8 u sizeof(float) * newn, newn);9 q* B( P* N* x4 c8 m
5 U: z; | Z8 j- Z
cudaMemset(ac, 0, pitch_a * newn);4 V1 B; h5 I; |+ S9 Z
cudaMemset(bc, 0, pitch_b * newn);8 N6 v3 Y% u; d+ U% G2 C$ {
这样一来，我们就可以把 kernel 中的 if 判断式都移除了：$ m0 y+ f) }# O+ N
__global__ static void matMultCUDA(const float* a, size_t lda,
, ~* b1 ~# I3 _1 n9 L( u e+ t const float* b, size_t ldb, float* c, size_t ldc, int n)0 l" `% b. u$ z$ P5 s# A) P
{9 s7 i" C: F4 P4 g4 n$ g4 q3 i
__shared__ float matA[BLOCK_SIZE][BLOCK_SIZE];
9 I% k, R. J% ] __shared__ float matB[BLOCK_SIZE][BLOCK_SIZE];0 b h, [. g! }1 Q+ m
const int tidc = threadIdx.x;: |* g0 @! D; G0 _% G s
const int tidr = threadIdx.y;& h) q, a6 Z6 z, g0 S
const int bidc = blockIdx.x * BLOCK_SIZE;* l; P' V7 F% j! @) K0 F
const int bidr = blockIdx.y * BLOCK_SIZE;# n( d5 o, ^0 Z( T/ C3 W* B
int i, j;( W+ b+ u2 c5 J9 Q( Z# L. d
) o/ \) \& j) A: ^! R1 {
float results = 0;% J" J# I) z) R. f& r1 R3 h: Y
float comp = 0;2 g( h+ Q+ W# ^. O" ]6 d. t& Y
( B* H$ D1 s2 R) e
for(j = 0; j < n; j += BLOCK_SIZE) {
: b$ O" ]( c: T6 M. P* {# |: K0 A matA[tidr][tidc] = a[(tidr + bidr) * lda + tidc + j];: u4 t4 u7 E Q+ V% A8 D
matB[tidr][tidc] = b[(tidr + j) * ldb + tidc + bidc];
( }1 M4 S/ g& ~9 a2 ^4 c/ ?, d( Q; X8 U" K# u
__syncthreads();' F" |3 Y" X1 B; ?3 }) R
- ~1 n; ^' {; v0 W: m for(i = 0; i < BLOCK_SIZE; i++) {, O- L$ j2 e4 r/ } @( r
float t;! x5 h8 b4 m0 I+ f. x
comp -= matA[tidr] * matB[tidc];
4 U4 }6 R2 |, X$ ~ t = results - comp;
- g1 R7 C; L+ _# v- f3 ` comp = (t - results) + comp;
6 k/ ~) \/ N0 |* S5 Z results = t;% M, y: X) s( H z* [
}/ p: L. N4 f/ u5 X3 V
|5 S! X3 w, S! m7 G2 H: [ __syncthreads();& e" T: N' g2 i; ^" a$ [" a
}) }9 N" w$ A' u" {/ i
7 @& A( q& J4 c# u' s- T- e! v c[(tidr + bidr) * ldc + tidc + bidc] = results;
; x! o ?2 ]" ^* j, [ }
8 o" f% W- p- I6 ~$ M# R5 w4 R" h这个版本的执行结果是：6 X+ ?& D( ?+ \! |
Max error: 1.19209e-007 Average error: 4.22751e-008
( d3 l3 R! h8 Q# ]2 Q Time used: 0.0780 (25.64 GFLOPS)
" O$ R: Z) M: _1 M7 C似乎没有改善。不过，实际上 kernel 的运行时间已经减少到 0.042 秒（约略相当于 48GFLOPS）。9 F: g( H# O4 u7 O; d0 d* j
结论
有些读者可能会想，如果把 block 再变得更大（例如 32x32）是否会有帮助呢？当然，由于最后的程序已经不再是受限于内存带宽（在 0.042 秒内存取 500MB 的数据约相当于 12GB/s 的带宽），所以把 block 再加大并不会有帮助了。而且，由于一个 block 内的 thread 数目最多只能到 512 个，将 block 变大也会造成很多额外负担。而且 shared memory 的大小也有限制（GeForce 8800GT 的 shared memory 大小限制是 16384 bytes），所以也不能任意增加 block 的大小。. e0 W; M# x4 Y
最后一版程序的完整档案可以从这里下载。 |
|
发表于 2008-9-9 20:35:00
