返回
Featured image of post CUDA学习

CUDA学习

CUDA学习笔记

CUDA C编程

基本概念

  • CPU:Host(主机)
  • GPU:device(设备)

一个典型的CUDA编程结构包含以下五个步骤:

  1. 分配GPU内存。
  2. 从CPU内存中拷贝数据到GPU内存。
  3. 调用CUDA内核函数来完成程序指定的运算。
  4. 将数据从GPU拷回CPU内存。
  5. 释放GPU内存空间。

调用GPU输出hello world

#include <stdio.h>

__global__ void helloFromGPU(void){
    printf("hello world form GPU!\n");
}

int main(){
    printf("hello world from CPU!\n");
    helloFromGPU<<<1,10>>>();
    cudaDeviceReset();
    return 0;
}

CUDA编程

内存管理

函数

CUDA内存管理函数如下:

标准C函数 CUDA C函数
malloc cudaMalloc
memcpy cudaMemcpy
memset cudaMemset
free cudaFree
具体解释如下:
cudaMalloc
cudaError_t cudaMalloc(void** devPtr, size_t size);

该函数负责向设备分配一定字节的线性内存,并以devPtr的形式返回指向所分配内存的指针。cudaMalloc与标准C语言中的malloc函数几乎一样,只是此函数在GPU的内存里分 配内存。通过充分保持与标准C语言运行库中的接口一致性,可以实现CUDA应用程序的轻松接入。

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

此函数从src指向的源存储区复制一定数量的字节到dst指向的目标存储区。复制方向由kind指定,其中的kind有以下几种。

  • cudaMemcpyHostToHost
  • cudaMemcpyHostToDevice
  • cudaMemcpyDeviceToHost
  • cudaMemcpyDeviceToDevice

这个函数以同步方式执行,因为在cudaMemcpy函数返回以及传输操作完成之前主应用程序是阻塞的。除了内核启动之外的CUDA调用都会返回一个错误的枚举类型cudaError_t。如果GPU内存分配成功,函数返回:

cudaSuccess
否则返回:
cudaErrorMemoryAllocation

可以使用以下CUDA运行时函数将错误代码转化为可读的错误消息:
char* cudaGetErrorString(cudaError_t error);

内存层次结构

CUDA编程模型从GPU架构中抽象出一个内存层次结构。GPU内存结构,它主要包含两部分:全局内存和共享内存. 在GPU内存层次结构中,最主要的两种内存是全局内存和共享内存。全局类似于CPU的系统内存,而共享内存类似于CPU的缓存。然而GPU的共享内存可以由CUDA C的内核直接控制。

example

下面,我们将通过一个简单的两个数组相加的例子来学习如何在主机和设备之间进行数据传输,以及如何使用CUDA C编程。数组a的第一个元素与数组b的第一个元素相加,得到的结果作为数组c的第一个元素,重复这个过程直到数组中的所有元素都进行一次运算。
首先我们用CPU实现如下:

#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <stdio.h>

void sumArraysOnHost(float* A, float* B, float* C, const int N){
    for(int idx=0; idx<N; ++idx){
        C[idx] = A[idx] + B[idx];
    }
}

void initData(float* ip, int size){
    time_t t;
    srand((unsigned int) time(&t));

    for(int i=0; i<size; ++i){
        ip[i] = (float)(rand() & 0xFF) / 10.0f;
    }
}

void printArray(float* ip, int size){
    for(int i=0; i<size; ++i){
        printf("%f ", ip[i]);
        if(i%10 + 1== 0) printf("\n");
    }
}
int main(int argc, char** argv){
    int nElem = 1024;
    size_t nBytes = nElem * sizeof(float);

    float* h_A, *h_B, *h_C;

    h_A = (float*)malloc(nBytes);
    h_B = (float*)malloc(nBytes);
    h_C = (float*)malloc(nBytes);
    
    initData(h_A, nElem);
    initData(h_B, nElem);
    initData(h_C, nElem);

    sumArraysOnHost(h_A, h_B, h_C, nElem);

    printf("\nh_A\n");
    printArray(h_A, nElem);
    printf("\nh_B\n");
    printArray(h_B, nElem);
    printf("\nh_C\n");
    printArray(h_C, nElem);


    free(h_A);
    free(h_B);
    free(h_C);
    return 0;

}

