【参加CUDA线上训练营】-Day3

【参加CUDA线上训练营】-Day3

Programming optimization for data copy

Before Unified Memory Model

在传统计算机架构中，CPU与GPU是分开独立工作的设备，数据传输需要经PCI-e通道传输。而在嵌入式GPU或片上系统SOC等集成环境。CPU/GPU/DSP/Modem被高度集成化。他们共享系统总线和内存。这就为减小IO瓶颈提供了可行的思路。

在何老师介绍进行今天的课程之前，我们将任务分配给GPU执行是这样的。

1. Malloc函数分别分配Host&Device端的内存空间并初始化
2. 执行cudaMemcpy()拷贝数据到Device
3. 编写核函数并执行运算
4. 再次调用cudaMemcpy()拷贝数据到Host
   我们在编程中必须在CPU和GPU两端都进行内存分配，并不断地进行手动copy，来保证两端的内存一致。

Unified Memory

Unified memory在程序员的视角中，维护了一个统一的内存池，在CPU与GPU中共享。使用了单一指针进行托管内存，由系统来自动地进行内存迁移。

How it Works

1. 使用cudaMallocManaged或__managed__修饰符来分配内存,这里定义的数据结构不需要分开处理Host&Device
2. 编写核函数并执行运算
3. 而因此在launch kernel后需要调用cudaDeviceSynchronize进行同步。

Obvious Advantages

通过运行比较代码


我们发现最快的零拷贝代码时间最快，是同样GPU运算但拷贝的时间的一半不到，CPU最慢。因此有效提高了运算速度。

Code Example

#include <stdio.h>
#include <math.h>
#include "error.cuh"

#define BLOCK_SIZE 16
__managed__ int a[1000 * 1000];
__managed__ int b[1000 * 1000];
__managed__ int c_gpu[1000 * 1000];
__managed__ int c_cpu[1000 * 1000];

__global__ void gpu_matrix_mult(int* a, int* b, int* c, int m, int n, int k)
{
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    int sum = 0;
    if (col < k && row < m)
    {
        for (int i = 0; i < n; i++)
        {
            sum += a[row * n + i] * b[i * k + col];
        }
        c[row * k + col] = sum;
    }
}
__global__ void gpu_matrix_mult_shared(int* d_a, int* d_b, int* d_result, int M, int N, int K)
{
    __shared__ int tile_a[BLOCK_SIZE][BLOCK_SIZE];
    __shared__ int tile_b[BLOCK_SIZE][BLOCK_SIZE];

    int row = blockIdx.y * BLOCK_SIZE + threadIdx.y;
    int col = blockIdx.x * BLOCK_SIZE + threadIdx.x;
    int tmp = 0;
    int idx;

    for (int sub = 0; sub <= N/BLOCK_SIZE; ++sub)
    {
        int r = row;
        int c = sub * BLOCK_SIZE + threadIdx.x;
        idx = r * N + c;

        if (r >= M || c >= N)
        {
            tile_a[threadIdx.y][threadIdx.x] = 0;
        }
        else
        {
            tile_a[threadIdx.y][threadIdx.x] = d_a[idx];
        }

        r = sub * BLOCK_SIZE + threadIdx.y;
        c = col;
        idx = r * K + c;
        if (c >= K || r >= N)
        {
            tile_b[threadIdx.y][threadIdx.x] = 0;
        }
        else
        {
            tile_b[threadIdx.y][threadIdx.x] = d_b[idx];
        }
        __syncthreads();
        for (int k = 0; k < BLOCK_SIZE; ++k)
        {
            tmp += tile_a[threadIdx.y][k] * tile_b[k][threadIdx.x];
        }
        __syncthreads();
    }
    if (row < M && col < K)
    {
        d_result[row * K + col] = tmp;
    }
}
void cpu_matrix_mult(int* a, int* b, int* h_result, int m, int n, int k) {
    for (int i = 0; i < m; ++i)
    {
        for (int j = 0; j < k; ++j)
        {
            int tmp = 0.0;
            for (int h = 0; h < n; ++h)
            {
                tmp += a[i * n + h] * b[h * k + j];
            }
            h_result[i * k + j] = tmp;
        }
    }
}

