CUDA编程模型

[TOC]

CUDA编程模型

数据传输有cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
常量有cudaMemcpyHostToDevice，cudaMemcpyDeviceToHost，cudaMemcpyDeviceToDevice。

CUDA 编程模型– 基本概念和数据类型

CUDA – C 程序

整合 host+device 两端执行的 C 程序，异构程序。

串行或者并行度不高的部分放在host 端
并行度高的部分放在device端
通过 SPMD kernel C 代码执行

CUDA 设备和线程(Threads)

一个计算设备
（1）对于host或者CPU来说是一个协作处理器
（2）具有自己 DRAM (设备内存)
（3）并行运行许多个线程
（4）通常典型来说是针对于GPU，但是也可以是其它类型的并行处理设备
一个应用的数据并行部分可以在CUDA设备通过kernels调用许多个线程来执行
GPU 和CPU 线程之间的差异
GPU 线程是相当轻量级的
非常小的创建开销（由硬件负责）
GPU 通过需要至少创建1000个线程以上来达到高效率
多核 CPU 仅仅需要创建几个线程

CUDA扩展的C

声明规范

global, device, shared, local, constant

__device__ float filter[N]; 

__global__ void convolve (float *image)  {

__shared__ float region[M];//共享内存变量，只针对一个block内的共享

关键字（内建变量）

threadIdx, blockIdx

1	region[threadIdx] = image[i];

固有的

__syncthreads

1	__syncthreads()

运行时API

Memory, symbol, execution management

1 2	// Allocate GPU memory void *myimage = cudaMalloc(bytes)

函数调用

1 2	// 100 blocks, 10 threads per block convolve<<<100, 10>>> (myimage);

并行的线程数组

一个 CUDA kernel的执行是通过一个线程数组实现的。
（1）所有的线程执行相同的代码(SPMD)
（2）每一个线程有自己的ID，通过线程ID来计算访问内存地址和做出控制决定

线程块(Thread Block): 可扩展的合作

把整体的线程数组划分到多个块中
（1）一个Block中的Threads通过shared memory, atomic operations 和 barrier synchronization来进行合作
（2）不同Block中的Threads 不能合作

Block IDs 和 Thread IDs

每一个thread用IDs 来决定对那块数据进行操作
（1）Block ID: 1D 或者 2D
（2）Thread ID: 1D, 2D, 或者3D

这样对于多维数据，可以简化其内存访问模式
（1）图像处理
（2）在体元范围内解PDE方程…

CUDA 内存模型概述※

全局内存
（1）host 和 device端主要的数据R/W交流手段
（2）其内容对于所有线程都是可见的
（3）较长的访问延迟
当前，首先关注全局内存的特点

CUDA 应用程序编程接口API-基本

这个 API 是对于 ANSI C 程序语言的一个扩展
低学习曲线
其专门设计的硬件可以启用轻量级的运行时API以及驱动
高性能

CUDA 设备内存分配

cudaMalloc()

在设备端的 Global Memory 上分配内存
需要两个参数
（1）指向待分配对象的指针的地址
（2）待分配对象的规模大小

cudaFree()

从设备端的Global Memory释放内存
（1）指向待释放对象的指针

代码示例

分配一个64单精度浮点数组
将分配空间绑定到aD
“D” 通常用来指示一个设备端的数据或数据结构

N = 64;
float *aD = NULL;
int size = N * sizeof(float);

cudaMalloc((void**)&aD, size);
cudaFree(Md);

CUDA Host-Device数据传输

cudaMemcpy()

内存数据传输
需要四个参数
（1）指向目标数据的指针
（2）指向源数据的指针
（3）需要拷贝的数据量(bytes)
（4）传输类型（方向）
- Host to Host
- Host to Device
- Device to Host
- Device to Device

异步传输

代码示例

传输一个64 单精度浮点数组
a位于host 内存， aD 位于device内存
cudaMemcpyHostToDevice 和 cudaMemcpyDeviceToHost 是符号常量

1 2	cudaMemcpy(aD, a, size, cudaMemcpyHostToDevice); cudaMemcpy(a, aD, size, cudaMemcpyDeviceToHost);

CUDA 函数声明

	Executed on the:	Only callable from the:
`__device__ float DeviceFunc()`	device	devic
`__global__ void KernelFunc()`	device	host
`__host__ float HostFunc()`	host	host

__global__ 定义一个 kernel 函数必须返回void
__device__ 和 __host__ 能一起使用，编译器会把这个函数编译为host和device通用的函数
__device__ 函数没有函数地址，也没有指向它的函数指针
在device端执行的函数有下面的限制
（1）没有递归
（2）函数内部没有静态变量
（3）参数的数量是固定的

调用Kernel函数 – 线程创建

一个kernel函数在调用之前，必须在执行之前进行线程配置:

__global__ void KernelFunc(...);
dim3   DimGrid(100, 50);    // 5000 thread blocks 
dim3   DimBlock(4, 8, 8);   // 256 threads per block 
size_t SharedMemBytes = 64; // 64 bytes of shared memory
KernelFunc<<< DimGrid, DimBlock, SharedMemBytes >>>(...);

一个SM有64K，
1个block，SharedMemBytes=64K
2个block，SharedMemBytes=32K
4个block，SharedMemBytes=16K

NVIDA最新卡的信息

Kernel函数是支持异步执行的

一个简单的可运行的例子-矢量相加

一个简单的矢量相加的例子演示了在CUDA程序中的内存基本特性和线程管理方法
（1）shared memory 的使用将放在后面的章节
（2）局部及寄存器的使用
（3）Thread ID 的使用
（4）host 和device之间内存数据传输API的使用

通用代码↑

简单版本代码↑

CUDA Syntax – Unified Memory 统一内存(不建议使用)

CUDA 6 introduced Unified Memory
（1）Use managed memory instead of explicitly declaring memory for the host and device
代替host和device内存概念
（2）cudaMallocManaged(void **devPtr, size_t size)

矢量相加Unified Memory实现

在CPU和GPU都创建内存函数

仅仅使用一个Thread Block处理矢量

一个thread block来计算矢量相加
（1）每一个线程计算矢量中的一个元素
每一个线程
（1）装载矢量a的一个数据
（2）装载矢量b的一个数据
（3）对于每一对a和b中的元素执行一次加法操作
同时矢量的尺寸受限于一个thread block中规定的线程数

对于c=a+b，a有四个线程，b有四个线程，运行受限制，这个就要考虑多个block联用。

处理任意大小的矢量

__global__ void SumArrays(float* const a, float* const b, float* const c, int const N)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    if (i < N)
        c[i] = a[i] + b[i];
}