下面代码展示了CPU与GPU之间的数据通信,让运算在GPU上执行:

#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <stdio.h>

__global__ void sumArraysOnHost(float* A, float* B, float* C, const int N){
    for(int idx=0; idx<N; ++idx){
        C[idx] = A[idx] + B[idx];
    }
}

void initData(float* ip, int size){
    time_t t;
    srand((unsigned int) time(&t));

    for(int i=0; i<size; ++i){
        ip[i] = (float)(rand() & 0xFF) / 10.0f;
    }
}

void printArray(float* ip, int size){
    for(int i=0; i<size; ++i){
        printf("%f ", ip[i]);
        if(i%10 + 1== 0) printf("\n");
    }
}
int main(int argc, char** argv){
    int nElem = 1024;
    size_t nBytes = nElem * sizeof(float);


    // CPU上分配内存
    float* h_A, *h_B, *h_C;

    h_A = (float*)malloc(nBytes);
    h_B = (float*)malloc(nBytes);
    h_C = (float*)malloc(nBytes);

    // GPU上分配内存
    float* d_A, *d_B, *d_C;

    cudaMalloc((float**)&d_A,nBytes);
    cudaMalloc((float**)&d_B,nBytes);
    cudaMalloc((float**)&d_C,nBytes);

    // 给CPU数据赋值
    initData(h_A, nElem);
    initData(h_B, nElem);
    initData(h_C, nElem);

    // 将CPU数据拷贝给GPU
    cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, nBytes, cudaMemcpyHostToDevice);


    sumArraysOnHost<<<1,1>>>(d_A, d_B, d_C, nElem);

    cudaMemcpy(h_C, d_C, nBytes, cudaMemcpyDeviceToHost);


    printf("\nh_A\n");
    printArray(h_A, nElem);
    printf("\nh_B\n");
    printArray(h_B, nElem);
    printf("\nh_C\n");
    printArray(h_C, nElem);


    free(h_A);
    free(h_B);
    free(h_C);
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    return 0;

}

线程管理

当核函数在主机端启动时,它的执行会移动到设备上,此时设备中会产生大量的线程并且每个线程都执行由核函数指定的语句。

CUDA明确了线程层次抽象的概念以便于你组织线程。这是一个两层的线程层次结构,由线程块和线程块网格构成,如下图.

由一个内核启动所产生的所有线程统称为一个 网格(grid) 。同一网格中的所有线程共享相同的全局内存空间。一个网格由多个 线程块 (block)构成,一个线程块包含一组线程,同一线程块内的线程协作可以通过以下方式来实现:

  • 同步
  • 共享内存

不同块内的线程不能协作。
线程依赖以下两个变量来区分彼此:

  1. blockIdx(线程块在线程格内的索引)
  2. threadIdx(块内的线程索引)

该坐标变量是基于uint3定义的CUDA内置的向量类型,是一个包含3个无符号整数的结构,可以通过x、y、z三个字段来指定。

blockIdx.x
blockIdx.y
blockIdx.z
threadIdx.x
threadIdx.y
threadIdx.z

CUDA可以组织三维的网格和块。网格和块的维度由下列两个内置变量指定。

  1. blockDim(线程块的维度,用每个线程块中的线程数来表示)
  2. gridDim(线程格的维度,用每个线程格中的线程数来表示)
    它们是dim3类型的变量,是基于uint3定义的整数型向量,用来表示维度。当定义一个dim3类型的变量时,所有未指定的元素都被初始化为1。dim3类型变量中的每个组件可以通过它的x、y、z字段获得。

在CUDA程序中有两组不同的网格和块变量:手动定义的dim3数据类型和预定义的uint3数据类型。在主机端,作为内核调用的一部分,你可以使用dim3数据类型定义一个网格和块的维度。当执行核函数时,CUDA运行时会生成相应的内置预初始化的网格、块和线程变量,它们在核函数内均可被访问到且为unit3类型。手动定义的dim3类型的网格和块变量仅在主机端可见,而unit3类型的内置预初始化的网格和块变量仅在设备端可见。

对于一个给定的数据大小,确定网格和块尺寸的一般步骤为:

  • 确定块的大小
  • 在已知数据大小和块大小的基础上计算网格维度

要确定块尺寸,通常需要考虑:

  • 内核的性能特性
  • GPU资源的限制