int main(int argc, char const* argv[])
{
    int m = 1000;
    int n = 1000;
    int k = 1000;

    cudaEvent_t start, stop_cpu, stop_gpu;
    CHECK(cudaEventCreate(&start));
    CHECK(cudaEventCreate(&stop_cpu));
    CHECK(cudaEventCreate(&stop_gpu));


    for (int i = 0; i < m; ++i) {
        for (int j = 0; j < n; ++j) {
            a[i * n + j] = 0*rand() % 1024+1;
        }
    }

    for (int i = 0; i < n; ++i) {
        for (int j = 0; j < k; ++j) {
            b[i * k + j] = 0 * rand() % 1024 +1;
        }
    }

    CHECK(cudaEventRecord(start));
    cudaEventQuery(start);

    unsigned int grid_rows = (m + BLOCK_SIZE - 1) / BLOCK_SIZE;
    unsigned int grid_cols = (k + BLOCK_SIZE - 1) / BLOCK_SIZE;
    dim3 dimGrid(grid_cols, grid_rows);
    dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);


    gpu_matrix_mult_shared << <dimGrid, dimBlock >> > (a, b, c_gpu, m, n, k);

    CHECK(cudaEventRecord(stop_gpu));
    CHECK(cudaEventSynchronize(stop_gpu));

    cpu_matrix_mult(a, b, c_cpu, m, n, k);
    CHECK(cudaEventRecord(stop_cpu));
    CHECK(cudaEventSynchronize(stop_cpu));
    float elapsed_time_cpu, elapsed_time_gpu;
    CHECK(cudaEventElapsedTime(&elapsed_time_gpu, start, stop_gpu));
    CHECK(cudaEventElapsedTime(&elapsed_time_cpu, stop_gpu, stop_cpu));
    printf("GPU Time = %g ms.\n", elapsed_time_gpu);
    printf("CPU Time = %g ms.\n", elapsed_time_cpu);

    CHECK(cudaEventDestroy(start));
    CHECK(cudaEventDestroy(stop_cpu));
    CHECK(cudaEventDestroy(stop_gpu));

    

    int ok = 1;
    for (int i = 0; i < m; ++i)
    {
        for (int j = 0; j < k; ++j)
        {
            //printf("GPU: % d; CPU: %d; ", h_c[i * k + j], h_cc[i * k + j]);
            if (fabs(c_gpu[i * k + j] - c_cpu[i * k + j]) > (1.0e-10))
            {

                ok = 0;
            }
            //printf("\n");
        }
    }

    if (ok)
    {
        printf("Pass!!!\n");
    }
    else
    {
        printf("Error!!!\n");
    }

    return 0;
}

CUDA Atom Operate

What’s mean of atomicity

CUDA的原子操作可以理解为对一个变量进行“读取-修改-写入”这三个操作的一个最小单位的执行过程，这个执行过程不能够再分解为更小的部分，在它执行过程中，不允许其他并行线程对该变量进行读取和写入的操作。基于这个机制，原子操作实现了对在多个线程间共享的变量的互斥保护，确保任何一次对变量的操作的结果的正确性。

优势

保证了变量不会被两个及以上线程同时调用，确保了程序执行的顺序性

缺点

程序失去并发性流水性，降低了性能。

Where Need the Atom Operate

1. 数据写入过程中
2. 特殊中断处理操作
3. 数据读取时但有进程想写入
4. 其他有读写顺序性要求的场合

Atomicity Function

1. int atomicAdd(int* address, int val);
   读取位于全局或共享存储器中地址address 处的32 位或64 位字old，计算(old + val)，并将结果存储在存储器的同一地址中。这三项操作在一次原子事务中执行。该函数将返回old。只有全局存储器支持64 位字。

2. int atomicSub(int* address, int val);
   读取位于全局或共享存储器中地址address 处的32 位字old，计算(old - val)，并将结果存储在存储器的同一地址中。这三项操作在一次原子事务中执行。该函数将返回old。

3. int atomicExch(int* address, int val);
   读取位于全局或共享存储器中地址address 处的32 位或64 位字old，并将val 存储在存储器的同一地址中。这两项操作在一次原子事务中执行。该函数将返回old。只有全局存储器支持64 位字。

4. atomicMin()
   读取位于全局或共享存储器中地址address 处的32 位字old，计算old 和val 的最小值，并将结果存储在存储器的同一地址中。这三项操作在一次原子事务中执行。该函数将返回old。

5. int atomicMax(int* address, int val);
   读取位于全局或共享存储器中地址address 处的32 位字old，计算old 和val 的最大值，并将结果存储在存储器的同一地址中。这三项操作在一次原子事务中执行。该函数将返回old。

6. int atomicXor(int* address, int val);
   读取位于全局或共享存储器中地址address 处的32 位字old，计算 (old 异或val)，并将结果存储在存储器的同一地址中。这三项操作在一次原子事务中执行。该函数将返回old。

Array SUM example

#include <stdio.h>    
#include <stdlib.h>   
 
 
#define SIZE 1024
 
__global__ void kernel(int size, unsigned int *histo)
{
    int i = threadIdx.x + blockIdx.x * blockDim.x;
    if (i < size)
    {
        //*histo+=i;
        atomicAdd(histo, i);
    }
}
 
int main(void)
{
    int threadSum = 0;
 
    //分配内存并拷贝初始数据
    unsigned int *dev_histo;
 
    cudaMalloc((void**)&dev_histo, sizeof(int));
    cudaMemcpy(dev_histo, &threadSum, sizeof(int), cudaMemcpyHostToDevice);
 
    // kernel launch - 2x the number of mps gave best timing  
   
 
    int blocks = prop.multiProcessorCount;
    //确保线程数足够
    kernel << <blocks * 2, (SIZE + 2 * blocks - 1) / blocks / 2 >> > (SIZE, dev_histo);
 
    //数据拷贝回CPU内存
    cudaMemcpy(&threadSum, dev_histo, sizeof(int), cudaMemcpyDeviceToHost);
    printf("Threads SUM：%d\n", threadSum);
    getchar();
    cudaFree(dev_histo);
    return 0;
}

Learning Feedback

通过更加深入的学习，我对CUDA背后的工程师由衷的感叹。CUDA帮我们自动化完成了同步、异步、互斥，数据传输、内存管理等很多编程问题，使得程序员专注于功能实现本身。让我学到很多受益匪浅！！