CUDA C 任意矩陣相乘
阿新 • • 發佈：2019-02-10
好久沒來了，繼續學習！留下腳印，留下回憶！目前在深圳實習！
/***************************** *A=M*I,B=I*N *求矩陣C=A*B,無論怎麼看，求任意矩陣內積，即使只用一個grid，在邊界上也存在很多問題。 *最主要的就是在邊界上，如果block尺寸大於劃分的資料小矩陣的尺寸，那麼只有部分執行緒使用，這倒沒有問題， *關鍵在於，對於部分執行緒塊而言，執行緒塊與矩陣資料塊大小不同，就不能做到一個執行緒取一個數，邏輯就有問題了 *NOTE:考慮這麼多，想想最簡單的解決辦法，還是將資料矩陣補齊為能整除32*32塊 *memset()函式因為是按位元組填充，所以只適合清零 ******************************/ #include<stdio.h> #include<stdlib.h> #include<iostream> #include<cuda_runtime.h> #include<device_launch_parameters.h> using namespace std; #define M 33 #define N 33 #define I 33 #define block_size 32 #define size_a M*I*sizeof(float) #define size_b I*N*sizeof(float) #define size_c M*N*sizeof(float) __global__ void matrix_kernel(float *d_a,float *d_b,float *d_c,int wA,int wB);//wA,wB是A,B矩陣的寬度，也就是列數 void showResult(float *c,int h,int Kw,int w); void init(float *a,int h,int Kw,int w,int value);//Kw為擴充後矩陣寬度，而w為實際矩陣寬度 //考慮將矩陣補齊，就是行列都變成32*32塊的整數倍,函式功能是補齊並初始化,還要分配空間 float* supply(int h,int w,int value);//引數為行，寬，就可以得到擴充後的大小,value為需要初始化的值 int main() { float *h_a,*h_b,*h_c;//宣告主機指標 float *d_a,*d_b,*d_c;//宣告裝置指標 //在考慮補齊的時候,實際行數列數都發生了變化,記住最後的矩陣大小為row_c*col_c,但是實際矩陣大小為M*N即可 int row_a=32*((M+32-1)/32);//補齊後a的行列，但是從[32*M/32,M)為初始化數值，而[M,row_a)則應初始化為零 int col_a=32*((I+32-1)/32);//[I,col_a)初始化為零 cout<<"a矩陣擴充後的行列分別為："<<row_a<<" "<<col_a<<" "<<(M+31)/32<<" "<<32*((M+32-1)/32)<<endl; int row_b=col_a;//[I,row_b)行，初始化為0 int col_b=32*((N+32-1)/32);//[N,col_b）列初始化為0 int row_c=row_a; int col_c=col_b; cout<<"c矩陣擴充後的行列分別為："<<row_c<<" "<<col_c<<" "<<(M+31)/32<<" "<<32*((M+32-1)/32)<<endl; cout<<"b矩陣擴充後的行列分別為："<<row_b<<" "<<col_b<<" "<<(M+31)/32<<" "<<32*((M+32-1)/32)<<endl; int size_a_a=row_a*col_a*sizeof(float);//補齊之後的矩陣所佔用空間 int size_b_b=row_b*col_b*sizeof(float); int size_c_c=row_c*col_c*sizeof(float); cout<<"a,b,c矩陣擴充後所佔用空間："<<size_a_a<<" "<<size_b_b<<" "<<size_c_c<<endl; h_a=supply(row_a,col_a,1);//分配空間並初始化 //顯示擴充後a矩陣 showResult(h_a,row_a,col_a,col_a); h_b=supply(row_b,col_b,2); showResult(h_b,row_b,col_b,col_b); h_c=supply(row_c,col_c,0); showResult(h_c,row_c,col_c,col_c); /* h_a=(float*)malloc(size_a);init(h_a,M,I,1);//主機上分配空間並初始化 h_b=(float*)malloc(size_b);init(h_b,I,N,2);showResult(h_a,M,I);cout<<"a矩陣輸出完成！"<<endl;showResult(h_b,I,N);cout<<"b矩陣輸出完成"<<endl; h_c=(float*)malloc(size_c);memset(h_c,0,size_c); */ if(h_a==NULL||h_b==NULL||h_c==NULL) { fprintf(stderr,"分配主機記憶體失敗!\n"); exit(EXIT_FAILURE); } //裝置上分配空間並檢查 cudaError_t err=cudaSuccess;//設定成功標誌碼 err=cudaMalloc((void**)&d_a,size_a_a); cudaMalloc((void**)&d_b,size_b_b); cudaMalloc((void**)&d_c,size_c_c);cudaMemset(d_c,0,size_c_c); /* err=cudaMalloc((void**)&d_a,size_a); cudaMalloc((void**)&d_b,size_b); cudaMalloc((void**)&d_c,size_c);cudaMemset(d_c,0,size_c);//只對c初始化為0 */ if(err!=cudaSuccess) { fprintf(stderr,"分配裝置記憶體a失敗！\n"); exit(EXIT_FAILURE); } //加入時間計時程式碼 cudaEvent_t start,stop;//宣告事件 cudaEventCreate(&start);//建立事件 cudaEventCreate(&stop); cudaEventRecord(start,0);//sample使用NULL流，應該和這裡一樣反正都是預設流的意思 //#################GPU工作部分######################矩陣大小變化了，後面這些東西都要變 cudaMemcpy(d_a,h_a,size_a_a,cudaMemcpyHostToDevice);//主機向device拷貝資料，對於錯誤的檢查，就不加了，太麻煩 cudaMemcpy(d_b,h_b,size_b_b,cudaMemcpyHostToDevice); dim3 dimsA(col_a,row_a);//這裡始終是第一個引數為列，第二個引數為行 dim3 dimsB(col_b,row_b); //dim3 grids((N+block_size-1)/block_size,(M+block_size-1)/block_size);//按照我的補齊思路其實也不用這樣寫，因為始終是齊的 //這裡一定要記住了，每次都犯錯誤，對於cuda中的dim3寫法，三個元素（x,y,z）對應的就是列，行，高， //所以說在這裡寫個grid（32,64）表示32列，64行 dim3 grids(col_c/block_size,row_c/block_size); dim3 threads(block_size,block_size); //呼叫kernel函式計算內積 matrix_kernel<<<grids,threads>>>(d_a,d_b,d_c,dimsA.x,dimsB.x); err=cudaMemcpy(h_c,d_c,size_c_c,cudaMemcpyDeviceToHost);//將結果資料拷貝回主機記憶體 if(err!=cudaSuccess) { fprintf(stderr,"拷貝資料回主機失敗！\n"); exit(EXIT_FAILURE); } cudaEventRecord(stop,0);//將記錄結束時間放入GPU工作佇列中 cudaEventSynchronize(stop);//CPU等待GPU完成所有工作後，再繼續執行 float totalTime=0.0f;//定義儲存時間的變數，使用float型 cudaEventElapsedTime(&totalTime,start,stop); //至此所有工作全部完成，剩下就是輸出和釋放了 printf("GPU計算時間為：%f\n",totalTime); showResult(h_c,M,col_c,N);//顯示輸出結果 cout<<M<<'\t'<<N<<endl; showResult(h_c,row_c,col_c,col_c);//顯示全部的大矩陣 cout<<row_c<<'\t'<<col_c<<endl; //釋放資源 free(h_a);free(h_b);free(h_c); cudaFree(d_a);cudaFree(d_b);cudaFree(d_c); cudaEventDestroy(start);cudaEventDestroy(stop); return 0; } //這裡核函式的執行緒排程為block（32,32） grid（20，10），即grid是20列，10行 __global__ void matrix_kernel(float *d_a,float *d_b,float *d_c,int wA,int wB) { //block index int bx=blockIdx.x,by=blockIdx.y; //thread index int tx=threadIdx.x,ty=threadIdx.y; //這裡只是考慮的是對齊的情況，對於每一個block而言，先得到block的初始指標位置aBegin實際上表示的是一種位置的偏移量，指標加上偏移量就得到值，本身不使用指標更方便 int aBegin=by*block_size*wA;//其實就是by*32*I，I就是A的寬度，這就很顯然了 int aEnd=aBegin+wA-1;//其實aBegin和aEnd就是一行的首尾位置 int aStep=block_size;//這裡注意a的步幅是在for迴圈中的跨度，因為每次只是載入一個塊，所以整行要分多個塊，就是這塊之間的距離 int bBegin=bx*block_size;//畫圖看的話，很清楚。分清楚，對於不同的block有不同的aBegin和bBegin，與for迴圈中跨度區別 int bStep=block_size*wB; float Csub=0;//每個執行緒都會計算出一個行列的相乘值出來，放在Csub中 for(int a=aBegin,b=bBegin;a<aEnd;a+=aStep,b+=bStep) { __shared__ float as[block_size][block_size]; __shared__ float bs[block_size][block_size]; as[ty][tx]=d_a[a+ty*wA+tx];//這裡其實本質上這樣看：a+tx其實就是在多少列上，這樣看也行，怎麼看都行，或者看為a為起始位置，加上ty*wA其實就是再下移ty行的位置，在加上tx就可以了 bs[ty][tx]=d_b[b+ty*wB+tx];//bs與as取值的規律本質上是一樣的。僅僅只是矩陣不同，而導致寬度不同而已 __syncthreads();//保證資料都載入shared //第二個for迴圈計算乘累加 for(int i=0;i<block_size;i++) { Csub+=as[ty][i]*bs[i][tx]; } __syncthreads();//做到這裡就表示一個分塊的小矩陣計算完成，再進入下一個分塊的小矩陣 } //當所有迴圈做完後，每個執行緒都得到一個Csub值，寫入d_c相應位置即可 int c=wB*by*block_size+bx*block_size;//c在這裡其實表示的是每個block的起始位置，要得到每個thread的在grid中的位置，還要加上在block中的相對位置 d_c[c+ty*wB+tx]=Csub; } //初始化資料函式,我採取的策略是對矩陣進行兩次賦值，每次賦值的行列不相同， void init(float *a,int h,int Kw,int w,int value)//w始終是擴充後矩陣的寬，也就是列數，h為高，也就是行數, { for(int i=0;i<h;i++) { for(int j=0;j<w;j++) a[i*Kw+j]=value; } } //顯示結果函式 void showResult(float *c,int h,int Kw,int w)//h高度就是行，w寬度就是列 { for(int i=0;i<h;i++) { for(int j=0;j<w;j++) { cout<<c[i*Kw+j]<<'\t'; } cout<<endl; } } //補齊矩陣函式 float* supply(int h,int w,int value) { float *p=(float*)malloc(h*w*sizeof(float)); init(p,h,w,w,0);//先將全矩陣初始化為0, init(p,M,w,I,value);//再講實際需要運算的部分初始化為value，這裡需要擴充矩陣的寬度，和需要實際賦值的矩陣寬度 return p; }