int blkNum = ((N + BLOCK_SIZE - 1) / BLOCK_SIZE);	//设置block数量
SumArrays<<<blkNum, BLOCK_SIZE>>>(aD, bD, cD, N);	//

写一个矩阵相加的并行计算程序，矩阵大小为35 * 35，请用任意数初始化两个矩阵A,B.

void MatAddonDevice(float * aH,float* bH, float *cH, int N)  // {
   float* aD, * bD, * cD;
   int N_BYTES = N * sizeof(float);
   
   // 分配内存
   cudaMalloc((void**)&aD, N_BYTES);
   cudaMalloc((void**)&bD, N_BYTES);
   cudaMalloc((void**)&cD, N_BYTES);
   
   // host -> Device
   cudaMemcpy(aD, aH, N_BYTES, cudaMemcpyHostToDevice);
   cudaMemcpy(bD, bH, N_BYTES, cudaMemcpyHostToDevice);
   
   dim3 threadsPerblock(16, 16);
   dim3 blocksPerGrid((N - 1) / 16 + 1,(N - 1) / 16 + 1);
   
   vectorAdd<<<blocksPerGrid, threadsPerblock>>>(aD, bD, cD, N);
   
   cudaMemcpy(cH, cD, size, cudaMemcpyDeviceToHost);
   
   cudaFree(aD);
   cudaFree(bD);
   cudaFree(cD);
   
   free(aH);
   free(bH);
   free(cH);
}

