CUDA编程基础

典型的CUDA程序的执行流程如下：

分配host内存，并进行数据初始化；
分配device内存，并从host将数据拷贝到device上；
调用CUDA的核函数在device上完成指定的运算；
将device上的运算结果拷贝到host上；
释放device和host上分配的内存。
kernel
kernel是CUDA中一个重要的概念，kernel是在device上线程中并行执行的函数，核函数用__global__符号声明，在调用时需要用<<<grid, block>>>来指定kernel要执行的线程数量，在CUDA中，每一个线程都要执行核函数，并且每个线程会分配一个唯一的线程号thread ID，这个ID值可以通过核函数的内置变量threadIdx来获得。

如上图：
1，一个kernal启动的所有线程叫做一个网格（grid），同一个网格上的线程共享的是一个相同的全局内存空间
2 一个grid又分为很多的线程块（block），一个线程块里面包含很多的线程，这是第二个层次
gird和block都是定义为dim3类型的变量。kernal通过执行配置<<<grid,block>>>来指定线程数和结构。

dim3 grid(3, 2);
dim3 block(5, 3);
kernel_fun<<< grid, block >>>(prams...);

• 1.使用N个线程块，每一个线程块只有一个线程，即

  dim3 dimGrid(N);
  dim3 dimBlock(1);

此时的线程号的计算方式就是

threadId = blockIdx.x;
其中threadId的取值范围为0到N-1。对于这种情况，我们可以将其看作是一个列向量，列向量中的每一行对应一个线程块。列向量中每一行只有1个元素，对应一个线程。

• 2.使用M×N个线程块，每个线程块1个线程

由于线程块是2维的，故可以看做是一个M*N的2维矩阵，其线程号有两个维度，即：

dim3 dimGrid(M,N);
dim3 dimBlock(1);

其中

blockIdx.x 取值0到M-1
blcokIdx.y 取值0到N-1
这种情况一般用于处理2维数据结构，比如2维图像。每一个像素用一个线程来处理，此时需要线程号来映射图像像素的对应位置，如

pos = blockIdx.y * blcokDim.x + blockIdx.x; //其中gridDim.x等于M

• 3.使用一个线程块，该线程具有N个线程，即

  dim3 dimGrid(1);
  dim3 dimBlock(N);

此时线程号的计算方式为

threadId = threadIdx.x;

其中threadId的范围是0到N-1，对于这种情况，可以看做是一个行向量，行向量中的每一个元素的每一个元素对应着一个线程。

• 4.使用M个线程块，每个线程块内含有N个线程，即

  dim3 dimGrid(M);
  dim3 dimBlock(N);

这种情况，可以把它想象成二维矩阵，矩阵的行与线程块对应，矩阵的列与线程编号对应，那线程号的计算方式为

threadId = threadIdx.x + blcokIdx*blockDim.x;

上面其实就是把二维的索引空间转换为一维索引空间的过程。

• 5.使用M×N的二维线程块，每一个线程块具有P×Q个线程，即

  dim3 dimGrid(M, N);
  dim3 dimBlock(P, Q);

这种情况其实是我们遇到的最多情况，特别适用于处理具有二维数据结构的算法，比如图像处理领域。

其索引有两个维度

threadId.x = blockIdx.x*blockDim.x+threadIdx.x;
threadId.y = blockIdx.y*blockDim.y+threadIdx.y;

上述公式就是把线程和线程块的索引映射为图像像素坐标的计算方法。

所以，一个线程需要两个内置的坐标变量（blockIdx，threadIdx）来唯一标识，它们都是dim3类型变量，其中blockIdx指明线程所在grid中的位置，而threaIdx指明线程所在block中的位置，如图中的Thread (1,1)满足：

threadIdx.x = 1
threadIdx.y = 1
blockIdx.x = 1
blockIdx.y = 1

