CUDA編程－（2）其實寫個矩陣相乘並不是那么難

 程序代碼及圖解析：

#include <iostream>
#include "book.h"
__global__ void add( int a, int b, int *c ) {
*c = a + b;
}
int main( void ) {
int c;
int *dev_c;
HANDLE_ERROR( cudaMalloc( (void**)&dev_c, sizeof(int) ) );
add<<<1,1>>>( 2, 7, dev_c );
HANDLE_ERROR( cudaMemcpy( &c,
dev_c,
sizeof(int),
cudaMemcpyDeviceToHost ) );
printf( "2 + 7 = %d\n", c );
cudaFree( dev_c );
return 0;
}

　　

函數原型：__host__cudaError_t cudaMemcpy (void *dst, const void *src, size_t count, cudaMemcpyKind kind)

作用：在設備端和主機端拷貝數據。

參數：dst 目的地址 src 源地址 count 拷貝字節大小kind 傳輸的類型

返回值：

cudaSuccess, cudaErrorInvalidValue, cudaErrorInvalidDevicePointer, cudaErrorInvalidMemcpyDirection

說明:

從源地址拷貝設定數量的字節數至目的地址，kind類型有四種，分別為：

cudaMemcpyHostToHost, cudaMemcpyHostToDevice,  cudaMemcpyDeviceToHost, cudaMemcpyDeviceToDevice，

通過指定方向進行拷貝。存儲器區域不可重疊。如若產生未定義拷貝方向的行為，dst和src將不匹配。

 正文

前面的圖是最簡單的一個CUDA程序，它引出了Grid Block Thread概念。很多threads組成1維，2維or3維的thread block. 為了標記thread在block中的位置（index），我們可以用上面講的threadIdx。threadIdx是一個維度<=3的vector。還可以用thread index（一個標量）表示這個位置。

thread的index與threadIdx的關系：

	Thread index
1	T
2	T.x + T.y * Dx
3	T.x+T.y*Dx+z*Dx*Dy

其中T表示變量threadIdx。（Dx, Dy, Dz）為block的size（每一維有多少threads）。

因為一個block內的所有threads會在同一處理器內核上共享內存資源，所以block內有多少threads是有限制的。目前GPU限制每個 block最多有1024個threads。但是一個kernel可以在多個相同shape的block上執行，效果等效於在一個有N*#thread per block個thread的block上執行。

Block又被組織成grid。同樣，grid中block也可以被組織成1維，2維or3維。一個grid中的block數量由系統中處理器個數或待處理的數據量決定。（來自這里)

 下圖中描述了Thread、Block、Grid內存的訪問機制。

每個thread有自己的local-memory。每一個block有自己的共享內存、grid和grid之間可以同時訪問全局內存。這里要注意：block和block之間不能訪問同一個共享內存，他們只能訪問自己的共享內存。

cudaGetDeviceCount( &count )查詢服務器的CUDA信息.

#include <stdio.h>
#include <cuda_runtime.h>
int main()
{
        int deviceCount;
        cudaGetDeviceCount(&deviceCount);
        int device;
        for(device = 0; device < deviceCount; ++device)
        {
                cudaDeviceProp deviceProp;
                cudaGetDeviceProperties(&deviceProp,device);
                printf("Device %d has compute capability %d.%d.\n",device,deviceProp.major,deviceProp.minor);
        }
}

　　結果：