___global___ void MatAddKernel(float *a, float*b, float *c, int N)
{
   int col = threadIdx.x + blockIdx.x * blockDim.x;
   int row = threadIdx.y + blockIdx.y * blockDim.y;
   if (Col < WIDTH && Row < HEIGHT)
	{
		c[Row * N + Col] = a[Row * N + Col] + b[Row * N + Col];
	}
}

任意可执行的CUDA代码需要两个动态库

（1）CUDA 运行时库 (cudart)
（2）CUDA 核心库 (cuda)

查看机器性能

执行线程块（Thread Blocks）

在Block层面上的资源分配

grid->blocks->threads(256)

在Warp层面上的调度形式

数据存储读取问题

SIMD执行的 SIMT实现

单指令多数据-SIMD在CUDA中是单指令多线程-SIMT

在`block`粒度上，线程是被分配到流多处理器(`Streaming Multiprocessors`)[特定架构]

每个SM最多运行放8 blocks
每个Fermi SM能接收1536 threads
可以装 256 (threads/block) 6 blocks
或者 512 (threads/block) 3 blocks，等。

线程的执行是并发同时执行

SM 维护thread/block id
SM 管理/调度线程执行

线程调度

每个Block的执行如32-thread Warps

Warps非CUDA编程模型的一部分，是一个硬件层的设施。
grid/block/thread是编程模型。
Warps是SM中的调度单元。

如果3 Blocks被分配到一个SM中，每个Block有 256 threads, 那么在1个SM有多少个Warps?

每个Block被划分成 256 / 32 = 8 Warps。
这里将会有 8 * 3 = 24 Warps 。

线程块如何切分

Thread blocks 被切进warps

在一个Warp中的Thread IDs 是连续递增的。
Warp 0开始于Thread ID 0。

切分总是相同的

在控制流中，需要用到这部分内容。
然而，下一代的硬件设备可能改变Warps的规模大小。
一些内容后面将涉及。

然而，warps之间不存在任何顺序的依赖

如果线程之间出现任何依赖关系，必须用__syncthreads()来得到正确的结果。

控制流指令

影响并行程序的性能主要考量的分支是发散的

在单一的warp中Threads出现了不同的分支路径。
不同的执行路径在当前的GPUs中的执行时串行的。
在一个warp中Threads所采用的控制路径将会遍历执行，直到再没有更多的路径。

通常处理方法: 当这个分支条件是thread ID的函数的时候，要避免这种发散情况的发生

出现发散情况的例子

1	If (threadIdx.x > 2) { }

这里将会在一个block 中，为所有线程创建两条不同的控制路径。
分支间隔尺寸< warp size；在第一个warp中，相对于其它线程，threads 0, 1 和 2 将会遵循不同的路径。

没有发散情况的例子:

1	If (threadIdx.x / WARP_SIZE > 2) { }

同样为所有线程创建两条不同的控制路径。
分支间隔尺寸是整个warp的大小；在第一个warp中，所有的线程将会遵循相同的路径。

块尺寸考虑

对于矩阵相乘，需要要多个Blocks，那么到底用8X8, 16X16 还是 32X32 的块尺寸?

对于 8 × 8，每一个Block有 64 threads。因为每个SM能接收1536 threads，这里就会有24 Blocks。然而，每一个SM只能接收 8 Blocks，那么仅仅 512 threads 将会进入到每个SM，仅仅能利用到 1 / 3 的SM的线程能力。
对于 16 × 16，每一个Block有 256 threads。因为每个SM能接收 1536 threads，这里就会有 6 Blocks，刚好能获得最大的SM利用率。
对于 32 × 32，每一个Block有 1024 threads。因为每个SM能接收 1536 threads，这里仅仅就只能用到 2 / 3 的SM的线程能力。

课堂练习

已知费米架构（英伟达显卡架构之一）下每个SM能接收 1536 threads，每一个SM能接收 8 Blocks。求下列线程块配置时SM的占用率（利用率）。