启动一个CUDA核函数

<<< >>>运算符内是和核函数的执行配置

kernel_name <<<grid, block>>>(argument list);

。执行配置的第一个值是网格维度,也就是启动块的数目。第二 个值是块维度,也就是每个块中线程的数目。通过指定网格和块的维度,你可以进行以下配置:

  • 内核中线程数目
  • 内核中使用的线程布局

同一个块中的线程之间可以相互协作,不同块内的线程不能协作。 。对于一个给定的问题,可以使用不同的网格和块布局来组织你的线程。例如,假设你有32个数据元素用于计算,每8个元素一个块,需要启动4个块:

kernel_name<<<4,8>>>(argument list);

其线程布局如下图所示:

由于数据在全局内存中是线性存储的,因此可以用变量blockIdx.x和threadId.x来进行以下操作。

  • 在网格中标识一个唯一的线程
  • 建立线程和数据元素之间的映射关系

如果把所有32个元素放到一个块里,那么只会得到一个块:

kernel_name<<<1,32>>>(argument list);

如果每个块只含有一个元素,那么会有32个块:

kernel_name<<<32,1>>>(argument list);

核函数的调用与主机线程是 异步 的。核函数调用结束后,控制权立刻返回给主机端。你可以调用以下函数来强制主机端程序等待所有的核函数执行结束:

cudaError_t cudaDeviceSynchronize(void);

一些CUDA运行时API在主机和设备之间是 隐式同步 的。当使用cudaMemcpy函数在主机和设备之间拷贝数据时,主机端隐式同步,即主机端程序必须等待数据拷贝完成后才能 继续执行程序。

编写核函数

核函数是在设备端执行的代码。在核函数中,需要为一个线程 规定要进行的计算以及要进行的数据访问 。当核函数被调用时,许多不同的CUDA线程并行执行同一个计算任务。以下是用 __global__声明定义核函数:

__global__ void kernel_name(argument list);

核函数必须有一个void返回类型。

下总结了了CUDA C程序中的函数类型限定符。函数类型限定符指定一个函数在主机上执行还是在设备上执行,以及可被主机调用还是被设备调用。

限定符 执行 调用 备注
global 设备端 可以主机端调用,也可以从计算能力为3的设备中调用 返回类型必须为void
device 设备端 设备端
host 主机端 主机端 可以省略

__device__和__host__限定符可以一齐使用,这样函数可以同时在主机和设备端进行编译。

以下限制适用于所有核函数:

  • 只能访问设备内存
  • 必须具有void返回类型
  • 不支持可变数量的参数
  • 不支持静态变量
  • 显示异步行为

验证核函数

验证核函数有三种方法:

  • 编写主机函数验证核函数
  • 在Fermi及更高版本的设备端的核函数中使用printf函数。
  • ,可以将执行参数设置为«<1,1»>,因此强制用一个块和一个线程执行核函 数,这模拟了串行执行程序。这对于调试和验证结果是否正确是非常有用的,而且,如果 你遇到了运算次序的问题,这有助于你对比验证数值结果是否是按位精确的。

处理错误

可以定义如下宏进行错误检查

#define CHECK(call){                                                        \
    const cudaError_t error = call;                                         \
    if(error!= cudaSuccess){                                                \
        printf("Error: %s:%d ", __FILE__, __LINE__);                        \
        printf("code: %d, reason: %s\n", error, cudaGetErrorString(error)); \
        exit(1);                                                            \
    }                                                                       \
}                                                                           \

例如,你可以在以下代码中使用宏:

CHECK(cudaMemcpy(d_C,gpuRef, nBytes, cudaMemcpyHostToDevice));

如果内存拷贝或之前的异步操作产生了错误,这个宏会报告错误代码,并输出一个可读信息,然后停止程序。也可以用下述方法,在核函数调用后检查核函数错误:

kernel_function<<<grid, block>>>(argument list);
CHECK(cudaDeviceSynchronize())

CHECK(cudaDeviceSynchronize()) 会阻塞主机端线程的运行直到设备端所有的请求任务都结束,并确保最后的核函数启动部分不会出错。 以上仅是以调试为目的的 ,因为在核函数启动后添加这个检查点会阻塞主机端线程,使该检查点成为全局屏障。