一个线程块中的线程是放在同一个流式多处理器（sm）上的，但是单个sm的资源是有有限的，就导致线程块中的线程数是有限的，现在GPU的线程块可以支持的线程数量是1024个对于一个2-dim的block(Dx,Dy)(Dx,Dy)，线程(x,y)(x,y)的ID值为(x+y∗Dx)(x+y∗Dx)，如果是3-dim的block(Dx,Dy,Dz)(Dx,Dy,Dz)，线程(x,y,z)(x,y,z)的ID值为(x+y∗Dx+z∗Dx∗Dy)(x+y∗Dx+z∗Dx∗Dy)。另外线程还有内置变量gridDim，用于获得网格块各个维度的大小。
kernel的这种线程组织结构天然适合vector,matrix等运算，如我们将利用上图2-dim结构实现两个矩阵的加法，每个线程负责处理每个位置的两个元素相加，代码如下所示。线程块大小为(16, 16)，然后将N*N大小的矩阵均分为不同的线程块来执行加法运算。

// Kernel定义
__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; 
    if (i < N && j < N) 
        C[i][j] = A[i][j] + B[i][j]; 
}
int main() 
{ 
    ...
    // Kernel 线程配置
    dim3 threadsPerBlock(16, 16); 
    dim3 numBlocks(N / threadsPerBlock.x, N / threadsPerBlock.y);
    // kernel调用
    MatAdd<<<numBlocks, threadsPerBlock>>>(A, B, C); 
    ...
}

内存结构

每一个线程都有自己的私有本地内存，local memory，每个线程块都有共享内存，（share Memory），可以被线程块中的所有线程共享，生命周期和线程块一致。所有的线程都可以访问全局内存，（global Memory）。还可以访问一些只读内存，常量内存（constant）和纹理内存（texture ）SM的基本执行单元是包含32个线程的线程束，所以block大小一般要设置为32的倍数。
GPU的配置属性：
在进行CUDA编程前，可以先检查一下自己的GPU的硬件配置，这样才可以有的放矢，可以通过下面的程序获得GPU的配置属性：

    int dev = 0;
    cudaDeviceProp devProp;
    CHECK(cudaGetDeviceProperties(&devProp, dev));
    std::cout << "使用GPU device " << dev << ": " << devProp.name << std::endl;
    std::cout << "SM的数量：" << devProp.multiProcessorCount << std::endl;
    std::cout << "每个线程块的共享内存大小：" << devProp.sharedMemPerBlock / 1024.0 << " KB" << std::endl;
    std::cout << "每个线程块的最大线程数：" << devProp.maxThreadsPerBlock << std::endl;
    std::cout << "每个EM的最大线程数：" << devProp.maxThreadsPerMultiProcessor << std::endl;
    std::cout << "每个EM的最大线程束数：" << devProp.maxThreadsPerMultiProcessor / 32 << std::endl;

    // 输出如下
    使用GPU device 0: GeForce GT 730
    SM的数量：2
    每个线程块的共享内存大小：48 KB
    每个线程块的最大线程数：1024
    每个EM的最大线程数：2048
    每个EM的最大线程束数：64

CUDA编程中内存管理API
首先是在device上分配内存的cudaMalloc函数：

cudaError_t cudaMalloc(void** devPtr, size_t size);

负责host和device之间数据通信的cudaMemcpy函数：

cudaError_t cudaMemcpy(void* dst, const void* src, size_t count, cudaMemcpyKind kind)

其中src指向数据源，而dst是目标区域，count是复制的字节数，其中kind控制复制的方向：cudaMemcpyHostToHost, cudaMemcpyHostToDevice, cudaMemcpyDeviceToHost及cudaMemcpyDeviceToDevice，如cudaMemcpyHostToDevice将host上数据拷贝到device上。
引用了一个实现二维数组的编程
矩阵乘法：
CPU版本：

#include <iostream>
#include <stdlib.h>
#include <sys/time.h>

#define ROWS 1024
#define COLS 1024

using namespace std;

void matrix_mul_cpu(float* M, float* N, float* P, int width)
{
    for(int i=0;i<width;i++)
        for(int j=0;j<width;j++)
        {
            float sum = 0.0;
            for(int k=0;k<width;k++)
            {
                float a = M[i*width+k];
                float b = N[k*width+j];
                sum += a*b;
            }
            P[i*width+j] = sum;
        }
}

int main()
{
    struct timeval start, end;
    gettimeofday( &start, NULL );
    float *A, *B, *C;
    int total_size = ROWS*COLS*sizeof(float);
    A = (float*)malloc(total_size);
    B = (float*)malloc(total_size);
    C = (float*)malloc(total_size);

    //CPU一维数组初始化
    for(int i=0;i<ROWS*COLS;i++)
    {
        A[i] = 80.0;
        B[i] = 20.0;
    }

    matrix_mul_cpu(A, B, C, COLS);

    gettimeofday( &end, NULL );
    int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
    cout << "total time is " << timeuse/1000 << "ms" <<endl;

    return 0;
}