(1) int threadPerBlock = 128;

128 * 8 = 1024
1024 / 1536 = 66.67%

(2) int threadPerBlock = 256;

256 * 6 = 1536
1536 / 1536 = 100%

(3) int threadPerBlock = 512;

512 * 3 = 1536
1536 / 1536 = 100%

(4) int threadPerBlock = 1024；

1024 * 1 = 1024
1024 / 1536 = 66.67%

费米架构：英伟达显卡架构之一，不同的架构下，每个SM的性能是不一样的。

多维线程应用：方矩阵相乘

P = M N，尺寸为 WIDTH HEIGHT；
并行计算可以采用一个线程计算 P 中的一个元素。
M 与 N 将会从 global memory 被装载 WIDTH 次

并行计算不能有数据依赖，采用不同线程计算P中间的不同的元素。

矩阵相乘CPU端代码

矩阵相乘时间复杂度为O^3，对应位置相乘最后相加。

// Matrix multiplication on the (CPU) host in double precision，N*N
void MatrixMulOnHost(float* M, float* N, float* P, int Width)
{
    for (int i = 0; i < Width; ++i)
        for (int j = 0; j < Width; ++j) {
            double sum = 0;
            for (int k = 0; k < Width; ++k) {
                double a = M[i * width + k];
                double b = N[k * width + j];
                sum += a * b;
            }
            P[i * Width + j] = sum;
        }
}

输入矩阵数据传输(Host端代码)

void MatrixMulOnDevice(float* M, float* N, float* P, int Width)
{
    int size = Width * Width * sizeof(float); 
    float* Md, Nd, Pd;
    …
//1.//  Allocate and Load M, N to device memory 
    cudaMalloc(&Md, size);
    cudaMemcpy(Md, M, size, cudaMemcpyHostToDevice);

    cudaMalloc(&Nd, size);
    cudaMemcpy(Nd, N, size, cudaMemcpyHostToDevice);

    // Allocate P on the device
    cudaMalloc(&Pd, size);

//2.// Kernel invocation code – 后面显示调用过程
     …
     
//3.// Read P from the device
    cudaMemcpy(P, Pd, size, cudaMemcpyDeviceToHost);

    // Free device matrices
    cudaFree(Md);
    cudaFree(Nd);
    cudaFree(Pd);
}

Kernel函数

// Matrix multiplication kernel – per thread code
__global__ void MatrixMulKernel(float* Md, float* Nd, float* Pd, int Width)
{
    // 2D Thread ID
    int tx = threadIdx.x;	//对应列 col
    int ty = threadIdx.y;	//对应行 row

    // Pvalue 用来存储结果矩阵的元素
    // 通过每一个线程来计算结果矩阵的每一个元素
    float Pvalue = 0;
    
    //用线程去计算
    for (int k = 0; k < Width; ++k)
    { 
         float Melement = Md[ty * Width + k];
         float Nelement = Nd[k * Width + tx];
         Pvalue += Melement * Nelement;
    }
    // Write the matrix to device memory;
    // each thread writes one element;
    Pd[ty * Width + tx] = Pvalue;
}

Kernel函数调用(host-side Code)

// 线程配置
dim3 dimBlock(Width, Width);
dim3 dimGrid(1, 1);

// 开始在device端执行该Kernel函数!
MatrixMulKernel<<<dimGrid, dimBlock>>>(Md, Nd, Pd);

仅适用一个Thread Block

一个thread block来计算矩阵 Pd

每一个线程计算Pd中的一个元素

每一个线程

装载矩阵Md的一行数据
装载矩阵Nd的一列数据
对于每一对Md和Nd中的元素执行一次乘法和加法操作
计算与内存访问比例接近1:1 (不是非常高)

同时矩阵的尺寸也受限于一个thread block中规定的线程数

处理任意大小的方阵

每一个 2D thread block 计算结果矩阵中一个大小为 (TILE_WIDTH)^2 的子矩阵 (瓷砖化)，每一个都有 (TILE_WIDTH)^2 threads。
产生一个2D GRID，其具有 (WIDTH/TILE_WIDTH)^2 blocks。