编译并执行

综上,对example的程序进行修改后,利用GPU并行计算数组和,最终程序如下。

#include <stdlib.h>
#include <string.h>
#include <time.h>
#include <stdio.h>


#define CHECK(call){                                                        \
    const cudaError_t error = call;                                         \
    if(error!= cudaSuccess){                                                \
        printf("Error: %s:%d ", __FILE__, __LINE__);                        \
        printf("code: %d, reason: %s\n", error, cudaGetErrorString(error)); \
        exit(1);                                                            \
    }                                                                       \
}                                                                           \

void sumArraysOnHost(float* A, float* B, float* C, const int N){
    for(int idx=0; idx<N; ++idx){
        C[idx] = A[idx] + B[idx];
    }
}

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

void checkResult(float* hostRef, float* gpuRef, const int N){
    double epsilon = 1.0E-8;
    bool match = 1;
    for(int i =0; i<N; ++i){
        if(abs(hostRef[i] - gpuRef[i]) > epsilon){
            match = 0;
            printf("Arrarys don't match!\n");
            printf("host %5.2f gpu %5.2f at current %d\n",hostRef[i], gpuRef[i], i);
            break;
        }
    }
    if(match) printf("Arrays match. \n\n");
}


void initData(float* ip, int size){
    time_t t;
    srand((unsigned int) time(&t));

    for(int i=0; i<size; ++i){
        ip[i] = (float)(rand() & 0xFF) / 10.0f;
    }
}


int main(int argc, char** argv){
    printf("%s Starting.....\n", argv[0]);

    // 设置device
    int dev = 0;
    cudaSetDevice(dev);

    // 设置数据大小
    int nElem = 32;
    printf("Vector size %d\n", nElem);
    size_t nBytes = nElem * sizeof(float);


    // CPU上分配内存
    float* h_A, *h_B, *hostRef, *gpuRef;

    h_A = (float*)malloc(nBytes);
    h_B = (float*)malloc(nBytes);
    hostRef = (float*)malloc(nBytes);
    gpuRef = (float*)malloc(nBytes);

    // 给CPU数据赋值
    initData(h_A, nElem);
    initData(h_B, nElem);

    // GPU上分配内存
    float* d_A, *d_B, *d_C;

    cudaMalloc((float**)&d_A,nBytes);
    cudaMalloc((float**)&d_B,nBytes);
    cudaMalloc((float**)&d_C,nBytes);


    // 将CPU数据拷贝给GPU
    cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, nBytes, cudaMemcpyHostToDevice);

    // invoke kernel at host side
    dim3 block(nElem);
    dim3 grid(nElem/block.x);


    sumArraysOnGPU<<<grid,block>>>(d_A, d_B, d_C);

    cudaMemcpy(gpuRef, d_C, nBytes, cudaMemcpyDeviceToHost);
    
    sumArraysOnHost(h_A, h_B, hostRef, nElem);

    checkResult(hostRef, gpuRef, nElem);



    free(h_A);
    free(h_B);
    free(hostRef);
    free(gpuRef);
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    return 0;
}

执行以下命名即可对程序进行编译:

nvcc sum.cu -o sum

程序计时

利用CPU计时

可以通过time.h库中clock()函数,利用CPU对程序进行计时。如下面程序所示:

double iStart, iElaps;
iStart = clock();
sumArraysOnGPU<<<grid,block>>>(d_A, d_B, d_C,nElem);
cudaDeviceSynchronize();
iElaps = clock() - iStart;
printf("sumArraysOnGPU<<<%d,%d>>> Time elapsed %f sec\n", grid.x, block.x, iElaps);

利用nvprof计时

nvidia提供了nvprof命令可以让我们方便的查看核函数耗时,指令格式如下:

nvprof ./sumgpu.exe

设备管理

NVIDIA提供了几个查询和管理GPU设备的方法。学会如何查询GPU设备信息是很重要的,因为在运行时你可以使用它来帮助设置内核执行配置。

CUDA运行时API

在CUDA运行时API中有很多函数可以帮助管理这些设备。可以使用以下函数查询关于GPU设备的所有信息:

cudaError_t cudaGetDeviceProperties(cudaDeviceProp* prop, int device);

cudaDeviceProp结构体返回GPU设备的属性,具体属性参见cudaDeviceProp