梳理一下CUDA求解矩阵乘法的思路：因为C=A×B，我们利用每个线程求解C矩阵每个(x, y)的元素，每个线程载入A的一行和B的一列，遍历各自行列元素，对A、B对应的元素做一次乘法和一次加法。

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <sys/time.h> 
#include <stdio.h>
#include <math.h>
#define Row  1024
#define Col 1024

 
__global__ void matrix_mul_gpu(int *M, int* N, int* P, int width)
{
    int i = threadIdx.x + blockDim.x * blockIdx.x;
    int j = threadIdx.y + blockDim.y * blockIdx.y;
                
    int sum = 0;
    for(int k=0;k<width;k++)
    {
        int a = M[j*width+k];
        int b = N[k*width+i];
        sum += a*b;
    }
    P[j*width+i] = sum;
}
 
int main()
{
    struct timeval start, end;
    gettimeofday( &start, NULL );

    int *A = (int *)malloc(sizeof(int) * Row * Col);
    int *B = (int *)malloc(sizeof(int) * Row * Col);
    int *C = (int *)malloc(sizeof(int) * Row * Col);
    //malloc device memory
    int *d_dataA, *d_dataB, *d_dataC;
    cudaMalloc((void**)&d_dataA, sizeof(int) *Row*Col);
    cudaMalloc((void**)&d_dataB, sizeof(int) *Row*Col);
    cudaMalloc((void**)&d_dataC, sizeof(int) *Row*Col);
    //set value
    for (int i = 0; i < Row*Col; i++) {
        A[i] = 90;
        B[i] = 10;
    }
                                                                
    cudaMemcpy(d_dataA, A, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);
    cudaMemcpy(d_dataB, B, sizeof(int) * Row * Col, cudaMemcpyHostToDevice);
    dim3 threadPerBlock(16, 16);
    dim3 blockNumber((Col+threadPerBlock.x-1)/ threadPerBlock.x, (Row+threadPerBlock.y-1)/ threadPerBlock.y );
    printf("Block(%d,%d)   Grid(%d,%d).\n", threadPerBlock.x, threadPerBlock.y, blockNumber.x, blockNumber.y);
    matrix_mul_gpu << <blockNumber, threadPerBlock >> > (d_dataA, d_dataB, d_dataC, Col);
    //拷贝计算数据-一级数据指针
    cudaMemcpy(C, d_dataC, sizeof(int) * Row * Col, cudaMemcpyDeviceToHost);
                                                                                             
    //释放内存
    free(A);
    free(B);
    free(C);
    cudaFree(d_dataA);
    cudaFree(d_dataB);
    cudaFree(d_dataC);

    gettimeofday( &end, NULL );
    int timeuse = 1000000 * ( end.tv_sec - start.tv_sec ) + end.tv_usec - start.tv_usec;
    printf("total time is %d ms\n", timeuse/1000);

    return 0;
}

参考文献：

1 CUDA编程入门极简教程 - 小白将 - CSDN博客
https://blog.csdn.net/xiaohu2022/article/details/79599947 2
2【CUDA并行编程之五】计算向量的欧式距离 - 懂幸福，爱生活 - CSDN博客
https://blog.csdn.net/lavorange/article/details/42125029 3
3 cuda编程入门示例1—两个向量对应元素相乘 - Gone_HuiLin的博客 - CSDN博客
https://blog.csdn.net/Gone_HuiLin/article/details/53228386 4
4 GPU/CUDA程序初体验 向量加法 - 旭东的博客 - 博客园
https://www.cnblogs.com/xudong-bupt/p/3461163.html 5
5(1条消息)CUDA编程–并行矩阵向量乘法【80+行代码】 - 肥宅Sean - CSDN博客
https://blog.csdn.net/a19990412/article/details/85222642 6 CUDA编程之快速入门
-6 Madcola - 博客园 https://www.cnblogs.com/skyfsm/p/9673960.html 7 https://blog.csdn.net/u014030117/article/details/45952971