CUDA编程学习（一）

CUDA编程学习

一、CPU和GPU的区别

1.1 应用场景不同

CPU： CPU需要很强的通用性来处理各种不同的数据类型，同时又要逻辑判断又会引入大量的分支跳转和中断的处理。这些都使得CPU的内部结构异常复杂。
GPU： GPU面对的则是类型高度统一的、相互无依赖的大规模数据和不需要被打断的纯净的计算环境。

1.2 结构不同

• Cache（缓存）, local memory（局部内存）： CPU > GPU
• Threads(线程数): GPU > CPU
• Registers: GPU > CPU 多寄存器可以支持非常多的Thread, thread需要用到register, thread数目大，register也必须得跟着很大才行。
• SIMD Unit(单指令多数据流,以同步方式，在同一时间内执行同一条指令): GPU > CPU。
  总结：
  CPU：有强大的ALU（算术运算单元），缓存大（cache），复杂的逻辑控制单元（Control）。
  GPU：很多的ALU，很少的cache（ 缓存的目的不是保存后面需要访问的数据的，而是为thread提高服务的），数据存于Dram。用于处理重复大规模运算。

二、GPU的结构

硬件部分与软件部分的对应关系：

• Thread------SP (每个线程由每个线程处理器（SP）执行)
• Thread block-------SM (线程块由多核处理器（SM）执行)
• Grid------Device (一个kernel其实由一个grid来执行，一个kernel一次只能在一个GPU上执行)

block是软件概念，一个block只会由一个sm调度，程序员在开发时，通过设定block的属性，告诉GPU硬件，我有多少个线程，线程怎么组织。而具体怎么调度由sm的warps scheduler负责，block一旦被分配好SM，该block就会一直驻留在该SM中，直到执行结束。一个SM可以同时拥有多个blocks，但需要序列执行。

2.1 硬件结构

硬件角度从小到大：
SP（streaming processor）： 最基本的处理单元，也称为CUDA core。最后具体的指令和任务都是在SP上处理的。GPU进行并行计算，也就是很多个SP同时做处理。
SM（streaming multiprocessor）： 多个SP加上其他的一些资源组成一个SM。也叫GPU大核，其他资源如：warp scheduler，register，shared memory等。SM可以看做GPU的心脏（对比CPU核心），register和shared memory是SM的稀缺资源。CUDA将这些资源分配给所有驻留在SM中的threads。因此，这些有限的资源就使每个SM中active warps有非常严格的限制，也就限制了并行能力。
Device： GPU运算端

2.2 CUDA线程模型

从小到大：
Thread： 线程，并行的基本单位。(自己的寄存器、局部内存)
Thread Block： 线程块，互相合作的线程组。（共享内存，所有线程共享）
a) 线程组中，线程允许彼此同步
b) 可以通过共享内存快速交换数据
c) 以1维、2维或3维组织
Grid： 一组线程块（全局内存、常量内存、纹理内存）
a) 以1维、2维组织
b) 共享全局内存
Kernel: 在GPU上执行的核心程序，这个kernel函数是运行在某个Grid上的。Kernel和Grid一一对应。

• 每一个block和每个thread都有自己的ID，我们通过相应的索引找到相应的线程和线程块。
  threadIdx，blockIdx
  Block ID: 1D or 2D
  Thread ID: 1D, 2D or 3D
• 速度：register（线程） > local memory（线程） >shared memory （线程块）> global memory（Grid）
• <<<grid, block>>>指定kernel运算的线程数

三、CUDA编程

术语：

• Device=GPU
• Host=CPU
• Kernel：在device上运行的函数

3.1 CUDA编程要点

3.1.1 要点1：让函数运行在不同设备

指令	执行位置	调用位置
__device__ float DeviceFunc()	device	device
__global__ void KernelFunc()	device	host
__host__ float HostFunc()	host	host

3.1.2 要点2：host和device数据传输

指令	作用
cudaMalloc()	cuda分配内存
cudaFree()	释放cuda内存
cudaMemcpy(void *dst, void *src, size_t nbytes, enum cudaMemcpyKind direction)	拷贝函数

• cudaMemcpyHostToDevice（CPU到GPU）
• cudaMemcpyDeviceToHost（GPU到CPU）
• cudaMemcpyDeviceToDevice（GPU到GPU）

3.1.3 要点3: 用代码表示线程组织模型

• 使用Grid创建1维排列或二维排列的Block，坐标指向Block。
• 使用Block创建1、2、3维排列的Thread，坐标指向Thread。

格式：

dim3 DimGrid(100, 50);  //5000个线程块，维度是100*50
dim3 DimBlock(4, 8, 8);  //每个线层块内包含256个线程，线程块内的维度是4*8*8

1、使用N个线程块，每一个线程块只有一个线程

dim3 dimGrid(N);
dim3 dimBlock(1);
//则线程号为：
threadId = blockIdx.x;

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