如果 WIDTH/TILE_WIDTH 大于最大的grid大小 (64K)!，需要用多个循环来处理。

一个稍大些的例子

WIDTH = 8；TILE_WIDTH = 2；
每一个block有 2 * 2 = 4 threads

WIDTH / TILE_WIDTH = 4
使用 4* 4 = 16 blocks

一个再大一些的例子

WIDTH = 8；TILE_WIDTH = 4；
每一个block有 4 * 4 = 16 threads

WIDTH / TILE_WIDTH = 2
使用 2* 2 = 4 blocks

C/C++的“行主”分布

Blcok(0,0)

Blcok(0,1)

一个简单的矩阵相乘Kernel

__global__ void MatrixMulKernel(float* d_M, float* d_N, float* d_P, int Width)
{
    // 计算输入数据d_M和输出数据d_P的行号以确定处理的元素
    int Row = blockIdx.y * blockDim.y + threadIdx.y;
    // 计算输入数据d_N和输出数据d_P的行号以确定处理的元素
    int Col = blockIdx.x * blockDim.x + threadIdx.x;

    if ((Row < Width) && (Col < Width)) {
        float Pvalue = 0;
        // 每一个线程处理输出数据d_P的一个元素
        for (int k = 0; k < Width; ++k)
            Pvalue += d_M[Row * Width + k] * d_N[k * Width + Col];
        d_P[Row * Width + Col] = Pvalue;
    }
}

高性能计算第4次课-PartA.mp4
Lecture3B_CUDA编程模型.ppt 15/29

CUDA编程模型

CUDA 编程模型– 基本概念和数据类型

CUDA – C 程序

CUDA 设备和线程(Threads)

CUDA扩展的C

声明规范

关键字（内建变量）

固有的

运行时API

函数调用

并行的线程数组

线程块(Thread Block): 可扩展的合作

Block IDs 和 Thread IDs

CUDA 内存模型概述※

CUDA 应用程序编程接口API-基本

CUDA 设备内存分配

cudaMalloc()

cudaFree()

代码示例

CUDA Host-Device数据传输

cudaMemcpy()

异步传输

代码示例

CUDA 函数声明

调用Kernel函数 – 线程创建

一个简单的可运行的例子-矢量相加

CUDA Syntax – Unified Memory 统一内存(不建议使用)

矢量相加Unified Memory实现

仅仅使用一个Thread Block处理矢量

处理任意大小的矢量

写一个矩阵相加的并行计算程序，矩阵大小为35 * 35，请用任意数初始化两个矩阵A,B.

任意可执行的CUDA代码需要两个动态库

查看机器性能

执行线程块（Thread Blocks）

在block粒度上，线程是被分配到流多处理器(Streaming Multiprocessors)[特定架构]

线程的执行是并发同时执行

线程调度

每个Block的执行如32-thread Warps

如果3 Blocks被分配到一个SM中，每个Block有 256 threads, 那么在1个SM有多少个Warps?

线程块如何切分

Thread blocks 被切进warps

切分总是相同的

然而，warps之间不存在任何顺序的依赖

控制流指令

影响并行程序的性能主要考量的分支是发散的

通常处理方法: 当这个分支条件是thread ID的函数的时候，要避免这种发散情况的发生

块尺寸考虑

对于矩阵相乘，需要要多个Blocks，那么到底用8X8, 16X16 还是 32X32 的块尺寸?

课堂练习

多维线程应用：方矩阵相乘

矩阵相乘CPU端代码

输入矩阵数据传输(Host端代码)

Kernel函数

Kernel函数调用(host-side Code)

仅适用一个Thread Block

处理任意大小的方阵

一个稍大些的例子

一个再大一些的例子

C/C++的“行主”分布

Blcok(0,0)

Blcok(0,1)

一个简单的矩阵相乘Kernel

在`block`粒度上，线程是被分配到流多处理器(`Streaming Multiprocessors`)[特定架构]