CUDA實現矩陣相加的平行計算
（一）目的
熟悉基本的CUDA程式架構以及如何呼叫相應的API進行CUDA程式設計
（二）內容
完成矩陣相加的並行程式的實現（不用share memory實現）
要求：
- 實現2個矩陣（32*32）的相加，M矩陣的初始值全為2，N矩陣的初始值全為5。同時用CPU程式碼實現，比較兩個程式碼的執行時間
- 實現2個矩陣（1024*1024）的相加，M矩陣的初始值全為2，N矩陣的初始值全為5。同時用CPU程式碼實現，比較兩個程式碼的執行時間
- 實現2個矩陣（5000*5000）的相加，M矩陣的初始值全為2，N矩陣的初始值全為5。同時用CPU程式碼實現，比較兩個程式碼的執行時間
實驗步驟一 軟體設計分析:
一．資料型別：
根據實驗任務，M矩陣的初始值為2，N矩陣的初始值為5，相加後得到的矩陣應該為值全為7的矩陣。因此，把實驗的資料型別定義為short int型別比較合適。
二．對於方陣的儲存，一共有2種方式：
矩陣在記憶體中的儲存按照行列先後可以分為兩種方式，一種是行優先的儲存方式，一種是列優先的方式。
這兩種儲存方式是在訪問對應位置的資料的時候有很大的差別。在CULA內部，矩陣預設是按列優先存放的，如果要使用CULA device 函式，就必須考慮儲存方式的問題，有的時候可能需要我們對儲存方式進行轉換。但是無論哪種儲存方式，最終在記憶體中是順序儲存的。
三．GPU程式的block和threads的相關設定：
我們cuda實驗平臺的每一個Grid可以按照一維或者二維的方式組織，每一個Block可以按照一維，二維或者三維的方式進行儲存，每一個block最多隻能有1536個執行緒。對於上面三個任務GPU程式的block和threads的設定是不一樣的。
(1)第一個任務，實現2個矩陣（32*32）的相加的時候，需要32*32=1024（<1536）個執行緒，這時的block和threads的設定為：
DimGrid(1,1)，DimBlock(32,32);
即Grid裡面的Block是按照1x1的方式組織，每一個Block裡面的執行緒按照32x32組織。
(2)第二個任務，實現2個矩陣（1024*1024）的相加的時候一個Block肯定是不夠的，這時候我的的block和threads的設定為：
DimGrid(32,32)，DimBlock(32,32);
即Grid裡面的Block是按照32x32的方式組織，每一個Block裡面的執行緒按照32x32組織。
(3)第三個任務，實現2個矩陣（5000*5000）的相加的時候我的的block和threads的設定為：
DimGrid(157,157)，DimBlock(32,32);
即Grid裡面的Block是按照157x157的方式組織，每一個Block裡面的執行緒按照32x32組織。
實驗步驟二 實驗裝置:
本地裝置：PC機+Windows10作業系統
Putty遠端連線工具
遠端裝置：NVIDIA-SMI 352.79
Driver Version: 352.79
實驗步驟三 CPU計算程式碼:
#include<stdio.h>
#include<stdlib.h>
#include<time.h>
const int DIM = 5000;
short int M[DIM][DIM];
short int N[DIM][DIM];
short int C[DIM][DIM];
void main()
{
clock_t start, end;
//生成M矩陣
for (int i = 0; i < DIM; i++)
for (int j = 0; j < DIM; j++)
M[i][j] = 2;
//生成N矩陣
for (int i = 0; i < DIM; i++)
for (int j = 0; j < DIM; j++)
N[i][j] = 5;
//求矩陣的和
start = clock();//計時開始
for (int i = 0; i < DIM; i++)
for (int j = 0; j < DIM; j++)
C[i][j] = M[i][j] + N[i][j];
end = clock();//計時結束
//列印結果
printf("矩陣M+矩陣N\n");
for (int i = 0; i < DIM; i++)
{
for (int j = 0; j < DIM; j++)
{
printf("%d ", C[i][j]);
}
printf("\n");
}
float time1 = (float)(end - start) / CLOCKS_PER_SEC;
printf("執行時間為：%f\n", time1);
system("pause");
}
實驗步驟四 GPU計算程式碼:
#include<cuda_runtime_api.h>
#include<device_launch_parameters.h>
#include<stdio.h>
#include<time.h>
static const int M = 32;
static const int N = 32;
__global__ void add(int a[M][N], int b[M][N], int c[M][N])
{
int i = threadIdx.x + blockIdx.x * blockDim.x;
int j = threadIdx.y + blockIdx.y * blockDim.y;
if (i < M && j < N)
{
c[i][j] = a[i][j] + b[i][j];
}
}
int main()
{
clock_t start,end;
int (*a_h)[N] = new int[M][N];
int (*b_h)[N] = new int[M][N];
int (*c_h)[N] = new int[M][N];
for(int i = 0;i < M;i++)
for(int j = 0;j < N;j++)
a_h[i][j] = 2;
for(int i = 0;i < M;i++)
for(int j = 0;j < N;j++)
b_h[i][j] = 5;
int (*a_d)[N];
int (*b_d)[N];
int (*c_d)[N];
cudaMalloc((void **)&a_d, sizeof(int)* M*N);
cudaMalloc((void **)&b_d, sizeof(int)* M*N);
cudaMalloc((void **)&c_d, sizeof(int)* M*N);
start=clock();//開始計時
cudaMemcpy(a_d, a_h, sizeof(int)* M*N, cudaMemcpyHostToDevice);
cudaMemcpy(b_d, b_h, sizeof(int)* M*N, cudaMemcpyHostToDevice);
dim3 DimGrid(1, 1);
dim3 DimBlock(32, 32);
add <<<DimGrid, DimBlock>>>(a_d, b_d, c_d);
cudaMemcpy(c_h, c_d, sizeof(int)* M*N, cudaMemcpyDeviceToHost);
end=clock();//計時結束
for (int i = 0; i < M; i++)
{
for (int j = 0; j < N; j++)
printf("%10d", c_h[i][j]);
printf("\n");
}
float time1=(float)(end-start)/CLOCKS_PER_SEC;
printf("執行時間為：%f\n",time1);
return 0;
}
實驗步驟五 觀察輸出結果:
- 第一個任務，實現2個矩陣（32*32）的相加的時候：
- CPU輸出結果：
-
- GPU輸出結果：
- 第二個任務，實現2個矩陣（1024*1024）的相加的時候：
- CPU輸出結果：
-
- GPU輸出結果：
- 第三個任務，實現2個矩陣（5000*5000）的相加的時候：
- CPU輸出結果：
-
- GPU輸出結果：
實驗結論：
cpu程式計算所需時間：
第一個任務， cpu程式計算所需時間:0.000000s
第二個任務， cpu程式計算所需時間:0.009000s
第三個任務， cpu程式計算所需時間:0.117000s
gpu程式計算所需時間：
第一個任務， gpu程式計算所需時間: 0.000000s
第二個任務， gpu程式計算所需時間:0.0100000s
第三個任務， gpu程式計算所需時間:0.1600000s
總結：平行計算是一種以空間換時間的方式來提高資料的執行效率。按照實驗預期，計算相同時間複雜度的程式，GPU在時間效率上應該優於CPU，但是由我們的實驗來看，三個任務中，GPU都沒有在時間效率上超過CPU，本人認為主要的原因是我在編寫程式的時候把host端與device端之間的資料傳輸所花費的時間也計算在內，而在資料的傳輸還受到網路延遲的影響。這使得計算較小資料量的時候，GPU體現不出明顯的優勢，而在大資料量平行計算的時候，GPU將會大大節省我們的計算時間。