struct cudaDeviceProp {
    char name[256]; //識別設備的ASCII字符串（例如，“GeForce GTX 280”）
    size_t totalGlobalMem; //全局內存大小
    size_t sharedMemPerBlock; //每個block內共享內存的大小
    int regsPerBlock; //每個block32位寄存器的個數
    int warpSize; // warp大小
    size_t memPitch; //內存中允許的最大間距字節數
    int maxThreadsPerBlock; //每個Block中最大的線程數是多少
    int maxThreadsDim[3]; // 一個塊中每個維度的最大線程數
    int maxGridSize[3]; //一個網格的每個維度的塊數量
    size_t totalConstMem; //可用恆定內存量
    int major; //該設備計算能力的主要修訂版號
    int minor; //設備計算能力的小修訂版本號
    int clockRate; //時鍾速率
    size_t textureAlignment; //該設備對紋理對齊的要求
    int deviceOverlap; //一個布爾值，表示該裝置是否能夠同時進行cudamemcpy()和內核執行
    int multiProcessorCount; //設備上的處理器的數量
    int kernelExecTimeoutEnabled; //一個布爾值，該值表示在該設備上執行的內核是否有運行時的限制
    int integrated; //返回一個布爾值，表示設備是否是一個集成的GPU（即部分的芯片組、沒有獨立顯卡等）
    int canMapHostMemory; //表示設備是否可以映射到CUDA設備主機內存地址空間的布爾值
    int computeMode; //一個值，該值表示該設備的計算模式：默認值，專有的，或禁止的
    int maxTexture1D; //一維紋理內存最大值
    int maxTexture2D[2]; //二維紋理內存最大值
    int maxTexture3D[3]; //三維紋理內存最大值
    int maxTexture2DArray[3]; //二維紋理陣列支持的最大尺寸
    int concurrentKernels; //一個布爾值，該值表示該設備是否支持在同一上下文中同時執行多個內核
｝

 

矩陣相乘也非常簡單，難在如何在這個基礎上提高速率。比如：引入sharememory。

代碼：

#include <stdio.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <time.h>
#include <stdlib.h>

__global__ void MatrixMuiOnDevice(int *M,int *N, int *P, int width)
{       
        int x = threadIdx.x; 
        int y = threadIdx.y; //獲取該線程的位置
        
        float Pervalue = 0;
        
        for (int i = 0; i < width; i++)
        {       
                float Mdlement = M[y * width + i];
                float Ndlement = N[width * i + x];
                
                Pervalue += Mdlement * Ndlement;
        }
        
        P[y * width + x] = Pervalue;
}
int main()
{       
        int a[30][30],b[30][30],c[30][30];
        int *M, *N, *P;
        int width = 30;
        int NUM = 900;
        dim3 dimBlock(30,30);
        cudaEvent_t start,stop;
        float elapsedTime;
        cudaEventCreate(&start);
        cudaEventCreate(&stop);


        cudaMalloc((void**)&M, 900*sizeof(int));
        cudaMalloc((void**)&N, 900*sizeof(int));
        cudaMalloc((void**)&P, 900*sizeof(int));
        //初始化
        for(int i = 0; i < 30; i++)
                for(int j = 0; j < 30; j++)
                {
                        a[i][j] = 2;
                        b[i][j] = 3;
                }

        cudaMemcpy(M,a,NUM*sizeof(int),cudaMemcpyHostToDevice);
        cudaMemcpy(N,b,NUM*sizeof(int),cudaMemcpyHostToDevice);
        cudaMemcpy(c,P,NUM*sizeof(int),cudaMemcpyDeviceToHost);
        cudaEventRecord(start,0);
        MatrixMuiOnDevice<<<1,dimBlock>>>(M,N,P,width);
        cudaThreadSynchronize();
        cudaEventRecord(stop,0);
        cudaEventSynchronize(stop);
        cudaEventElapsedTime(&elapsedTime,start,stop);

        printf("%f\n",elapsedTime);
        for(int i = 0; i < 30; i++)
                for(int j = 0; j < 30; j++)
                {
                        printf("%d \n",c[i][j]);
                }

        cudaFree(M);
        cudaFree(N);
        cudaFree(P);
        return 0;
}   

　　share memory 改進。加入同步機制 __syncthreads(),即 等待之前的所有線程執行完畢后再接下去執行。

#include <stdio.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <time.h>
#include <stdlib.h>

#define TILE_WIDTH 25

__global__ void MatrixMuiOnDevice(int *M,int *N, int *P, int width)
{
       
        __shared__ float Mds[TILE_WIDTH][TILE_WIDTH];
        __shared__ float Nds[TILE_WIDTH][TILE_WIDTH];

    
        int bx = blockIdx.x;
        int by = blockIdx.y;
        int tx = threadIdx.x;
        int ty = threadIdx.y;
        int Col = bx * TILE_WIDTH + tx;
        int Row = by * TILE_WIDTH + ty; //獲取該線程的位置
    
        int Pervalue = 0;
        
        for (int i = 0; i < width / TILE_WIDTH; i++)
        {       
            Mds[ty][tx] = Md[Row * width+(i * TILE_WIDTH + tx)];
            Nds[ty][tx] = Nd[Col + (i * TILE_WIDTH + ty) * width];
            __syncthreads();

            for (int  k = 0; k < width / TILE_WIDTH; k++)
                Pervalue += Mds[ty][k] * Nds[k][tx];
            __syncthreads();         
        }

        P[Row * width + Col] = Pervalue;
}

int main()
{
        int WID = 100;
        int a[WID][WID],b[WID][WID],c[WID][WID];
        int *M, *N, *P;
        int width = WID / 4 ;;
        int NUM = WID*WID;
        dim3 dimGrid(WID/width,WID/width);
        dim3 dimBlock(width,width);
        cudaEvent_t start,stop;
        float elapsedTime;
        cudaEventCreate(&start);
        cudaEventCreate(&stop);


        cudaMalloc((void**)&M, NUM*sizeof(int));
        cudaMalloc((void**)&N, NUM*sizeof(int));
        cudaMalloc((void**)&P, NUM*sizeof(int));
        //初始化
        for(int i = 0; i < 100; i++)
                for(int j = 0; j < 100; j++)
                {
                    a[i][j] = 2;
                    b[i][j] = 3;
                }

        cudaMemcpy(M,a,NUM*sizeof(int),cudaMemcpyHostToDevice);
        cudaMemcpy(N,b,NUM*sizeof(int),cudaMemcpyHostToDevice);
        cudaMemcpy(c,P,NUM*sizeof(int),cudaMemcpyDeviceToHost);
        cudaEventRecord(start,0);
        MatrixMuiOnDevice<<<dim,dimBlock>>>(M,N,P,width);
        cudaThreadSynchronize();
        cudaEventRecord(stop,0);
        cudaEventSynchronize(stop);
        cudaEventElapsedTime(&elapsedTime,start,stop);

        printf("%f\n",elapsedTime);

        cudaFree(M);
        cudaFree(N);
        cudaFree(P);
        return 0;
}

　　

小結

第一個執行時間：

share memory執行時間：

 注意，核函數內不是所有線程一起進去執行，這個概念模糊不清。我們需要理解成，所有的線程並行執行核函數里面的程序，即每一個線程都會執行該函數，所有線程執行完，即結束。這個簡單的概念，我一開始想了很久。

注：轉載請注明出處。