由矩阵运算理解CUDA kernel函数的运行

一、kernel函数介绍

kernel在cuda中指的是一个函数，当一个kernel被调用的时候，gpu会同时启动很多个线程来执行这一个kernel，这样就实现了并行化；每个线程执行这一kernel将通过线程号来对应输入数据的下标，这样保证每个thread执行的kernel一样，但是处理的数据不一样。

核函数以下几个的前缀分别代表：

__global__:在GPU上执行，可以在CPU上被调用，也可以在GPU上被调用
__device__:在GPU上执行，只能在GPU上被调用
__host__:在CPU上执行，只能在CPU上被调用
一个kernel在cuda中可以这么定义：

二、矩阵的加法运算

1、一维矩阵的加法运算

__global__ void VecAdd(float* A,float* B, float C){
    int i = threadIdx.x;
    C[i] = A[i] + B[i]
}

// 调用
VecAddd<<<1,N>>>(A,B,C);

上面代码中调用的时候通过指定<<<1,N>>>前者1表明使用1个block，每个block有N个threads，这样在函数调用内部通过取内置变量threadIdx来得到thread的索引，以此对应数组下标，实现并行化。

2、二维矩阵的加法运算

上面的例子中只是针对thredIdx是一维的情况，实际上，为了方便，我们可以直接将threadIdx变成多维的，以方便我们对矩阵或者更高维的数据的处理，这跟C语言里的二维数组等概念是一样的，就是用x和y分别代表一个维度的下标。

__global__ void MatAdd(float A[N][N], float B[N][N], float C[N][N]){
    int x = threadIdx.x;
    int y = threadIdx.y;
    C[x][y] = A[x][y] + B[x][y];
}

int main(){
    int num_blocks = 1;
    dim3 threads_per_block(N,N);
    MatAdd<<<num_blocks,threads_per_block>>>(A,B,C);
}

这个例子是相对很容易理解的了。需要知道的是，每个block的threads数目是有最大限制的，不能超过一定的值，书上说现在的GPU每个block最多是1024个thread。所以上面的例子其实并不能处理很大的矩阵，那么要处理很大的矩阵的时候该怎么办呢？

可以用上block：

___global___ void MatAdd(float A[N][N], float B[N][N], float C[N][N]){
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int j = blockIdx.y * blockDim.y + threadIdx.y;
    C[i][j] = A[i][j] + B[i][j];
}
int main(){
    dim3 threads_per_block(16,16);
    dim3 num_blocks(N/num_blocks.x,N/num_blocks.y);
    MatAdd(A,B,C);
}

这个例子充分利用了block数量，这样每个block是16*16也就是256个threads，假设N = 32，那么num_blocks = (2,2)也就是用了4个block。上面blockIdx.x取值是[0,1]，blockDim.x取值是16,threadIdx.x取值是[0,15].
 

3、一个完整的矩阵加法案例

之前的代码没有考虑A、B、C的内存申请，这一节讨论的就是如何把host的内存(cpu data)转到device里(gpu data)。

如果要申请一块连续的GPU内存空间，直接可以通过调用cudaMalloc()来实现，调用cudaFree()释放空间，那么把host上内存copy到gpu上是通过调用cudaMemcpy()来实现。
 

__global__ void VecAdd(float* A, float* B, float* C){
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    // 根据下面的分析，i是可能大于等于N的，
    //当N>threadsPerBlock且为threadsPerBlock的非整数倍的时候，
    //这时候会引入多一个block来运算，其中多余的threads不应该参与运算，所以下面if需要判断。
    if (i < N){
        C[i] = A[i] + B[i];
    }

}

int main(){
    int N = 32;
    size_t size = N * sizeof(float);
    float* h_A = (float*)malloc(size);
    float* h_B = (float*)malloc(size);
    float* h_C = (float*)malloc(size);

    float* d_A;
    float* d_B;
    float* d_C;
    
    cudaMalloc(&d_A,size);
    cudaMalloc(&d_B,size);
    cudaMalloc(&d_C,size);
    
    cudaMemcpy(d_A, h_A,size,cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B,size,cudaMemcpyHostToDevice);

    int threadsPerBlock = 256;
    // 由于threadsPerBlock通常是定的，所以我们需要调整的一般是Blocks的数量
    // 考虑到如果N < threadsPerBlock,那么N/threadsPerBlock就是0，但这是是少需要一个block
    // 如果N = threadsPerBlock，这时候如果(N + threadsPerBlock)/threadsPerBlock就是2，但其实这时候只需要一个block
    // 综上，应该是下式才满足条件
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;

    VecAdd<<<blocksPerGrid,threadsPerBlock>>>(d_A,d_B,d_C,N);
    
    cudaMemcpy(h_C,d_C,cudaMemcpyDeviceToHost);
    
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
}

三、矩阵的乘法运算