确定最优GPU

一些系统支持多GPU。在每个GPU都不同的情况下,选择性能最好的GPU运行核函数是非常重要的。通过比较GPU包含的多处理器的数量选出计算能力最佳的GPU。如果你有一个多GPU系统,可以使用以下代码来选择计算能力最优的设备:

int main(){
    int numDevices = 0;
    cudaGetDeviceCount(&numDevices);
    if(numDevices > 1){
        int maxMultiprocessors = 0, maxDevice = 0;
        for(int device=0; device<numDevices; ++device){
            cudaDeviceProp props;
            cudaGetDeviceProperties(&props, device);
            if(maxMultiprocessors < props.multiProcessorCount){
                maxMultiprocessors = props.multiProcessorCount;
                maxDevice = device;
            }
        }
        cudaSetDevice(maxDevice);
    }
}

NVIDIA系统管理界面(nvidia-smi)命令行实用程序

shell里面输入以下命令:

nvidia-smi

CUDA执行模型

理解线程束本质

线程束:把32个线程划分到一个执行单元中

线程束和线程块

线程束是SM中基本的执行单元。当一个线程块的网格被启动后,网格中的线程块分 布在SM中。一旦线程块被调度到一个SM上,线程块中的线程会被进一步划分为线程束。 一个线程束由32个连续的线程组成,在一个线程束中,所有的线程按照单指令多线程 (SIMT)方式执行;也就是说,所有线程都执行相同的指令,每个线程在私有数据上进 行操作。下图展示了线程块的逻辑视图和硬件视图之间的关系。

然而,从硬件的角度来看,所有的线程都被组织成了一维的线程块可以被配置为一维、二维或三维的 。在一个块中,每个线程都有一个唯一的ID。用x维度作为最内层的维度,y维度作为第二个维度,z作为最外层的维度,则二维或三维线程块的逻辑布局可以转化为一维物理布局。ID如下:

  • 对于一维的线程块,唯一的线程ID被存储在CUDA的内置变量threadIdx.x中,并且,threadIdx.x中拥有连续值的线程被分组到线程束中
  • 对于一个给定的二维线程块,在一个块中每个线程的独特标识符都可以用内置变量threadIdx和blockDim来计算: $$threadIdx.y * blockDim.x + threadIdx.x $$
  • 对于一个三维线程块,计算如下: $$ threadIdx.z * blockDim.y * blockDim.x + threadIdx.y * blockDim.x + threadIdx.x $$

注意:以上均为块内ID

一个线程块的线程束的数量可以根据下式确定: $$ 一个线程块内线程束的数量 = \lceil \frac{一个线程块中线程的数量}{线程束大小} \rceil $$

如果线程块的大小不是线程束大小的偶数倍,那么在最后的线程束里有些线程就不会活跃。

线程块:逻辑角度与硬件角度

从逻辑角度来看,线程块是线程的集合,它们可以被组织为一维、二维或三维布局。
从硬件角度来看,线程块是一维线程束的集合。在线程块中线程被组织成一维布局,每32个连续线程组成一个线程束。

线程束分化

所谓线程束分化也就是指同一个线程束内不同线程在一次执行中分化出不同命令,这与SIMT是相悖的,将导致GPU性能显著下降。具体解释如下:

控制流是高级编程语言的基本构造中的一种。GPU支持传统的、C风格的、显式的控制流结构,例如,if…then…else、for和while。

CPU拥有复杂的硬件以执行分支预测,也就是在每个条件检查中预测应用程序的控制流会使用哪个分支。如果预测正确,CPU中的分支只需付出很小的性能代价。如果预测不正确,CPU可能会停止运行很多个周期,因为指令流水线被清空了。我们不必完全理解为什么CPU擅长处理复杂的控制流。这个解释只是作为对比的背景。

GPU是相对简单的设备,它没有复杂的分支预测机制。一个线程束中的所有线程在同一周期中必须执行相同的指令,如果一个线程执行一条指令,那么线程束中的所有线程都必须执行该指令。如果在同一线程束中的线程使用不同的路径通过同一个应用程序,这可能会产生问题。例如,思考下面的语句:

if(cond) {
    ...
} else {
    ...
}