dim3 dimGrid(M,N);
dim3 dimBlock(1);
//其中
//blockIdx.x 取值0到M-1
//blcokIdx.y 取值0到N-1
//
pos = blockIdx.y * blcokDim.x + blockIdx.x; //其中blcokDim.x等于M

3、使用一个线程块，该线程具有N个线程

dim3 dimGrid(1);
dim3 dimBlock(N);
//则线程号为：
threadId = threadIdx.x;

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;

3.1.3 一些简单的例子

例子：查看设备信息

#include "device_launch_parameters.h"
#include <iostream>

int main()
{
    int deviceCount;
    cudaGetDeviceCount(&deviceCount);
    for(int i=0;i<deviceCount;i++)
    {
        cudaDeviceProp devProp;
        cudaGetDeviceProperties(&devProp, i);
        std::cout << "使用GPU device " << i << ": " << devProp.name << std::endl;
        std::cout << "设备全局内存总量： " << devProp.totalGlobalMem / 1024 / 1024 << "MB" << 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 << "设备上一个线程块（Block）种可用的32位寄存器数量： " << devProp.regsPerBlock << std::endl;
        std::cout << "每个EM的最大线程数：" << devProp.maxThreadsPerMultiProcessor << std::endl;
        std::cout << "每个EM的最大线程束数：" << devProp.maxThreadsPerMultiProcessor / 32 << std::endl;
        std::cout << "设备上多处理器的数量： " << devProp.multiProcessorCount << std::endl;
        std::cout << "======================================================" << std::endl;     
        
    }
    return 0;
}

nvcc test1.cu -o test1

使用GPU device 0: TITAN X (Pascal)
设备全局内存总量： 12189MB
SM的数量：28
每个线程块的共享内存大小：48 KB
每个线程块的最大线程数：1024
设备上一个线程块（Block）种可用的32位寄存器数量： 65536
每个EM的最大线程数：2048
每个EM的最大线程束数：64
设备上多处理器的数量： 28

例子：将两个元素数目为1024×1024的float数组相加。
CPU版本：逐个元素相加

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

using namespace std;

int main()
{
    struct timeval start, end;
    gettimeofday( &start, NULL );
    float*A, *B, *C;
    int n = 1024 * 1024;
    int size = n * sizeof(float);
    A = (float*)malloc(size);
    B = (float*)malloc(size);
    C = (float*)malloc(size);

    for(int i=0;i<n;i++)
    {
        A[i] = 90.0;
        B[i] = 10.0;
    }
    for(int i=0;i<n;i++)
    {
        C[i] = A[i] + B[i];
    }
    float max_error = 0.0;
    for(int i=0;i<n;i++)
    {
        max_error += fabs(100.0-C[i]);
    }
    cout << "max_error is " << max_error << endl;
    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;
}

max_error is 0
total time is 22ms

GPU版本多线程：

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

using namespace std;

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

int main()
{
    struct timeval start, end;
    gettimeofday( &start, NULL );
    float*A, *Ad, *B, *Bd, *C, *Cd;
    int n = 1024 * 1024;
    int size = n * sizeof(float);

    // CPU端分配内存
    A = (float*)malloc(size);
    B = (float*)malloc(size);
    C = (float*)malloc(size);

    // 初始化数组
    for(int i=0;i<n;i++)
    {
        A[i] = 90.0;
        B[i] = 10.0;
    }

    // GPU端分配内存
    cudaMalloc((void**)&Ad, size);
    cudaMalloc((void**)&Bd, size);
    cudaMalloc((void**)&Cd, size);

    // CPU的数据拷贝到GPU端
    cudaMemcpy(Ad, A, size, cudaMemcpyHostToDevice);
    cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);
    cudaMemcpy(Bd, B, size, cudaMemcpyHostToDevice);

    // 定义kernel执行配置，（1024*1024/512）个block，每个block里面有512个线程
    dim3 dimBlock(512);
    dim3 dimGrid(n/512);

    // 执行kernel
    Plus<<<dimGrid, dimBlock>>>(Ad, Bd, Cd, n);

    // 将在GPU端计算好的结果拷贝回CPU端
    cudaMemcpy(C, Cd, size, cudaMemcpyDeviceToHost);

    // 校验误差
    float max_error = 0.0;
    for(int i=0;i<n;i++)
    {
        max_error += fabs(100.0 - C[i]);
    }

    cout << "max error is " << max_error << endl;

    // 释放CPU端、GPU端的内存
    free(A);
    free(B);
    free(C);
    cudaFree(Ad);
    cudaFree(Bd);
    cudaFree(Cd);
    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;
}

max error is 0
total time is 1278ms

因为cpu和gpu之间的拷贝很多。

例子：两矩阵相乘
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;
}

total time is 7617ms

GPU版本：

#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;
}

Block(16,16)   Grid(64,64).
total time is 506 ms

GPU速度：0.5s，CPU：7.6s。
参考：https://www.cnblogs.com/skyfsm/p/9673960.html