该矩阵相乘的例子，不仅帮助你更好的理解kernel函数是怎样操作矩阵的，内容还包括通过block内shared memory来减少从global memory读取和写入数据的次数，从而加速矩阵相乘的速度的实现方法。

首先是没有shared memory的代码实现：

typedef struct{
    int width;
    int height;
    float* elements;
} Matrix;

#define BLOCK_SIZE 16

__global__ void MatMulKernel(const Matrix, const Matrix, Matrix);

void MatMul(const Matrix A, const Matrix B, Matrix C){
    Matrix d_A;
    d_A.width = A.width;
    d_A.height = A.height;
    size_t size = d_A.width * d_A.height * sizeof(float);
    cudaMalloc(&d_A.elements,size);
    cudaMemcpy(d_A.elements,A.elements,size,cudaMemcpyHostToDevice);

    Matrix d_B;
    d_B.width = B.width;
    d_B.height = B.height;
    size_t size = d_B.width * d_B.height * sizeof(float);
    cudaMalloc(&d_B.elements,size);
    cudaMemcpy(d_B.elements,B.elements,size,cudaMemcpyHostToDevice);

    Matrix d_C；
    d_C.width = C.width;
    d_C.height = C.height;
    size_t size = d_C.width * d_C.height * sizeof(float);
    cudaMalloc(&d_C.elements,size);
    
    dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
    // 列索引和行索引划分，结果C的形状应该是A.height*B.width；
	// 因为grid的形状是(x,y)类型，也即第一个代表的是列数，所以应该写成
	// (B.width / dimBlock.x,A.height / dimBlock.y)，表示C有x列，y行
    dim3 dimGrid(B.width / dimBlock.x,A.height / dimBlock.y);
    MatMulKernel<<<dimBlock,dimGrid>>>(d_A,d_B,d_C);

    cudaMemcpy(C.elements,d_C.elements,size,cudaMemcpyDeviceToHost);
    cudaFree(d_A.elements);
    cudaFree(d_B.elements);
    cudaFree(d_C.elements);
}

__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C){
    float Cvalue = 0;
    // 这里求对应A矩阵的第row行以及B矩阵的第col列
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
	//这里是根据矩阵相乘的运算法则，也即A的一行乘以B的一列，然后求和得到一个值。
	//例如，C[1,2]就是A的第一行和B的第二列每个元素相乘再相加，得到的一个值。
    for(int i=0;i<A.width;i++){
        Cvalue += A.elements[row * A.width + i] * B.elements[i * B.width + col];
    }
    C.elements[row * C.width + col] = Cvalue;
}

 代码很好理解，每个thread负责通过一个循环来计算输出矩阵C中的一个数值，需要注意的是，C矩阵中的每个数值计算，都需要一个thread从global memory中读取数据，计算完成，然后写回。下图表达的就是这样的一个计算实例：

              

 所以总的读取和写入次数是非常频繁的，为了减少访问global memory的次数，可以对C矩阵分块求取，C矩阵中的每个分块矩阵由一个block负责计算，这block一次将数据读入block内的共享内存，然后thread在共享内存上计算，最后写入到global memory，这大大减少了访问global memory的次数，加速了计算。
 

// 加的这个stride就是为了表示原矩阵width，
// 这里的width和height可以是子矩阵的width和height，elements是子矩阵的首地址
// M(row,col) = *(M.elements + row * M.stride + col)
// 当然要说的一句是，既然是子矩阵，那么就像子集一样，有可能就是自己本身。
// 所以，其实width和height也可以是自己实际的宽高，那么此时elements的值就是第一个元素的地址。
typedef struct{
    int width;
    int height;
    int stride;
    float* elements;
} Matrix;

#define BLOCK_SIZE 16

__device__ void SetElement(Matrix A, int row, int col, float value){
    // 虽然A是一个局部变量，但是A里面的elements指向的东西没被深拷贝
    A.elements[row * A.stride + col] = value; 
}
__device__ float GetElement(Matrix A,int row, int col){
    return A.elements[row * A.stride + col];
}

// 这里的获取子矩阵，其实就是修改矩阵的elements，使其指向一个新的元素。
// 这么说吧，假如子矩阵都是(3,3)的，那么第一个子矩阵elements的值是A[0][0]的地址，
// 那么A矩阵最左上角的3*3的部分就是第一个子矩阵。那以此类推，第二个子矩阵假如是第一个子矩阵往右平移3个元素，
// 那么此时只需要将elements的值修改为A[0][3]这个元素的地址，也即指向该元素。由于我们设置了hw，再加上独特的
// 访问方式，我们就访问不到其他的值，所以间接的就认为我们截取了原矩阵的部分元素，也即子矩阵。
// 所以，子矩阵的获取，就是使elements属性的值等于你想要截取的子矩阵最左上角那个元素的地址即可。
__device__ Matrix getSubMatrix(Matrix A,int row, int col){
    Matrix Asub;
    Asub.width = BLOCK_SIZE;
    Asub.height = BLOCK_SIZE;
    Asub.stride = A.stride;
	// 这里的row和col在下面使用的时候，其实指的是block的索引位置
    Asub.elements = A.elements + (BLOCK_SIZE * row) * stride + (BLOCK_SIZE*col) // 相当于乘以BLOCK_SIZE
    return Asub;
}