设在一个线程束中有16个线程执行这段代码,cond为true,但对于其他16个来说cond为false。一半的线程束需要执行if语句块中的指令,而另一半需要执行else语句块中的指令。在同一线程束中的线程执行不同的指令,被称为线程束分化。我们已经知道,在一个线程束中所有线程在每个周期中必须执行相同的指令,所以线程束分化似乎会产生一个悖论。

如果一个线程束中的线程产生分化,线程束将连续执行每一个分支路径,而禁用不执行这一路径的线程 。线程束分化会导致性能明显地下降。在前面的例子中可以看到,线程束中并行线程的数量减少了一半:只有16个线程同时活跃地执行,而其他16个被禁用了。条件分支越多,并行性削弱越严重。

注意,线程束分化只发生在同一个线程束中。在不同的线程束中,不同的条件值不会引起线程束分化。

example

例如,假设有两个分支,下面展示了简单的算术内核示例。我们可以用一个偶数和奇数线程方法来模拟一个简单的数据分区,目的是导致线程束分化。该条件 $(tid%2==0)$使偶数编号的线程执行if子句,奇数编号的线程执行else子句。

其中mathKernel1以线程为粒度,mathKernel2以线程束为粒度。

一维网格和块

#include <stdio.h>

// 以线程为粒度
__global__ void mathKernel1(float *c) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float ia, ib;
    ia = ib = 0.0f;

    if (tid % 2 == 0){
        ia = 100.0f;
    } else {
        ib = 200.0f;
    }
    c[tid] = ia + ib;
}

// 以线程束为粒度
__global__ void mathKernel2(float *c) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    float ia, ib;
    ia = ib = 0.0f;

    if ((tid / warpSize) % 2 == 0) {
        ia = 100.0f;
    } else {
        ib = 200.0f;
    }
    c[tid] = ia + ib;
}

int main(int argc, char **argv) {
    // set up device
    int dev = 0;
    cudaDeviceProp deviceProp;
    cudaGetDeviceProperties(&deviceProp, dev);
    printf("%s using Device %d: %s\n", argv[0], dev, deviceProp.name);

    // set up data size
    int size = 64;
    int blocksize = 64;

    if(argc > 1) blocksize = atoi(argv[1]);

    if(argc > 2) size      = atoi(argv[2]);

    printf("Data size %d ", size);

    // set up execution configuration
    dim3 block (blocksize, 1);
    dim3 grid  ((size + block.x - 1) / block.x, 1);
    printf("Execution Configure (block %d grid %d)\n", block.x, grid.x);

    // allocate gpu memory
    float *d_C;
    size_t nBytes = size * sizeof(float);
    cudaMalloc((float**)&d_C, nBytes);


    // call kernel 1
    mathKernel1<<<grid, block>>>(d_C);

    // call kernel 3
    mathKernel2<<<grid, block>>>(d_C);

    // free gpu memory and reset divece
    cudaFree(d_C);
    cudaDeviceReset();
    return EXIT_SUCCESS;
}

利用nvprof命令得到两种情况下GPU执行时间如下,可以看出,以线程束为可以明显提高GPU性能。

分支效率

对分支效率定义如下: $$ 分支效率 = 100 * \frac{分支数-分化分支数}{分支数} $$

运行vprof --metrics branch_efficiency ./branch命令可查看分支效率,如下:

奇怪的是,没有报告显示出有分支分化(即分支效率是100%)。这个 奇怪的现象是CUDA编译器优化导致的结果 ,它将短的、有条件的代码段的断定指令取代了分支指令(导致分化的实际控制流指令)。

在分支预测中,根据条件,把每个线程中的一个断定变量设置为1或0。这两种条件流路径被完全执行,但只有断定为1的指令被执行。断定为0的指令不被执行,但相应的线程也不会停止。这和实际的分支指令之间的区别是微妙的,但理解它很重要。只有在条件语句的指令数小于某个阈值时,编译器才用断定指令替换分支指令。因此,一段很长的代码路径肯定会导致线程束分化。

使用nvcc -G branch.cu -o branch命令重新编译,可以让编译器关闭优化

–device-debug (-G)
Generate debug information for device code. Turns off all optimizations.
Don’t use for profiling; use -lineinfo instead.

我们可以看出mathKernel1的分支效率已经下降到了 $83.3 %$ !关闭优化后的分支效率

