2024/05/12
Professional CUDA C Programming
Chapter02 CUDA Programming Model
2.1 CUDA 编程模型概述
CUDA编程模型提供了一个计算机架构抽象作为应用程序和其可用硬件之间的桥梁。
CUDA编程模型利用GPU架构的计算能力提供的几个特有功能
- 一种通过层次结构在GPU中组织线程的方法
- 一种通过层次结构在GPU中组织内存的方法
CUDA 编程结构
主机:CPU及其内存(主机内存)
设备:GPU及其内存(设备内存)
代码规范
主机内存中的变量名以h_为前缀,设备内存中的变量名以d_为前缀。
一个典型的CUDA程序实现流程遵循以下模式:
- 把数据从CPU内存拷贝到GPU内存;
- 调用核函数对存储在GPU内存中的数据进行操作;
- 将数据从GPU内存传送回CPU内存;
内存管理
表 2-1 主机和设备内存函数
标准的C函数 | CUDA C函数 | 标准的C函数 | CUDA C函数 |
---|---|---|---|
malloc | cudaMalloc | memset | cudaMemset |
memcpy | cudaMemcpy | free | cudaFree |
用于执行GPU内存分配的cudaMalloc函数,其函数原型为:
cudaError_t cudaMalloc(void** devPtr,size_t size)
该函数负责向设备分配一定字节的线性内存,并以devPtr的形式返回指向所分配内存的指针。
cudaMemcpy函数负责主机和设备之间的数据传输,其函数原型为:
cudaError_t cudaMemcpy(void* dst,void* src,size_t count,cudaMemcpyKind lind)
此函数从src指向的源存储区复制一定数量的字节到dst指向的目标存储区。复制方向由kind执行,其中kind有以下几种
cudaMemcpyHostToHost
cudaMemcpyHostToDevice
cudaMemcpyDeviceToHost
cudaMemcpyDeviceToDevice
这个函数以同步方式执行,因为在cudaMemcpuy函数返回以及完成传输操作之前主机应用程序是阻塞的。除了内核启动之外的CUDA调用都会返回一个错误的枚举类型cudaError_t。如果GPU内存分配成功。函数返回:
cudaSuccess
否则返回
cudaErrorMemoryAllocation
可使用以下CUDA运行时函数将错误信息转化为可读的错误信息:
char* cudaGetErrorString(cudaError_t error)
内存层次结构
CUDA编程模型最显著的一个特点就是揭示了内存层次结构。每一个GPU设备都有用于不同用途的存储类型。
在GPU内存层次中,最主要的两种内存是全局内存和共享内存。全局类似于CPU的系统内存,而共享内存类似于CPU的缓存。然而GPU的共享内存可以由CUDA C的内核直接控制。
//第2章代码公用函数头文件
//ch02_Header.cuh
#if !defined __CH02_HEADER_H__
#define __CH02_HEADER_H__
#include <cuda_runtime.h>
#include <time.h>
#include <stdlib.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); \
} \
}
static clock_t cpuSecond() {
return clock();
}
static void initialData(float* p_A, const int p_size) {
for (size_t i = 0; i < p_size; i++)
{
p_A[i] = i;
}
}
static void printResult(float *p_A,float *p_B, float*p_hostRef,float* p_gpuRef, const int p_size) {
for (size_t i = 0; i < p_size; i++)
{
printf("A %5f B %5f C %5f SUM %5f \n", p_A[i], p_B[i], p_hostRef[i], p_gpuRef[i]);
}
}
static void initialData_random(float* p_ip, int p_size) {
time_t t;
srand((unsigned int)time(&t));
for (size_t i = 0; i < p_size; i++)
{
p_ip[i] = (float)(rand() & 0xFF) / 10.0F;
}
}
static void initilaInt(int* p_ip, int p_size) {
for (size_t i = 0; i < p_size; i++)
{
p_ip[i] = i;
}
}
static void sumArrayOnHost(float* p_A, float* p_B, float* p_C, const int p_N) {
for (size_t i = 0; i < p_N; i++)
{
p_C[i] = p_A[i] + p_B[i];
}
}
static void checkResult(float* p_hostRef, float* p_gpuRef, const int p_N) {
double epsilon = 1.0E-8;
int match = 1;
for (size_t i = 0; i < p_N; i++)
{
if (abs(p_hostRef[i] - p_gpuRef[i]) > epsilon) {
match = 0;
printf("Array do not match\n");
printf("host %.2f , gpu %.2f , at current %3d\n", p_hostRef[i], p_gpuRef[i],i);
break;
}
}
if (match)
printf("Array match.\n\n");
return;
}
#endif
//sumArrayOnHost_GPU.cu
//两个数组相加
#include "ch02_Header.cuh"
int invokeKernel();
//int main() {
// return invokeKernel();
//}
__global__ void sumArraysOnGPU(float* p_A, float* p_B, float* p_C, const int p_N) {
int i = blockIdx.x * blockDim.x + threadIdx.x;
if (i < p_N)
p_C[i] = p_A[i] + p_B[i];
printf("%f", p_C[i]);
}
static int invokeKernel() {
printf("%s Starting...\n");
int dev = 0;
cudaDeviceProp deviceProp;
CHECK(cudaGetDeviceProperties(&deviceProp, dev));
printf("Using Device %d: %s\n", dev, deviceProp.name);
CHECK(cudaSetDevice(0));
//不是可以无限大
int nElem = 1 << 24;
printf("Vector size %d\n", nElem);
size_t nBytes = nElem * sizeof(float);
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);
long iStart, iElaps;
iStart = cpuSecond();
printf("iStart %d\n", iStart);
initialData_random(h_A, nElem);
initialData_random(h_B, nElem);
iElaps = cpuSecond() - iStart;
memset(hostRef, 0, nBytes);
memset(gpuRef, 0, nBytes);
iElaps = cpuSecond() - iStart;
printf("memset Time elapsed %d ms\n", iElaps);
iStart = cpuSecond();
sumArrayOnHost(h_A, h_B, hostRef, nElem);
printf("iStart %d\n", cpuSecond());
iElaps = cpuSecond() - iStart;
printf("sumArrayOnHost Time elapsed %20d ms\n", iElaps);
float* d_A, * d_B, * d_C;
//用cudaMalloc在GPU上申请内存
cudaMalloc((float**)&d_A, nBytes);
cudaMalloc((float**)&d_B, nBytes);
cudaMalloc((float**)&d_C, nBytes);
//使用cudaMemcpy函数把数据从主机内存拷贝到GPU的全局内存中,参数cudaMemcpyHostToDevice指定数据拷贝方向。
cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, nBytes, cudaMemcpyHostToDevice);
dim3 block(nElem);
dim3 grid(nElem / block.x);
iStart = cpuSecond();
sumArraysOnGPU << <grid, block >> > (d_A, d_B, d_C, nElem);
iElaps = cpuSecond() - iStart;
printf("sumArraysOnGPU <<<%d, %d>>> Time elapsed %d ms\n", grid.x, block.x, iElaps / CLOCKS_PER_SEC);
printf("Execution configuration <<<%d, %d>>>\n", grid.x, block.x);
//使用cudaMemcpy函数把结果从GPU复制到主机的数组gpuRef中
CHECK(cudaMemcpy(gpuRef, d_C, nBytes, cudaMemcpyDeviceToHost));
checkResult(hostRef, gpuRef, nElem);
printResult(h_A, h_B, hostRef, gpuRef, nElem);
//调用cudaFree释放GPU的内存
cudaFree(d_A);
cudaFree(d_B);
cudaFree(d_C);
free(h_A);
free(h_B);
free(hostRef);
free(gpuRef);
return 0;
}
在Unity3D中用ComputeShader实现的相同功能
using System.Collections;
using System.Collections.Generic;
using UnityEngine;
public class sumArrayOnHost_GPU : MonoBehaviour
{
[SerializeField] ComputeShader computeShader;
float[] h_A;
float[] h_B;
float[] gpuRef;
ComputeBuffer d_A;
ComputeBuffer d_B;
ComputeBuffer d_C;
int sumArraysOnGPU;
// Start is called before the first frame update
void Start()
{
int nElem = 1 << 8;
h_A = new float[nElem];
h_B = new float[nElem];
initData(h_A, nElem);
initData(h_B, nElem);
gpuRef = new float[nElem];
d_A = new ComputeBuffer(nElem, sizeof(float));
d_B = new ComputeBuffer(nElem, sizeof(float));
d_C = new ComputeBuffer(nElem, sizeof(float));
computeShader.SetInt("N", nElem);
d_A.SetData(h_A);
d_B.SetData(h_B);
//d_C.SetData(gpuRef);
sumArraysOnGPU = computeShader.FindKernel("sumArraysOnGPU");
computeShader.SetBuffer(sumArraysOnGPU, "d_A", d_A);
computeShader.SetBuffer(sumArraysOnGPU, "d_B", d_B);
computeShader.SetBuffer(sumArraysOnGPU, "d_C", d_C);
int grid_x = nElem / 2;
Debug.Log(grid_x);
computeShader.Dispatch(sumArraysOnGPU, grid_x, 2, 3);
d_C.GetData(gpuRef);
for (int i = 0; i < gpuRef.Length; i++)
{
Debug.Log(string.Format($"A :{h_A[i]} + B {h_B[i]} = C {gpuRef[i]}"));
}
}
void initData(float[] p_f, int p_size)
{
for (int i = 0; i < p_size; i++)
{
p_f[i] = i;
}
}
// Update is called once per frame
void Update()
{
}
private void OnApplicationQuit()
{
d_A.Dispose();
d_B.Dispose();
d_C.Dispose();
}
}
// Each #kernel tells which function to compile; you can have many kernels
#pragma kernel sumArraysOnGPU
#pragma enable_d3d11_debug_symbols
// Create a RenderTexture with enableRandomWrite flag and set it
// with cs.SetTexture
int N;
StructuredBuffer<float> d_A;
StructuredBuffer<float> d_B;
RWStructuredBuffer<float> d_C;
[numthreads(8,1,1)]
void sumArraysOnGPU(uint3 id : SV_DispatchThreadID,uint3 block:SV_GroupID)
{
if(id.x<N)
d_C[id.z]=d_A[id.x] + d_B[id.x];
}
对于其中的参数SV_DispatchThreadID,SV_GroupID的含义在学习后会进一步说明,目前的问题是SV_DispathThreadID貌似也是线性展开的。
2024/05/15
线程管理
由一个内核启动所产生的所有线程称为一个网格。同一网格中的所有线程共享相同的全局内存空间。一个网格由多个线程块构成,一个线程块包含一组线程,同一线程块内的线程协作可以通过以下方式来实现。
- 同步
- 共享内存
不同块内的线程不能协作
线程依靠以下两个坐标变量来区分彼此。
blockIdx(线程块在线程网格内的索引)
threadIdx(线程在线程块内的索引)
CUDA可以组织三维的网格和块。
网格和块的维度由下列两个内置变量指定。
blockDim(线程块的维度,用每个线程块中的线程数来表示)
gridDim(线程网格的维度,用每个线程网格中的线程块数来表示(中文版似乎是翻译错了))
//checkDimension.cu
//检查网格和块的索引和维度
#include<cuda_runtime.h>
#include<stdio.h>
__global__ void checkIndex(void) {
printf("threadIdx:(%d, %d, %d) blockIdx:(%d, %d, %d) blockDim:(%d, %d, %d) ,gridDim(%d, %d, %d)\n", threadIdx.x, threadIdx.y, threadIdx.z, blockIdx.x, blockIdx.y, blockIdx.z, blockDim.x, blockDim.y, blockDim.z, gridDim.x, gridDim.y, gridDim.z);
}
int main() {
int nElem = 7;
dim3 block(3);
//功能类似于向上取整
dim3 grid((nElem + block.x - 1) / block.x);
printf("grid.x %d, grid.y %d, grid.z %d)\n", grid.x, grid.y, grid.z);
printf("block.x %d, block.y %d, block.z %d)\n", block.x, block.y, block.z);
checkIndex << <grid, block >> > ();
cudaDeviceReset();
return 0;
}
输出
grid.x 3, grid.y 1, grid.z 1)
block.x 3, block.y 1, block.z 1)
threadIdx:(0, 0, 0) blockIdx:(0, 0, 0) blockDim:(3, 1, 1) ,gridDim(3, 1, 1)
threadIdx:(1, 0, 0) blockIdx:(0, 0, 0) blockDim:(3, 1, 1) ,gridDim(3, 1, 1)
threadIdx:(2, 0, 0) blockIdx:(0, 0, 0) blockDim:(3, 1, 1) ,gridDim(3, 1, 1)
threadIdx:(0, 0, 0) blockIdx:(1, 0, 0) blockDim:(3, 1, 1) ,gridDim(3, 1, 1)
threadIdx:(1, 0, 0) blockIdx:(1, 0, 0) blockDim:(3, 1, 1) ,gridDim(3, 1, 1)
threadIdx:(2, 0, 0) blockIdx:(1, 0, 0) blockDim:(3, 1, 1) ,gridDim(3, 1, 1)
threadIdx:(0, 0, 0) blockIdx:(2, 0, 0) blockDim:(3, 1, 1) ,gridDim(3, 1, 1)
threadIdx:(1, 0, 0) blockIdx:(2, 0, 0) blockDim:(3, 1, 1) ,gridDim(3, 1, 1)
threadIdx:(2, 0, 0) blockIdx:(2, 0, 0) blockDim:(3, 1, 1) ,gridDim(3, 1, 1)
从主机端和设备端访问网格/块变量
区分主机端和设备端的网格和块变量的访问是很重要的。例如,声明一个主机端的块变量,你按如下定义它的坐标并对其进行访问:
block.x,block.y,block.z
在设备端,你已经预定义了内置块变量的大小:
blockDim.x,blockDim.y,blockDim.z
总之,在启动内核之前就定义了主机端的网格和块变量,并从主机端通过由x,y,z三个字段决定的矢量结构来访问它们。当内核启动时,可以使用内核中预初始化的内置变量。
对于一个给定的数据大小,确定网格和块尺寸的一般步骤为:
- 确定块的大小
- 在已知数据大小和块大小的基础上计算网格维度
要确定块尺寸,通常需要考虑
- 内核的性能特性
- CPU资源的限制
#include <cuda_runtime.h>
#include <stdio.h>
int main() {
//define total data elements
int nElem = 1024;
//define grid and block structure
dim3 block(1024);
dim3 grid((nElem + block.x - 1) / block.x);
//reset block
printf("grid.x %d block.x %d \n", grid.x, block.x);
block.x = 512;
grid.x = ((nElem + block.x - 1) / block.x);
printf("grid.x %d block.x %d \n", grid.x, block.x);
//reset block
block.x = 256;
grid.x = ((nElem + block.x - 1) / block.x);
printf("grid.x %d block.x %d \n", grid.x, block.x);
//reset block
block.x = 128;
grid.x = ((nElem + block.x - 1) / block.x);
printf("grid.x %d block.x %d \n", grid.x, block.x);
cudaDeviceReset();
}
输出
grid.x 1 block.x 1024
grid.x 2 block.x 512
grid.x 4 block.x 256
grid.x 8 block.x 128
线程层次结构
CUDA的特点之一就是通过编程模型揭示了一个两层的线程层次结构。由于一个内核启动的网格和块的维数会影响性能,这一结构为程序员优化程序提供了一个额外的途径。
网格和块的维度存在几个限制因素,对于块大小的一个主要限制因素就是可利用的计算资源,如寄存器,共享内存等。Some limits can be retrieved by querying the GPU device.
网格和块从逻辑上代表了一个核函数的线程层次结构,在第3章中,你会发现这种线程组织方式能使你在不同的设备上有效执行相同的程序代码,而且每一个线程组织具有不同数量的计算和内存资源。
启动一个CUDA核函数
CUDA内核调用是对C语言函数调用语句的延伸,<<<>>>运算符内是核函数的执行配置。
kernel_name <<<grid,block>>>(argument list);
CUDA编程模型揭示了线程层次结构。利用执行配置可以指定线程在GPU上调度运行的方式,执行配置的第一个值是网格维度,也就是启块的数目。第二个值是块维度,也就是每个块中线程的数目。通过指定网格和块的维度,你可以进行以下配置:
- 内核中的线程的数目
- 内核中使用的线程布局
同一个块中的线程之间可以相互协作,不同块内的线程不能协作。对于一个给定的问题,可以使用不同的网格和块布局来组织你的线程。
由于数据在全局内存中是线性存储的,因此可以用变量blockIdx.x和threadIdx来进行操作。
- 在网格中标识一个唯一的线程
- 建立线程和数据元素之间的映射关系
核函数的调用与主机线程是异步的。核函数调用结束后,控制权立刻返回给主机端。
你可以调用以下函数来强制主机端程序等待所有的核函数执行结束:
cudaError_t cudaDiveceSynchornize(void);
异步行为
不同于C语言的函数调用,所有的CUDA核函数的启动都是异步的。CUDA内核调用完成后,控制权立刻返回给CPU。
2024/05/16
编写核函数
核函数是在设备端执行的代码。在核函数中,需要为一个线程规定要进行的计算以及要进行的数据访问。当核函数被调用时,许多不同的CUDA线程并行执行同一个计算任务。
以下是__global__声明定义该函数:
__gloabal__ void kernel_name(argument list);
表 2-2 函数类型限定符
限定符 | 执行 | 调用 | 备注 |
---|---|---|---|
__global__ | 在设备端执行 | 可从主机端调用 | 必须有一个void返回类型 |
__device__ | 在设备端执行 | 仅能从设备端调用 | |
__host__ | 在主机端执行 | 仅能从主机端调用 | 可以省略 |
__device__和 __host__限定符可以一起使用,这样函数可以同时在主机和设备端进行编译。
CUDA核函数的限制
以下限制使用与所有核函数:
- 只能访问设备内存
- 必须具有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(cuda_Synchronize());
给核函数计时
用CPU计时器计时
#include <time.h>
static clock_t cpuSecond() {
return clock();
}
clock_t iStart=cpuSecond();
kernel_name<<<grid,block>>>(argument list);
cudaDeviceSynchronize();
//由于核函数调用与主机端程序是异步的,你需要用cudaDeviceSynchronize函数来等待所有的GPU线程运行结束。
clock_t iElaps=cpuSecond()-iStart;
了解自身局限性
在调整执行配置时需要了解一个关键点是对网格和块维度的限制。线程层次结构中每个层级的最大尺寸取决于设备。
CUDA提供了通过查询GPU来了解这些限制的能力。在本章的2.4节有详细的介绍。
对于Fermi设备,每个块的最大线程数是1024,且网格的x、y、z三个方向上的维度的最大值是65535
用nvprof工具计时
$ nvprof ./sumArrayOnGPU
nvprof E:\C_CPP_CUDA\CUDA\x64\Release\CUDATest.exe
组织并行线程
NOTES:Release模式下,书中的几种不同的布局并未得到不同的性能表现,可能是现在GPU架构对代码优化更为强大了吧,通过本节的学习可以掌握如何使用线程块和线程与数据索引之间建立联系
使用块和线程建立矩阵索引
//checkThreadIndex.cu
//检查块和线程索引
#include "ch02_Header.cuh"
int invokeKernel();
int main() {
invokeKernel();
}
void printMatrix(int* p_C, const int p_nx, const int p_ny) {
int* ic = p_C;
printf("\nMatrix: (%d.%d)\n ", p_nx, p_ny);
for (size_t iy = 0; iy < p_ny; iy++)
{
for (size_t ix = 0; ix < p_nx; ix++)
{
printf("%3d", ic[ix]);
}
ic += p_nx;
printf("\n");
}
printf("\n");
}
__global__ void printThreadIndex(int* p_A, const int p_nx, const int p_ny) {
int ix = threadIdx.x + blockIdx.x * blockDim.x;
int iy = threadIdx.y + blockIdx.y * blockDim.y;
unsigned int idx = ix * p_nx + ix;
printf("thread_id (%d,%d) block_id (%d,%d) coordinate (%d,%d)"
"global index %3d ival %3d\n", threadIdx.x, threadIdx.y, blockIdx.x, blockIdx.y, ix, iy, idx, p_A[idx]);
}
static int invokeKernel() {
printf("%s Starting ...\n");
//get device information
int dev = 0;
cudaDeviceProp deviceProp;
CHECK(cudaGetDeviceProperties(&deviceProp, dev));
printf("Using Device %d: %s\n", dev, deviceProp);
CHECK(cudaSetDevice(dev));
//set matrix dimension
int nx = 8;
int ny = 6;
int nxy = nx * ny;
int nBytes = nxy * sizeof(float);
//malloc host memory
int* h_A;
h_A = (int*)malloc(nBytes);
//initialize host matrix with integer
initilaInt(h_A, nxy);
printMatrix(h_A, nx, ny);
int* d_MatA;
//malloc device memory
cudaMalloc((void**)&d_MatA, nBytes);
//transfer data from host to device
cudaMemcpy(d_MatA, h_A, nBytes, cudaMemcpyHostToDevice);
//set up execution configuration
dim3 block(4, 2);
dim3 grid((nx + block.x - 1) / block.x, (ny + block.y - 1) / block.y);
//invoke the kernel
printThreadIndex << <grid, block >> > (d_MatA, nx, ny);
//free host and device memory
cudaFree(d_MatA);
free(h_A);
//reset device
cudaDeviceReset();
return 0;
}
使用二维网格