__global__ void MatMulKernel(const Matrix, const Matrix, Matrix);

void MatMul(const Matrix A, const Matrix B, Matrix C){
    Matrix d_A;
    d_A.width = A.width;
    d_A.height = A.height;
    size_t size = d_A.width * d_A.height * sizeof(float);
    cudaMalloc(&d_A.elements,size);
    cudaMemcpy(d_A.elements,A.elements,size,cudaMemcpyHostToDevice);

    Matrix d_B;
    d_B.width = B.width;
    d_B.height = B.height;
    size_t size = d_B.width * d_B.height * sizeof(float);
    cudaMalloc(&d_B.elements,size);
    cudaMemcpy(d_B.elements,B.elements,size,cudaMemcpyHostToDevice);

    Matrix d_C;
    d_C.width = C.width;
    d_C.height = C.height;
    size_t size = d_C.width * d_C.height * sizeof(float);
    cudaMalloc(&d_C.elements,size);
    
    dim3 dimBlock(BLOCK_SIZE,BLOCK_SIZE);
    dim3 dimGrid(B.width / dimBlock.x,A.height / dimBlock.y);
    MatMulKernel<<<dimBlock,dimGrid>>>(d_A,d_B,d_C);
	
	
    cudaMemcpy(C.elements,d_C.elements,size,cudaMemcpyDeviceToHost);
    cudaFree(d_A.elements);
    cudaFree(d_B.elements);
    cudaFree(d_C.elements);
}

__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C){
    
	int block_row = blockIdx.y;
    int block_col = blockIdx.x;
	
	// 根据每个block所处的位置，截取该block最左上角的元素，它将作为子矩阵的最左上角的元素。
    Matrix Csub = GetSubMatrix(C,block_row,block_col);
    float Cvalue = 0;
    int row = threadIdx.y;
    int col = threadIdx.x;

    for(int i=0;i<(A.width / BLOCK_SIZE);++i){
        
		// 这里结合下面的图比较好理解。
		// 我们之前定义了和C大小一致的grid，其中每个线程其实对应的就是C的每个元素。
		// 这里就是从C的每个元素出发，根据矩阵相乘的原理，来选择对应的A和B的元素进行计算得到结果。
		
		// 这里之所以有一个i的for循环，其实看图就很清晰了，就是为了遍历A的一整行和B的一整列，
		// 然后分别相乘求和，这样得到的数据才是最终C中一个元素的值。所以，这个i就是为了遍历多个block，
		// 因为一个block并没有包含A或B的一整行或一整列。看图就可以知道，这里就需要两个block才包含整行或者整列。
        Matrix Asub = GetSubMatrix(A,block_row,i);
        Matrix Bsub = GetSubMatrix(B,i,block_col);
        
		// 定义shared内存变量
        __shared__ float As[BLOCK_SIZE][BLOCK_SIZE];
        __shared__ float Bs[BLOCK_SIZE][BLOCK_SIZE];
        
		// 这里其实有一点非常需要注意的就是，同时会有很多线程执行这个核函数，也即很多线程会执行下面这两句赋值语句
		// 所以，表面上C中的一个位置，只对应一对(col,row)的值，下面As和Bs也只会给各自的一个位置赋值，但其实不是的。
        As[row][col] = GetElement(Asub,row,col);
        Bs[row][col] = GetElement(Bsub,row,col);
        
		// 保证整个子矩阵都被读入
        // 也就是一个block的所有thread都执行到这里
        __syncthreads();

        // 因为有循环，这里会执行多次，现在假设是两次，那么第二次和第一次不同之处在哪呢？
		// 其实就在As和Bs这两个矩阵的赋值部分，因为因为i变了，所以获取的子矩阵Asub和Bsub也变了，所以赋值也变了
		// 那下面这里的语句虽然没有改变，但是参与计算的数值已经发生变化，所以Cvalue的值遍历了A的一整行和B的一整列。
        for(int e=0; e<BLOCK_SIZE;++e){
            Cvalue += As[row][e] * Bs[e][col];
        }
        // 保证这个block算完，然后下一次训练读入下一轮的sub matrix
        __syncthreads();
    }
    SetElement(Csub,row,col,Cvalue)
}

 假设A矩阵是m行n列的，B矩阵是p行q列的，如果用前一种方法，A和B矩阵从global memory中被读取的次数是n + p次，后一种方法是(n+ p)/block_size次。