另外可通过nvprof --events branch,divergent_branch ./branch命令获得分支和分化分支的事件计数器。

资源分配

线程束的本地执行上下文主要由以下资源组成:

  • 程序计数器
  • 寄存器
  • 共享内存

线程束上下文切换没有损失。

每个SM都有32位的寄存器组,它存储在寄存器文件中,并且可以在线程中进行分配,同时固定数量的共享内存用来在线程块中进行分配。对于一个给定的内核,同时存在于同一个SM中的线程块和线程束的数量取决于在SM中可用的且内核所需的寄存器和共享内存的数量

有下列两种情况:

  • 若每个线程消耗的寄存器越多,则可以放在一个SM中的线程束就越少。如果可以减少内核消耗寄存器的数量,那么就可以同时处理更多的线程束。
  • 若一个线程块消耗的共享内存越多,则在一个SM中可以被同时处理的线程块就会变少。如果每个线程块使用的共享内存数量变少,那么可以同时处理更多的线程块。

资源可用性通常会限制SM中常驻线程块的数量。每个SM中寄存器和共享内存的数量因设备拥有不同的计算能力而不同。如果每个SM没有足够的寄存器或共享内存去处理至少一个块,那么内核将无法启动。一些关键的限度如下表所示。

计算资源(如寄存器和共享内存)已分配给线程块时,线程块被称为活跃的块。它所包含的线程束被称为活跃的线程束。活跃的线程束可以进一步被分为以下3种类型(阻塞,执行,就绪):

  • 选定的线程束
  • 阻塞的线程束
  • 符合条件的线程束

如果同时满足以下两个条件则线程束符合执行条件。

  • 32个CUDA核心可用于执行
  • 当前指令中所有的参数都已就绪

由于 计算资源是在线程束之间进行分配 的,而且在线程束的整个生存期中都保持在芯片内,因此线程束上下文的切换是非常快的

延迟隐蔽

利用率与常驻线程束的数量直接相关。在指令发出和完成之间的时钟周期被定义为 指令延迟。当每个时钟周期中所有的线程调度器都有一个符合条件的线程束时,可以达到计算资源的完全利用。这就可以保证,通过在其他常驻线程束中发布其他指令,可以隐藏每个指令的延迟。

考虑到指令延迟,指令可以被分为两种基本类型:

  • 算术指令(计算)
  • 内存指令(I/O)

算术指令延迟是一个算术操作 从开始到它产生输出之间的时间 。内存指令延迟是指 发送出的加载或存储操作和数据到达目的地之间的时间。对于每种情况,相应的延迟大约为:

  • 算术操作为10~20个周期
  • 全局内存访问为400~800个周期

你可能想知道如何估算隐藏延迟所需要的活跃线程束的数量。利特尔法则(Little’s Law)可以提供一个合理的近似值。它起源于队列理论中的一个定理,它也可以应用于GPU中: $$ 所需线程束数量=延迟×吞吐量 $$

假设在内核里一条指令的平均延迟是5个周期。为了保持在每个周期内执行6个线程束的吞吐量,则至少需要30个未完成的线程束。

内存模型

寄存器

使用以下命令查看当前.cu文件使用了多少寄存器,全局内存,常量内存和共享内存

nvcc -Xptxas -v test.cu

local memory

可通过__lsLocal()函数判断变量是否在local memory上

shared memory

通过__shared__修饰指定声明共享的变量

动态指定shared memeory大小

见下面例子

static __global__ void kernel(){
    // 使用extern声明外部的动态大小共享内存,由启动核函数的第三个参数指定
    extern __shared__ int shared_array[];

    if(threadIdx.x == 0){
        shared_array[0]=blockIdx.x;
    }

    // 线程同步
    __syncthreads();

    printf("%d,%d, shared_array[0]= %d\n", blockIdx.x, threadIdx.x, shared_array[0]);
}


int main(){
    kernel<<<2,2,sizeof(int)*5>>>();
    coudaDeviceSynchronize();
}

global memory

Pinned Memory(页锁定内存)

  • Pageable Memory
    通过mallocnew得到的内存
  • Pinned Memory 通过cudaMallocHostcudaHostAlloc可以显式的分配Pinned Memory. PM会常驻在物理内存上, 不会被交换。 (用cudaFreeHost释放内存);

Zero Copy Memory(零拷贝内存)

