共享數(shù)據(jù)及課件gpu程序設(shè)計(jì)初步

GPU CUDA程序設(shè)計(jì)初步第2講

CUDA程序優(yōu)化

1計(jì)算流體力學(xué)與湍流培訓(xùn)系列教程李新亮中國科學(xué)院力學(xué)研究所

2021-07參考書目：CUDA官方網(wǎng)站資料：攀哲勇，《CUDA編程基礎(chǔ)與實(shí)踐》，清華大學(xué)出版社，2020ShaneCook著，蘇統(tǒng)華等譯《CUDA并行程序設(shè)計(jì)GPU編程指南》，機(jī)械工業(yè)出版社，2020示例: 計(jì)算兩個(gè)矩陣的乘積??

=

???

???????

=

????????????已知矩陣A和B均為N*N矩陣，

N=1024, 利用GPU計(jì)算矩陣C=A*B, 并分析結(jié)果的正確性及計(jì)算效率分析：啟動(dòng)N*N個(gè)線程，每個(gè)線程計(jì)算一個(gè)Cij

點(diǎn)；在計(jì)算Cij

過程中，Bkj

涉及非連續(xù)數(shù)據(jù)的訪問，影響效率引入矩陣BT, 使得數(shù)據(jù)連續(xù)訪問?????

=

??????????????????

=

????????????#include<stdio.h>#include<stdlib.h>#include<time.h>void

matrix_mul_cpu(

int

n,

float

*a,

float*b,

float*c){for

(int

i

=

0;

i

<

n;

++i)for

(int

j=

0;

j

<

n;

++j){float

s

=

0.0;for

(int

k

=0;

k

<

n;

++k)s

+=

a[i*n

+

k]

*

b[k*n

+

j];

//A(i,k)*B(k,j)c[i*n

+

j]

=

s;}}

global

void

matrix_mul_gpu(int

n,

float

*a,

float*b,

float*c){const

int

bid=

blockIdx.x;const

int

tid

=threadIdx.x;float

s

=0.0;for

(int

k

=

0;

k

<

n;

++k)s

+=

a[bid*n

+

k]

*

b[k*n

+

tid];

//A(i,k)*B(k,j)c[bid*n+tid]

=

s;}CPU計(jì)算Cij（計(jì)算N*N個(gè),全部元素）GPU計(jì)算Cij（每個(gè)線程只計(jì)算1個(gè)元素）B[]數(shù)組非連續(xù)訪問

global

voidtranspose(intn,float*b,float*bt){bt[blockIdx.x*n

+

threadIdx.x]

=

b[threadIdx.x*n

+

blockIdx.x];}

global

void

matrix_mul_gpu2(int

n,

float

*a,

float*b,

float*c){const

int

bid

=

blockIdx.x;constint

tid

=

threadIdx.x;float

s

=0.0;for

(int

k

=

0;

k

<

n;++k)//A*BT連續(xù)訪問s+=

a[bid*n

+

k]

*

b[tid*n

+k];c[bid*n

+

tid]

=

s;}GPU計(jì)算CijA,B矩陣均連續(xù)訪問（需要首先將B矩形轉(zhuǎn)置）矩陣轉(zhuǎn)置void

check_data(int

n,

float

*c,

float

*c1){float

s

=0.0;for

(int

i

=

0;

i

<

n*n;

++i)s

+=

fabs(c[i]

-

c1[i]);printf("Total

error

is

%f

\n",

s);}測試計(jì)算是否正確intmain(){const

int

n

=

1024,

M

=

n*n*

sizeof(float);float

*a

=

(float*)malloc(M);float

*b

=

(float*)malloc(M);float

*c

=

(float*)

malloc(M);float

*c1

=(float*)malloc(M);float

*d_a,

*d_b,

*d_c,

*d_bT;for

(int

i

=0;

i

<n*n;++i){a[i]

=

(float)

(rand()%100);b[i]

=

(float)

(rand()%100);}clock_t

time1,

time2,time3,time4;time1

=clock();matrix_mul_cpu(n,

a,

b,

c);time2

=clock();printf("Time

for

CPUrun

is:

%f

seconds

\n",

(double)(time2

-

time1)

/

CLOCKS_PER_SEC);cudaMalloc(&d_a,

M);cudaMalloc(&d_b,

M);cudaMalloc(&d_c,

M);cudaMemcpy(d_a,

a,

M,

cudaMemcpyHostToDevice);cudaMemcpy(d_b,

b,

M,

cudaMemcpyHostToDevice);Time3=clock();賦初值CPU計(jì)算C=A*BGPU計(jì)算

C=A*B開辟數(shù)組A,B，C；matrix_mul_gpu<<<n,n>>>(n,d_a,d_b,d_c);//計(jì)算C=A*B,每個(gè)線程上計(jì)傳算一至個(gè)點(diǎn)GPU端cudaDeviceSynchronize();cudaMemcpy(c1,d_c,

M,

cudaMemcpyDeviceToHost);time4

=clock();printf("Time

for

GPU

run

is:

%f

seconds

\n",

(double)(time4

-

time3)

/

CLOCKS_PER_SEC);check_data(n,c,

c1);cudaFree(d_a);cudaFree(d_b);cudaFree(d_c);cudaFree(d_bT);//

by

using

transposetime1

=

clock();cudaMalloc(&d_a,M);cudaMalloc(&d_b,M);cudaMalloc(&d_c,

M);cudaMalloc(&d_bT,

M);cudaMemcpy(d_a,

a,

M,

cudaMemcpyHostToDevice);cudaMemcpy(d_b,

b,

M,

cudaMemcpyHostToDevice);transpose

<<<n,

n

>>>

(n,

d_b,

d_bT);cudaDeviceSynchronize();matrix_mul_gpu2

<<<n,

n

>>>

(n,

d_a,

d_bT,

d_c);cudaMemcpy(c1,

d_c,

M,

cudaMemcpyDeviceToHost);time2

=

clock();printf("Time

for

GPU

run

is:

%f

seconds

\n",

(double)(time2-

time1)

/CLOCKS_PER_SEC);check_data(n,

c,

c1);free(a)