Zero Copy Memory:即内存的复制过程,没有CPU参与,直接由GPU和内存条操作。(DMA)

ZCM实质上就是PM映射到device的空间地址,本质上的等价的.他由cudaMallocHostcudaHostAlloc实现分配, cudaFreeHost释放。cudaHostGetDevicePointer可以获取分配的Pinned Memory映射到device的地址.

Unified Memory(统一内存)

由统一内存管理的内存。Unified Virtual Addressing(统一虚拟地址)将CPU和GPU内存看作一个整体进行管理和使用。分配的 内存可以CPU/GPU直接访问.

cudaMallocManaged分配,cudaFree释放

例子

static __global__ void test_kernel(float* array){
    array[threadIdx.x] = threadIdx.x;
}

int main(){
    int num=5;
    float* array = nullptr;
    cudaMallocManaged(&array, sizeof(float)*num);
    
    test_kernel<<<1, num>>(array);
    cudaDeviceSynchronoze();

    for(int i=0; i<num;++i){
        printf("array[%d]= %f\n",i,array[i]);
    }

    cudaFree(array);

}

constant memory

__constant__修饰,例子如下:

__constant__ float warp_matrix[6] = {1,2,3,4,5,6};

static __global__ void test_kernel(){
    printf("wrap_matrix[%d]=%f\n",threadIdx.x, warp_matrix[threadIdx.x]);
}

int main(){
    float host_warp_matrix[6]={6,5,3,1,3,2};
    cudaMemcpyToSymbol(warp_matrix,host_warp_matrix, sizeof(host_warp_matrix));

    test_kernel<<<1,6>>>();
    cudaDeviceSynchronoze();
}

线程束

SIMD and SIMT

  • SIMD:单指令多数据,多在CPU的流指令集里。
  • SIMT:单指令多线程,CUDA的执行模型

warps和thread blocks的区别

任意时刻,任务的执行都是以warp为单位的,$size_{warp} $ 的大小通常是32,也就是启动$n$个线程,则需要 $ceil(n/size_{warp})$个线程束。

一个block需要的warp数量=$ceil(T/size_{warp})$,其中$T = block中的线程数$ .

对于$gradDim=m, blockDim=n$的任务,总共需要warp数为: $m*ceil(n/size_{warp})$

而实际上物理的warp数量有限,则GPU调度器会根据任务总数分别对warp进行调度。
选中条件:32个core空闲,并且指令参数准备就绪。

block内的warps共用shared memory,通常block内的所有warps需要一定程度上的同步与并行。

grid内的blocks,通常是独立、并行的。每个block之间使用相互独立的shared memory。

因此上级gridDim,blockDim时,考虑block内通常可以互通数据,并行,同步数据等。而block间相互独立。 block的大小也尽量是$size_{warp}$的倍数,尽量避免$size_{warp}$无法填满造成的inactive线程消耗。

warp的执行

在没有线程束分化时,一个warp内的32个线程,可以保证每个线程执行的同一行代码。一个block内的n个线程(比如512),此时512个线程是否执行同一行 代码不确定。

如果需要保证block内的每个线程都到特定行后再执行,在核函数中使用__syncthreaads().

性能优化总结

  1. 尽量warp内线程访问的内存地址是连续的
  2. 尽量少的使用分支(if, switch),造成部分线程inactive
  3. block的大小设置应当是$size_{warp}$的整数倍。
  4. block不宜太小,一般为256,512
  5. 尽量使用shared memory做缓存,避免频繁的与global memory交互
  6. pinned memory应该是内存复制到device的媒介,避免GPU直接访问
  7. 善用constant memory,对于常量性质的媒介,可以利用并加速
  8. 尽量使用pinned memory,而非pageable memory
  9. 当需简化代码内存管理时,可以使用cudaMallocManaged

阿姆达尔定律

Amdahl’s Law(加速比): $$ Speedup = \frac{1}{r_{s}+\frac{r_{p}}{N}} = \frac{t_{p}+t_{s}}{ \frac{t_{p}}{N} + t_{s}}$$

$r_{s}$不可加速部分占比,$r_{p}$可加速部分占比.$N$加速倍数

从公式中我们可以发现,程序中不可加速的部分是限制加速比的瓶颈所在

Licensed under CC BY-NC-SA 4.0