CUDA 学习笔记

推荐书目：

CUDA by Example. An Introduction to General-Purpose GPU Programming (Jason Sanders, Edward Kandrot)
CUDA by Example. GPU高性能编程CUDA实战 (桑德斯著；聂雪军等译)

下方笔记主要基于阅读 CUDA by Example 书目

1. 为什么需要 CUDA

2. CUDA 环境配置

3. CUDA C 简介

3.1. 本章目标

你将开始编写你的第一段 CUDA C 代码。
你将了解主机（host）代码与设备（device）代码之间的区别。
你将学习如何从主机端调用设备端代码。
你将学习 CUDA 支持的设备上，设备内存的不同使用方式。
你将学习如何查询系统中 CUDA 支持的设备信息。

3.2. 第一个程序

3.2.1. Hello, World!

#include "../common/book.h"

int main(void) {
    printf("Hello, World!\n");
    return 0;
}

CUDA C 与标准 C 在很大程度上没有区别。

CPU 以及系统的内存称为主机（host），而 GPU 及其内存称为设备（device），在 device 上执行的函数通常称为核函数（Kernel）

3.2.2. Kernel 调用

#include <iostream>

__global__ void kernel(void) {
}

int main(void) {
    // 启动 CUDA kernel，使用 1 个 block，1 个线程
    kernel<<<1, 1>>>();
    // 主机端输出
    printf("Hello, World!\n");
    return 0;
}

一个空函数 kernel()，并且带有修饰符 __global__
对这个空函数的调用，并且带有修饰字符 <<<1,1>>>

对于 1.2.1. 中的代码默认是由系统的标准 C 编译器编译的，但这节的例子 CUDA C 为标准 C 增加的 __global__ 修饰符，告诉编译器函数应该编译为在 device 而不是 host 上运行（kernel() 交给编译设备代码的编译器，而 main() 交给主机编译器）。

CUDA C 引入了一种特殊语法 <<<...>>> 将一个函数标记为“设备代码（Device Code）”，并调用设备端（GPU）函数（也称 kernel），这是从主机（CPU）调用设备（GPU）代码的关键机制。

__global__ 修饰的函数是设备代码，运行在 GPU 上；
<<<1, 1>>> 并不是传给函数本身的参数，而是传给 CUDA 运行时的配置参数，用来指定启动多少线程（block 数和线程数）；
小括号 () 内才是传递给 kernel 的实际参数；
这种调用好处是 CUDA C 让设备函数的调用 看起来和普通 C 函数类似；
编译器和运行时系统会处理背后的复杂调用过程。

3.2.3. 传递参数

#include <iostream>
#include "book.h"  // 提供 HANDLE_ERROR 宏，用于错误检查

// 设备端（GPU 上执行）的 kernel 函数
__global__ void add(int a, int b, int *c) {
    *c = a + b;  // 在 GPU 上执行加法，结果写入指针 c 指向的位置
}

int main(void) {
    int c;           // 主机端变量，用于存放结果
    int *dev_c;      // 设备端指针，指向 GPU 上的内存

    // 在设备（GPU）上分配一个整数大小的内存
    HANDLE_ERROR( cudaMalloc( (void**)&dev_c, sizeof(int) ) );

    // 调用 GPU 上的 kernel：add(2, 7, dev_c)
    // <<<1,1>>> 表示使用 1 个 block，1 个线程
    add<<<1, 1>>>(2, 7, dev_c);

    // 从设备内存复制结果到主机变量 c
    HANDLE_ERROR( cudaMemcpy(
        &c,             // 目标：主机端变量地址
        dev_c,          // 源：设备端内存
        sizeof(int),    // 拷贝大小
        cudaMemcpyDeviceToHost  // 拷贝方向：设备 → 主机
    ) );

    printf("2 + 7 = %d\n", c);
    // 释放设备内存
    cudaFree(dev_c);

    return 0;
}

可以像调用 C 函数那样将参数传递给核函数
当设备执行任何有用的操作时，都需要分配内存，如将计算值返回给主机

设备内存分配：cudaMalloc
- 类似于 C 的 malloc，但会在 GPU 上分配内存。
- 调用形式：cudaMalloc(void** ptr, size_t size)
- 注意：返回的设备指针不能在主机代码中解引用（读写）。
设备指针的使用规则
- ✅ 可以在设备代码中使用 cudaMalloc 的指针
- ✅ 可以在主机代码中使用该指针做指针传递
- ❌ 不允许在主机代码中解引用该指针
- ❌ 不允许在设备代码中解引用主机指针
- 主机指针只能在主机访问，设备指针只能在设备访问
内存拷贝：cudaMemcpy
- 用于在主机和设备之间拷贝数据
- 类似 C 的 memcpy，但需要指定方向：
  - cudaMemcpy(dest, src, size, cudaMemcpyHostToDevice); // 主机 → 设备
  - cudaMemcpy(dest, src, size, cudaMemcpyDeviceToHost); // 设备 → 主机
  - cudaMemcpy(dest, src, size, cudaMemcpyDeviceToDevice); // 设备 → 设备
- 主机 ↔ 主机用标准 memcpy 即可
释放设备内存：cudaFree
- 与 cudaMalloc 配套，不能使用标准的 free()。

3.3. 查询设备

使用 cudaGetDeviceCount() 可以获取设备数，cudaGetDeviceProperties() 可以迭代查询 cudaDeviceProp 包含设备相关属性。

3.4. 设备属性的使用

赋值 cudaDeviceProp 中的部分属性，可以使用 cudaChooseDevice() 自动查询某个满足条件的设备并返回其设备 ID，传给 cudaSetDevice() 后，所有的设备操作都将在该设备上执行。

4. CUDA C 并行编程

4.1. 本章目标

了解 CUDA 在实现并行性时采用的一种重要方式
用 CUDA C 编写并行代码

4.2. CUDA 并行编程

只在标准 C 函数前加上 __global__ 修饰符，并通过特殊尖括号语法来调用是很低效的，GPU 的并行运算才是高效的。

4.2.1. 矢量求和运算

CPU 的处理方式是使用循环，对于多线程则是分治，分配不同的索引给不同线程计算。

GPU的处理方式：

#include "../common/book.h"

#define N 10

// 设备端（GPU）上的 kernel 函数，每个线程处理数组中一个元素
__global__ void add(int *a, int *b, int *c) {
    int tid = blockIdx.x;  // 使用 block 索引作为数据索引
    if (tid < N) {
        c[tid] = a[tid] + b[tid];  // 逐元素相加
    }
}

int main(void) {
    int a[N], b[N], c[N];           // 主机端数组
    int *dev_a, *dev_b, *dev_c;     // 设备端指针

    // 在 GPU 上分配内存
    HANDLE_ERROR(cudaMalloc((void**)&dev_a, N * sizeof(int)));
    HANDLE_ERROR(cudaMalloc((void**)&dev_b, N * sizeof(int)));
    HANDLE_ERROR(cudaMalloc((void**)&dev_c, N * sizeof(int)));

    // 填充主机端数组 a 和 b
    for (int i = 0; i < N; i++) {
        a[i] = -i;
        b[i] = i * i;
    }

    // 将 a 和 b 拷贝到 GPU 设备内存中
    HANDLE_ERROR(cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice));
    HANDLE_ERROR(cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice));

    // 启动 kernel，使用 N 个线程，每个线程处理一个元素
    add<<<N, 1>>>(dev_a, dev_b, dev_c);

    // 从 GPU 拷贝结果数组 c 回主机
    HANDLE_ERROR(cudaMemcpy(c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost));

    for (int i = 0; i < N; i++) {
        printf("%d + %d = %d\n", a[i], b[i], c[i]);
    }

    // 释放 GPU 上分配的内存
    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);

    return 0;
}

主机函数 main 中使用了一些通用模式：

调用 cudaMalloc() 在设备上为三个数组分配内存，包含输入值和计算结果
避免内存泄漏，使用完 GPU 内存后通过 cudaFree() 释放
通过 cudaMemcpy() 完成了数据在设备和主机中的相互复制
通过尖括号语法，在主机代码中执行了设备代码

设备内核 add() 中使用了一种通用模式：

采用 C 来编写代码，并在函数名字前添加修饰符 __global__ 。

尖括号的参数与之前有些不同，第一个参数被指定为 N，表示设备在执行核函数时使用的并行线程块的数量。可以认为运行时将创建核函数的 N 个副本，并以并行方式运行，每个并行执行环境都称为一个线程块（Block）。这个 并行线程块的集合称为一个线程格（Grid）。

核函数通过内置变量（CUDA 运行时已经预先定义的变量） blockIdx 表示当前执行设备代码的线程块的索引。blockIdx.x 的原因：CUDA 支持二维的线程块数组，为了便利一些二维空间的计算。

需要注意硬件限制：启动线程块数组时，数组的每一维的最大数量都不能超过 65 535。

4.2.2. Julia 集示例

基于 GPU 的 Julia 集：

#include "../common/book.h"
#include "../common/cpu_bitmap.h"

#define DIM 1000  // 图像尺寸
// 复数类（设备端可用）
struct cuComplex {
    float r;
    float i;

    cuComplex(float a, float b) : r(a), i(b) {}

    __device__ float magnitude2(void) {
        return r * r + i * i;
    }

    __device__ cuComplex operator*(const cuComplex& a) {
        return cuComplex(r * a.r - i * a.i, i * a.r + r * a.i);
    }

    __device__ cuComplex operator+(const cuComplex& a) {
        return cuComplex(r + a.r, i + a.i);
    }
};

// 判断点 (x, y) 是否属于 Julia 集
__device__ int julia(int x, int y) {
    const float scale = 1.5;
    float jx = scale * (float)(DIM / 2 - x) / (DIM / 2);
    float jy = scale * (float)(DIM / 2 - y) / (DIM / 2);

    cuComplex c(-0.8, 0.156);
    cuComplex a(jx, jy);

    int i = 0;
    for (i = 0; i < 200; i++) {
        a = a * a + c;
        if (a.magnitude2() > 1000) return 0;  // 发散了
    }
    return 1;  // 属于 Julia 集
}

// CUDA kernel，每个线程负责一个像素点
__global__ void kernel(unsigned char *ptr) {
    int x = blockIdx.x;
    int y = blockIdx.y;
    int offset = x + y * gridDim.x;

    // 判断该点是否属于 Julia 集
    int juliaValue = julia(x, y);

    // 设置图像颜色（红色通道根据结果决定，其他为固定值）
    ptr[offset * 4 + 0] = 255 * juliaValue;  // R
    ptr[offset * 4 + 1] = 0;                 // G
    ptr[offset * 4 + 2] = 0;                 // B
    ptr[offset * 4 + 3] = 255;               // A
}

int main(void) {
    CPUBitmap bitmap(DIM, DIM);
    unsigned char *dev_bitmap;

    // 在 GPU 上分配图像内存
    HANDLE_ERROR(cudaMalloc((void**)&dev_bitmap, bitmap.image_size()));

    // 设置每个线程负责一个像素（DIM x DIM 的网格）
    dim3 grid(DIM, DIM);
    kernel<<<grid, 1>>>(dev_bitmap);

    // 将结果从 GPU 拷贝回主机
    HANDLE_ERROR(cudaMemcpy(bitmap.get_ptr(), dev_bitmap,
                            bitmap.image_size(),
                            cudaMemcpyDeviceToHost));

    // 显示图像
    bitmap.display_and_exit();

    // 释放 GPU 内存
    HANDLE_ERROR(cudaFree(dev_bitmap));
}

dim3 并不是标准 C 定义的类型，而是 CUDA 头文件中定义的辅助类型来封装多维数组，CUDA 运行时希望得到三维的 dim3 值（二维的话最后一维大小为 1）。

修饰符 __device__，表示代码将在 GPU 而不是主机上运行，且只能从其他 __device__ 函数或者从 __global__ 函数中调用它们。

5. 线程协作

5.1. 本章目标

了解 CUDA C 中的线程
了解不同线程之间的通信机制
了解并行执行线程的同步机制

5.2. 并行线程块的分解

add<<<N,1>>>( dev_a, dev_b, dev_c );

尖括号中，第一个参数表示设备在执行核函数时使用的并行线程块的数量，第二个参数表示 CUDA 运行时在每个线程块中创建的线程数量，上面代码表示在每个线程块中只启动一个线程，总共启动的线程数量可以按以下公式计算：

N 个线程块 × 1 个线程/线程块 = N 个并行线程

事实上，可以启动 N/2 个线程块，每个线程块包含 2 个线程等等分配方式。

5.2.1. 矢量求和：重新回顾

线程块中的并行线程能够完成并行线程块无法完成的工作。

5.2.1.1. 使用线程实现GPU上的矢量求和

若把 4.2.1. 中的示例，将使用并行线程块修改为使用并行线程需要做出两处修改

// add<<<N,1>>>( dev_a, dev_b, dev_c );
add<<<1,N>>>( dev_a, dev_b, dev_c );

// int tid = blockIdx.x;
int tid = threadIdx.x;

5.2.1.2. 在GPU上对更长的矢量求和

硬件将线程块的数量限制为不超过65535，同样，对于启动核函数时每个线程块中的线程数量，硬件也限制在族弟啊不能超过设备属性中 maxTreadsPerBlock 域的值。

将线程与线程块结合起来可以实现较长的矢量求和，需要改动两个地方：核函数中的索引计算方法和核函数的调用方式。

核函数中的索引计算方法：

此时需要 多个线程块且每个线程块包含多个线程，计算索引方法类似将二维索引空间转换为线性空间的标准算法。

int tid = threadIdx.x + blockIdx.x * blockDim.x;

blockDim 内置变量，对于所有线程块来说，该变量是常数，保存线程块中每一维的线程数量，回顾第 4 章，gridDim 中保存了类似的值，即 线程格中每一维的线程块数量。gridDim 是二维的，而 blockDim 实际上是三维的，即 CUDA 运行时允许启动一个二维线程格，线程格中每个线程块都是一个三维的线程数组，这是一种高维数组。

核函数的调用方式：

另一处修改是核函数调用本身，虽然仍需启动 N 个并行线程，但希望在多个线程块中启动，一种解决方案是，将线程块的大小设定为某个固定数值，计算 (N+127)/128（即大于或等于 N 的 128 的最小倍数）。

add<<< (N+127)/128, 128 >>>( dev_a, dev_b, dev_c );

上述代码为确保启动足够多的线程，当 N 不是 128 的整数倍时，将启动过多的线程。因此，在核函数中，访问输入和输出数组之前，必须检查线程的偏移是否位于 0 到 N 之间。

if (tid < N)
    c[tid] = a[tid] + b[tid];

因此当索引越过数组边界时，核函数将自动停止执行计算，不会对越过数组边界的内存进行读取或写入。

5.2.1.3. 在GPU上对任意长度的矢量求和

#include "../common/book.h"

#define N (33 * 1024)  // 33 × 1024 = 33792 元素

// CUDA kernel：多线程并行执行加法
__global__ void add(int *a, int *b, int *c) {
    // 全局线程索引
    int tid = threadIdx.x + blockIdx.x * blockDim.x;

    // **每个线程处理多个元素，步长为总线程数**
    while (tid < N) {
        c[tid] = a[tid] + b[tid];
        tid += blockDim.x * gridDim.x;
    }
}

int main(void) {
    int a[N], b[N], c[N];              // 主机端数组
    int *dev_a, *dev_b, *dev_c;        // 设备端指针

    // 在 GPU 上分配内存
    HANDLE_ERROR(cudaMalloc((void**)&dev_a, N * sizeof(int)));
    HANDLE_ERROR(cudaMalloc((void**)&dev_b, N * sizeof(int)));
    HANDLE_ERROR(cudaMalloc((void**)&dev_c, N * sizeof(int)));

    // 初始化主机端数据
    for (int i = 0; i < N; i++) {
        a[i] = i;
        b[i] = i * i;
    }

    // 拷贝 a 和 b 到 GPU
    HANDLE_ERROR(cudaMemcpy(dev_a, a, N * sizeof(int), cudaMemcpyHostToDevice));
    HANDLE_ERROR(cudaMemcpy(dev_b, b, N * sizeof(int), cudaMemcpyHostToDevice));

    // **启动 kernel：128 blocks，每 block 128 线程（共 16384 个线程）**
    add<<<128, 128>>>(dev_a, dev_b, dev_c);

    // 从 GPU 拷贝结果数组 c 回主机
    HANDLE_ERROR(cudaMemcpy(c, dev_c, N * sizeof(int), cudaMemcpyDeviceToHost));

    // 验证 GPU 确实完成了我们要求的工作
    bool success = true;
    for (int i = 0; i < N; i++) {
        if ((a[i] + b[i]) != c[i]) {
            printf("Error: %d + %d != %d\n", a[i], b[i], c[i]);
            success = false;
        }
    }

    if (success) {
        printf("We did it!\n");
    }

    // 释放设备内存
    cudaFree(dev_a);
    cudaFree(dev_b);
    cudaFree(dev_c);

    return 0;
}

核函数中做出的修改：

// **每个线程处理多个元素，步长为总线程数**
while (tid < N) {
    c[tid] = a[tid] + b[tid];
    tid += blockDim.x * gridDim.x;
}

核心思想：在 GPU 实现中，并行线程的数量可以看成是处理器的数量，认为每个线程在逻辑上都可以并行执行并可以被硬件调度。

当线程计算完当前索引上的任务后，需要对索引进行递增，其补偿为线程格中正在运行的线程数量，该数值等于每个线程块中的线程数量乘以线程格中线程块的数量，即 blockDim.x * gridDim.x 。

采取此方式的原因是因为先前当 (N+127)/128 大于 65 535 时，核函数调用 add<<< (N+127)/128, 128 >>>( dev_a, dev_b, dev_c ); 会失败，为确保不会启动过多的线程块，可以将线程块的数量固定为某个较小值。

add<<< 128, 128 >>>( dev_a, dev_b, dev_c );

5.2.2. 在 GPU 上使用线程实现波纹效果

// 数据结构：封装主机与设备图像数据
struct DataBlock {
    unsigned char *dev_bitmap;   // 指向设备（GPU）上的图像数据
    CPUAnimBitmap *bitmap;       // 指向主机上的 bitmap 对象
};

// 清理 GPU 分配的内存
void cleanup(DataBlock *d) {
    cudaFree(d->dev_bitmap);
}

// 主程序入口
int main(void) {
    DataBlock data;

    // 创建 CPU 上的 bitmap 并将 data 传入作为用户数据
    CPUAnimBitmap bitmap(DIM, DIM, &data);
    data.bitmap = &bitmap;

    // 在 GPU 上分配 bitmap 所需内存
    HANDLE_ERROR(cudaMalloc(
        (void**)&data.dev_bitmap,
        bitmap.image_size()
    ));

    // 启动动画与退出回调函数
    bitmap.anim_and_exit(
        (void (*)(void*, int))generate_frame,  // 每帧生成函数
        (void (*)(void*))cleanup               // 清理函数
    );
}

void generate_frame(DataBlock *d, int ticks) {
    dim3 blocks(DIM / 16, DIM / 16);
    dim3 threads(16, 16);

    // 启动 GPU 核函数进行图像生成
    kernel<<<blocks, threads>>>(d->dev_bitmap, ticks);

    // 将生成的图像数据从设备复制回主机
    HANDLE_ERROR(cudaMemcpy(
        d->bitmap->get_ptr(),
        d->dev_bitmap,
        d->bitmap->image_size(),
        cudaMemcpyDeviceToHost
    ));
}

对于 1920 × 1080 的动画，这种访华会创建超过 2 百万格线程，GPU 能实现 CPU 无法实现的，支持巨量的线程。

核函数 kernel 包含两个参数，第一个参数指针指向保存输出像素值的设备内存，内存是在 main() 中分配的，是一个全局变量，只是对于主机代码“全局”，需要作为参数传入确保设备代码能访问。其次，核函数需要知道当前动画时间以生成正确帧，第二个参数就是时间 ticks。

__global__ void kernel(unsigned char *ptr, int ticks) {
    // 将线程索引映射到图像中的像素位置
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    // 计算当前像素在图像数据中的偏移量
    int offset = x + y * blockDim.x * gridDim.x;

    // 计算当前像素到图像中心的距离
    float fx = x - DIM / 2;
    float fy = y - DIM / 2;
    float d = sqrtf(fx * fx + fy * fy);

    // 根据距离和时间 ticks 计算灰度值，形成动画波动效果
    unsigned char grey = (unsigned char)(128.0f + 127.0f *
        cos(d / 10.0f - ticks / 7.0f) /
        (d / 10.0f + 1.0f));

    // 设置 RGBA 四通道值（灰度图）
    ptr[offset * 4 + 0] = grey;  // R
    ptr[offset * 4 + 1] = grey;  // G
    ptr[offset * 4 + 2] = grey;  // B
    ptr[offset * 4 + 3] = 255;   // A（不透明）
}

kernel 的前三行代码是最重要的代码，每个线程都计算的到了其在图形中的唯一索引。例如，对于线程格中 (12, 8) 处的线程块的 (3, 5) 处的线程：

3 threads + 12 blocks * 16 threads/block = 195 threads to the left of it
5 threads + 8 blocks * 16 threads/block = 128 threads above it

5.3. 共享内存和同步

到目前为止，将线程块分解为线程的目的只是为了解决线程块数量的硬件限制，是较为勉强的动机，完全可以由 CUDA 运行时在幕后自动完成。

线程的动机不只是硬件限制
- 最初将 block 拆成多个 thread 只是为了解决硬件对 block 数量的限制，但这并不是唯一理由。
对 C 语言的另一个扩展，引入共享内存（__share__）
- __shared__ 关键字声明的变量存在于 GPU 上的物理共享内存中。
- 每个 block 拥有自己的共享内存副本，block 内的所有线程共享访问，但 不同 block 间互不可见。
共享内存的优势
- 位于 GPU 上，访问速度远快于全局内存（off-chip DRAM）。
- 适合作为 block 内线程协作的缓存或 中间计算缓冲区（scratchpad）。
线程间通信与同步
- 多线程共享内存后，可以通信与协作。
- 但需使用同步机制（如 __syncthreads()）来防止 竞态条件（race conditions）。
- 没有同步会导致程序执行结果取决于硬件的非确定性行为。

5.3.1. 点积运算

#include "../common/book.h"

#define imin(a, b) (a < b ? a : b)

const int N = 33 * 1024;             // 向量长度
const int threadsPerBlock = 256;     // 每个 block 中的线程数

__global__ void dot(float *a, float *b, float *c) {
    // 声明每个 block 使用的共享内存，供线程协作使用
    __shared__ float cache[threadsPerBlock];

    // 全局线程 ID（对应数据下标）
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    int cacheIndex = threadIdx.x;

    // 局部累加器
    float temp = 0;

    // 每个线程负责处理多个元素（跨 grid）
    while (tid < N) {
        temp += a[tid] * b[tid];
        tid += blockDim.x * gridDim.x;
    }

    // 将本线程计算结果存入共享内存
    cache[cacheIndex] = temp;

    // （后续添加同步和归约逻辑）
}

缓冲区 cache 将保存每个线程计算的加和值，声明一个 驻留在共享内存 中的变量相当于在标准 C 中声明一个 static 或 volatile 类型的变量。

分配全局内存是，为每个执行核函数的线程都分配了足够的内存；但对共享变量，编译器将为 每个线程块 生成该共享变量的一个 私有副本，因此只需根据线程块中线程的数量来分配内存。

当算法执行完上述过程后，需要将 cache 中所有值相加，需要一个线程来读取保存在 cache 中的值。由于该操作较为危险，需要某种方法来确保所有对共享数组 cache[] 的写入操作在读取 cache 之前完成。

// 对线程块中的线程进行同步
__syncthreads();

该函数调用将确保线程块中的每个线程都执行完 __syncthreads() 前的语句后，才会执行下一条语句。

相加过程可以抽象为更一般形式：对一个输入数组执行某种计算，然后产生一个更小的结果数组，也成为 归约（Reduction）。

归约运算的朴素方法是由一个线程在共享内存上进行迭代并计算出总和值，计算时间与数组长度成正比。但由于此处有数百个线程可用，因此可以以并行方式来执行归约运算，时间将与数组长度的对数成正比。

// 对共享内存中的数组进行归约操作
int i = blockDim.x / 2;
while (i != 0) {
    if (cacheIndex < i)
        cache[cacheIndex] += cache[cacheIndex + i];

    // 同步线程，确保所有更新完成
    __syncthreads();

    i /= 2;
}

// 若全部归约结束，拷贝全局内存
if (cacheIndex == 0)
    c[blockIdx.x] = cache[0];

每个线程块都需要一个线程来将归约值写入全局内存。

对于最后所有线程块归约得到的值的归约，对于 GPU 这种大规模并行机器，通常会浪费计算资源，因为数据集往往非常小。因此，将执行控制返回给主机，并且由 CPU 完成最后的归约。

const int blocksPerGrid = imin( 32, (N+threadsPerBlock-1) / threadsPerBlock );

dot<<<blocksPerGrid,threadsPerBlock>>>( dev_a, dev_b, dev_partial_c );

经验值选择 block 数量：
- 通常选择 32 个 blocks，足够让 GPU 忙碌，同时结果数量也便于 CPU 处理。
处理小数据量的技巧：
- 如果数据量 N 比线程总数还小，则只需启动 刚好够用的线程。
计算最小足够 block 数的通用公式：
- (N + threadsPerBlock - 1) / threadsPerBlock：向上取整
- min(32, ...)：不超过经验最大 block 数

5.3.2. （不正确的）点积运算优化

int i = blockDim.x / 2;
while (i != 0) {
    if (cacheIndex < i) {
        cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads(); // 若移入 if 判断中
    }

    i /= 2;
}

GPU 将停止响应，因为 cacheIndex 非小于 i 的线程无法执行 __syncthreads() 指令导致无法同步。当某些线程需要执行一条指令，而其他线程不需要执行时，这种情况称为 线程发散（Thread Divergence），会使得某些线程处于空闲状态。

使用 __syncthreads() 指令要注意不能出现在发散分支中。

5.3.3. 基于共享内存的位图

6. 常量内存与事件

6.1. 本章目标

了解如何在 CUDA C 中使用常量内存
了解常量内存的性能特性
学习如何使用 CUDA 事件来测量应用程序的性能

6.2. 常量内存

由于 GPU 上包含非常多的数学逻辑单元（ALU），因此性能瓶颈通常并不在芯片的数学计算吞吐量，而在于芯片的内存带宽，因此有时输入数据的速率甚至无法维持如此高的计算速率，需要手段来减少计算问题时的内存通信量。

常量内存 用于保存 在核函数执行期间不会发生变化的数据，NVIDIA 对常量内存采用了不同于标准全局内存的处理方式。在某些情况中，用常量内存来替换全局内存能有效地减少内存带宽。

6.2.1. 光线跟踪简介

光线追踪（ray tracing）的基本原理：与传统 GPU 的光栅化（rasterization）方法不同，光线追踪通过从每个像素“发射”射线进入三维场景，判断其与物体的交点，确定像素颜色。它可以自然扩展到处理反射（反光物体）、折射（半透明物体）等复杂效果，这些效果将生成二次射线（Secondary Ray）和三次射线（Tertiary Ray），从而生成更真实的图像。

6.2.2. 在 GPU 上实现光线跟踪

下方介绍一个基于 CUDA C 实现的简化光线追踪器：为了展示如何使用 常量内存（constant memory），该光线追踪器仅支持球体场景，摄像机固定在 z 轴朝向原点，不计算光照，仅判断每个像素发出的射线是否与球体相交，并选取最靠近相机的球体赋予颜色，主要用于演示射线与物体的遮挡关系。

#define INF 2e10f

struct Sphere {
    float r, g, b;        // 颜色分量（红、绿、蓝）
    float radius;         // 半径
    float x, y, z;        // 球心坐标

    // 在设备上计算光线与球的交点
    __device__ float hit(float ox, float oy, float *n) {
        float dx = ox - x;
        float dy = oy - y;

        // 判断射线是否命中该球（投影距离是否小于半径）
        if (dx * dx + dy * dy < radius * radius) {
            float dz = sqrtf(radius * radius - dx * dx - dy * dy);
            *n = dz / radius; // 法线信息归一化
            return dz + z;    // 返回击中点的深度（z 方向）
        }

        return -INF; // 没有命中，返回一个极小值
    }
};

在这段代码中并没有计算完整的法线向量（即 vec = (dx, dy, dz)）；而是只使用了 z 方向的信息来快速估算一个“垂直程度”；这是一个近似值，适用于简单着色场景而非真实光照模型。

main() 函数与之前大致相同

#include "cuda.h"
#include "../common/book.h"
#include "../common/cpu_bitmap.h"

#define rnd(x) (x * rand() / RAND_MAX)
#define SPHERES 20

Sphere *s;

int main(void) {
    // 创建并记录 CUDA 事件用于测量执行时间
    cudaEvent_t start, stop;
    HANDLE_ERROR(cudaEventCreate(&start));
    HANDLE_ERROR(cudaEventCreate(&stop));
    HANDLE_ERROR(cudaEventRecord(start, 0));

    // 创建位图对象
    CPUBitmap bitmap(DIM, DIM);
    unsigned char *dev_bitmap;

    // 为 GPU 输出图像分配显存
    HANDLE_ERROR(cudaMalloc((void**)&dev_bitmap, bitmap.image_size()));

    // 为球体数据分配显存
    HANDLE_ERROR(cudaMalloc((void**)&s, sizeof(Sphere) * SPHERES));

    // 分配临时内存，对其初始化，并复制到 GPU 上的内存，然后释放临时内存
    Sphere *temp_s = (Sphere*)malloc( sizeof(Sphere) * SPHERES );
    for (int i=0; i<SPHERES; i++) {
        temp_s[i].r = rnd( 1.0f );
        temp_s[i].g = rnd( 1.0f );
        temp_s[i].b = rnd( 1.0f );
        temp_s[i].x = rnd( 1000.0f ) - 500;
        temp_s[i].y = rnd( 1000.0f ) - 500;
        temp_s[i].z = rnd( 1000.0f ) - 500;
        temp_s[i].radius = rnd( 100.0f ) + 20;
    }

    HANDLE_ERROR( cudaMemcpy( s, temp_s, sizeof(Sphere) * SPHERES, cudaMemcpyHostToDevice ) );
    free( temp_s );

    // 从球面数据中生成一张位图
    dim3 grids(DIM/16,DIM/16);
    dim3 threads(16,16);
    kernel<<<grids,threads>>>( dev_bitmap );

    // 将位图从 GPU 复制回到 CPU 以显示
    HANDLE_ERROR( cudaMemcpy( bitmap.get_ptr(), dev_bitmap,
                              bitmap.image_size(),
                              cudaMemcpyDeviceToHost ) );
    bitmap.display_and_exit();

    // 释放内存
    cudaFree( dev_bitmap );
    cudaFree( s );
}

在临时内存上随机生成球面数组，再复制到 GPU，然后释放临时缓冲区，且为输出数据分配好空间，可以启动核函数。

__global__ void kernel(unsigned char *ptr) {
    // 映射线程到图像像素位置
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;

    // 将像素位置转换为以中心为原点的坐标系
    float ox = x - DIM / 2;
    float oy = y - DIM / 2;

    // 初始化颜色和最大深度值
    float r = 0, g = 0, b = 0;
    float maxz = -INF;

    // 遍历所有球体，查找当前像素射线与球体的最近相交
    for (int i = 0; i < SPHERES; i++) {
        float n;
        float t = s[i].hit(ox, oy, &n);
        if (t > maxz) {
            // 更新最大深度，并根据法线近似值 n 决定颜色强度
            float fscale = n;
            r = s[i].r * fscale;
            g = s[i].g * fscale;
            b = s[i].b * fscale;
            maxz = t;
        }
    }

    // 写入图像数据（RGBA）
    ptr[offset * 4 + 0] = (int)(r * 255);
    ptr[offset * 4 + 1] = (int)(g * 255);
    ptr[offset * 4 + 2] = (int)(b * 255);
    ptr[offset * 4 + 3] = 255;  // alpha 通道不透明
}

6.2.3. 通过常量内存来实现光线跟踪

常量内存 不能修改，对于上述示例只有一个输入数据，即球面数组，因此应当保存其至常量内存中。声明方法即在变量前面加上 __constant__ 修饰符：

// Sphere *s;
__constant__ Sphere s[SPHERES];

最初的示例中声明的指针，是通过 cudaMalloc() 来为指针分配 GPU 内存，当修改为常量内存时，同样要将声明修改为 在常量内存中静态地分配空间。但不需要再调用 cudaMalloc() 或 cudaFree() ，而是在编译时为该数组提交一个固定的大小。

    // HANDLE_ERROR( cudaMemcpy( s, temp_s, sizeof(Sphere) * SPHERES, cudaMemcpyHostToDevice ) );
    HANDLE_ERROR( cudaMemcpyToSymbol( s, temp_s, sizeof(Sphere) * SPHERES ) );

    // cudaFree( s );

当从主机内存复制到 GPU 上的常量内存时，需要使用这个特殊版本的 cudaMemcpyToSymbol() ，其与参数为 cudaMemcpyHostToDevice 的 cudaMemcpy() 之间的唯一差异在与，前者会复制到常量内存，而后者会复制到全局内存。

6.2.4. 常量内存带来的性能提升

将内存声明为 __constant__ 会限制其只能进行只读操作，但这种限制带来了性能上的优势。与从全局内存读取相比，从常量内存读取可以节省内存带宽，原因有以下两点：

广播机制：
如果一个“半warp”（即16个线程）中的所有线程请求常量内存中的同一地址，硬件只需发起一次读取请求，并将数据广播给所有线程。这意味着最多可以节省15次读取操作。
缓存机制：
常量内存是缓存的，因此对同一地址的连续读取不会产生额外的内存流量。

什么是“warp”和“半warp”？

在CUDA架构中，“warp”指的是一组32个线程，这些线程在执行程序时是同步的（即“锁步执行”，Lockstep），在程序的每一行，线程束中的每个线程都将在不同的数据上执行相同的指令。一个“半warp”（Half-Warp）是一个包含16个线程的子集。

常量内存的性能优势

广播机制：
当半warp中的每个线程从常量内存中读取相同的地址时，GPU只需发起一次读取请求，然后将数据广播给所有线程。相比从全局内存读取，这种方式可以将内存流量减少到1/16（约6%）。
缓存机制：
如果某个地址的数据被缓存，后续对该地址的读取几乎不会产生额外的内存流量。

例如，在光线追踪（ray tracing）应用中，所有线程需要读取同一个球体的数据来计算光线与球体的相交情况。将这些数据存储在常量内存中后，硬件只需发起一次读取请求，然后通过缓存和广播机制满足其他线程的读取需求，大大减少了内存流量。

常量内存的潜在劣势

尽管常量内存在某些情况下可以显著提高性能，但其广播机制也可能导致性能瓶颈。当半warp中的每个线程需要读取不同的地址时，请求会被串行化，即16个读取请求需要逐一完成。这会显著降低性能，甚至可能比从全局内存读取还慢。在全局内存中，这些读取请求可以同时发出。

6.3. 使用事件来测量性能

为了评估修改后的光线追踪程序是否因使用常量内存而提升了性能，可以通过比较两种实现的运行时间来判断。尽管可以使用CPU或操作系统的计时器，但这些方法会受到各种因素的影响（如操作系统线程调度、高精度CPU计时器的可用性等），并且可能包含主机端与设备端的异步计算时间。为此，最佳选择是使用 CUDA 事件（CUDA Event）API 来测量GPU执行任务的时间。

6.3.1. 什么是CUDA事件？

CUDA 事件本质上是 GPU 记录的时间戳，用于标记用户指定的时间点。由于时间戳是由 GPU 自身记录的，它避免了使用 CPU 计时器时可能遇到的问题，比如主机和设备的异步性。

6.3.2. 使用CUDA事件进行计时的基本步骤

cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
// 执行GPU上的任务
cudaEventRecord(stop, 0);

1. 创建和记录事件

在需要开始计时的代码段开始处，创建一个事件并记录当前时间：

cudaEventCreate(&start)：创建一个事件。
cudaEventRecord(start, 0)：记录当前时间。参数 0 暂时不需要关注（与CUDA流有关，后续会讨论）。

2. 记录结束时间

在代码段结束处，创建另一个事件并记录结束时间.

6.3.3. 异步性问题及解决方法

CUDA 中的许多操作（如内核启动）是 异步的，即 GPU 开始执行任务时，CPU 会继续执行程序的下一行代码，而不会等待 GPU 完成。这种异步性虽然有助于性能提升（可以同时在 GPU 和 CPU 上进行计算），但会使计时变得复杂。

cudaEventRecord() 实际上是将记录时间的指令放入GPU的任务队列中，只有在GPU完成了队列中之前的所有工作后，事件才会真正被记录。
如果不处理这种异步性，stop 事件的时间可能不准确。

解决方法：

使用cudaEventSynchronize()函数确保GPU完成所有工作后再读取事件时间戳：

cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start, 0);
// 执行GPU上的任务
cudaEventRecord(stop, 0);
cudaEventSynchronize(stop);

cudaEventSynchronize(stop)：让CPU等待，直到GPU完成 stop 事件之前的所有工作。只有在此函数返回后，才能安全地读取stop事件的时间戳。

注意：由于 CUDA 事件是直接在 GPU 上实现的，因此它们不适用于对同时包含设备代码和主机代码的混合代码计时。也就是说，如果试图通过 CUDA 事件对核函数和设备内存复制之外的代码进行计时，将得到不可靠的结果。

6.3.4. 测量光线跟踪器的性能

HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime, start, stop ) );
printf( "Time to generate: %3.1f ms\n", elapsedTime );
HANDLE_ERROR( cudaEventDestroy( start ) );
HANDLE_ERROR( cudaEventDestroy( stop ) );

调用了两个额外的函数，分别为 cudaEventElapsedTime() 和 cudaEventDestroy()
cudaEventElapsedTime() ：是一个工具函数，用来计算两个事件之间经历的时间。第一个参数为某个浮点变量的地址，在这个参数中将包含两次事件之间经历的时间，单位为毫秒。
cudaEventDestroy() ：当使用完事件后，需要调用该函数来销毁它们，相当于对 malloc() 分配的内存调用 free() ，因此每个 cudaEventCreate() 都对应一个 cudaEventDestroy() 同样是非常重要的。

7. 纹理内存

7.1. 本章目标

了解纹理内存的性能特性
了解如何在 CUDA C 中使用一维和二维纹理内存

7.2. 纹理内存简介

纹理内存（Texture Memory）和常量内存一样，是另一种类型的只读内存，同样缓存在芯片上，能够减少对内存的请求并提供更高效的内存带宽。纹理缓存是专门为那些在内存访问模式中存在大量空间局部性（Spatial Locality）的图形应用程序而设计的。比如说二维矩阵（若要访问局部 2×2 的小块，地址非连续），在这种情况下，纹理内存相对全局内存能获得性能提升。

7.3. 热传导模拟

7.3.1. 简单的传热模型

为展示如何有效使用 CUDA 的纹理内存，使用了一个二维热传导模拟：将一个房间划分为网格，并随机设置一些恒温加热器。每个时间步中，非加热器单元的温度根据上下左右邻居的温度差进行更新，模拟热量从高温单元向低温单元扩散的过程。

等式7.1

\[T_{\text{NEW}} = T_{\text{OLD}} + \sum\limits_{\text{NEIGHBORS}}k(T_{\text{NEIGHBORS}} – T_{\text{OLD}})\]

等式7.2

\[T_{\text{NEW}} = T_{\text{OLD}} +k(T_{\text{TOP}} + T_{\text{BOTTOM}} + T_{\text{LEFT}} + T_{\text{RIGHT}} – 4T_{\text{OLD}})\]

常量 k 表示模拟过程中热量的流动速率。k 值越大，表示系统会更快地达到温度，相反，则温度梯度将存在更长的时间。

该模拟的更新方程简单易计算，强调教学目的而非物理准确性，目的是为了说明如何在 GPU 上用纹理内存高效访问邻居单元数据，加速计算。

7.3.2. 温度更新的计算

详细代码见原书对应章节，热传导模拟的每一帧包含以下三个步骤：

复制加热器温度（copy_const_kernel）
将加热器非零的固定温度从常量网格 cptr 复制到输入温度网格 iptr，确保加热器单元不会随时间变化。
更新温度（blend_kernel）
每个线程负责一个网格单元，读取其本身和上下左右邻居的温度，依据公式进行温度更新。
交换缓冲区
更新后的 outSrc 成为下一时间步的输入 inSrc。

__global__ void blend_kernel(float *outSrc, const float *inSrc) {
    // 将线程和块索引映射为图像坐标
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;

    // 计算上下左右邻居的索引
    int left   = offset - 1;
    int right  = offset + 1;
    int top    = offset - DIM;
    int bottom = offset + DIM;

    // 边界条件处理
    if (x == 0)        left++;
    if (x == DIM - 1)  right--;
    if (y == 0)        top += DIM;
    if (y == DIM - 1)  bottom -= DIM;

    // 更新当前像素的值
    outSrc[offset] = inSrc[offset] + SPEED * (inSrc[top] + inSrc[bottom] + inSrc[left] + inSrc[right] - inSrc[offset] * 4);
}

7.3.3. 模拟过程动态演示

详细代码见原书对应章节，通过一个 anim_gpu() 函数，每帧执行 90 次上述温度更新步骤（1~3）。

void anim_gpu( DataBlock *d, int ticks ) {
    HANDLE_ERROR( cudaEventRecord( d->start, 0 ) );
    dim3 blocks(DIM/16,DIM/16);
    dim3 threads(16,16);
    CPUAnimBitmap *bitmap = d->bitmap;

    // 进行 90 个时间步的模拟更新
    for (int i=0; i<90; i++) {
        copy_const_kernel<<<blocks,threads>>>( d->dev_inSrc, d->dev_constSrc );
        blend_kernel<<<blocks,threads>>>( d->dev_outSrc, d->dev_inSrc );
        swap( d->dev_inSrc, d->dev_outSrc );
    }
    // 将浮点温度数据转换为颜色并写入输出图像
    float_to_color<<<blocks,threads>>>( d->output_bitmap, d->dev_inSrc );
    // 将结果从 GPU 拷贝到 CPU 显示
    HANDLE_ERROR( cudaMemcpy( bitmap->get_ptr(), d->output_bitmap, bitmap->image_size(), cudaMemcpyDeviceToHost ) );
    // 记录 GPU 结束事件并计算帧耗时
    HANDLE_ERROR( cudaEventRecord( d->stop, 0 ) );
    HANDLE_ERROR( cudaEventSynchronize( d->stop ) );
    float elapsedTime;
    HANDLE_ERROR( cudaEventElapsedTime( &elapsedTime, d->start, d->stop ) );
    // 累积时间统计并打印平均帧时间
    d->totalTime += elapsedTime;
    ++d->frames;
    printf( "Average Time per frame: %3.1f ms\n", d->totalTime/d->frames );
}

7.3.4. 使用纹理内存

首先，需要将输入的数据声明为 texture 类型的引用。

// 这些变量将位于 GPU 上
texture<float> texConstSrc;
texture<float> texIn;
texture<float> texOut;

在为这三个缓冲区分配了 GPU 内存后，需要通过 cudaBindTexture() 将这些变量绑定到内存缓冲区，相当于告诉 CUDA 运行时两件事：

希望将指定的缓冲区作为纹理来使用
希望将纹理引用作为纹理的“名字”

// 分配 GPU 内存
HANDLE_ERROR(cudaMalloc((void**)&data.dev_inSrc,    imageSize));
HANDLE_ERROR(cudaMalloc((void**)&data.dev_outSrc,   imageSize));
HANDLE_ERROR(cudaMalloc((void**)&data.dev_constSrc, imageSize));

// 绑定纹理内存
HANDLE_ERROR(cudaBindTexture(NULL, texConstSrc,  data.dev_constSrc, imageSize));
HANDLE_ERROR(cudaBindTexture(NULL, texIn,        data.dev_inSrc,    imageSize));
HANDLE_ERROR(cudaBindTexture(NULL, texOut,       data.dev_outSrc,   imageSize));

启动核函数后，当读取核函数的纹理时，需要通过特殊的函数来告诉 GPU 将读取请求转发到纹理内存而非标准的全局内存。因此，当读取内存时不再使用方括号从缓冲区中读取，而是在 blend_kernel() 中 // 更改当前像素的值 的部分改为使用 tex1Dfetch() 。

虽然 tex1Dfetch() 像函数，但其实是 编译器内置函数（Intrinsic）。由于纹理引用必须声明为文件作用域内的全局变量，因此不再将输入输出缓冲区作为参数传递给 blend_kernel() ，因为编译器需要在编译时知道 tex1Dfetch() 应该对应哪些纹理采样，而是将一个布尔标志 dstOut 传递给函数，该 symbol 会告诉我们使用哪个缓冲区作为输入输出。

__global__ void blend_kernel(float *dst, bool dstOut) {
    // ... 前面与 7.3.2.中代码相同

    // 更新当前像素的值
    // outSrc[offset] = inSrc[offset] + SPEED * (inSrc[top] + inSrc[bottom] + inSrc[left] + inSrc[right] - inSrc[offset] * 4);
    float t, l, c, r, b;
    if (dstOut) {
        t = tex1Dfetch(texIn,top);
        l = tex1Dfetch(texIn,left);
        c = tex1Dfetch(texIn,offset);
        r = tex1Dfetch(texIn,right);
        b = tex1Dfetch(texIn,bottom);
    } else {
        t = tex1Dfetch(texOut,top);
        l = tex1Dfetch(texOut,left);
        c = tex1Dfetch(texOut,offset);
        r = tex1Dfetch(texOut,right);
        b = tex1Dfetch(texOut,bottom);
    }
    dst[offset] = c + SPEED * (t + b + r + l - 4 * c);
}

由于 copy_const_kernel 将读取包含热源位置和温度的缓冲区，因此同样要修改为从纹理内存中读取。

// if (cptr[offset] != 0) iptr[offset] = cptr[offset];
float c = tex1Dfetch(texConstSrc,offset);
if (c != 0)
    iptr[offset] = c;

由于 blend_kernel() 的函数圆形被修改为接受一个标志，并且该标志表示在输入与输出缓冲区之间的切换，因此需要修改 anim_gpu() 函数，从交换缓冲区改为在每组调用之后通过设置 dstOut = !dstOut 来进行切换：仅修改 // 进行 90 个时间步的模拟更新 部分

void anim_gpu( DataBlock *d, int ticks ) {
    //  ... 前面与 7.3.3.中代码相同

    // 进行 90 个时间步的模拟更新
    // 由于 tex 是全局并且有界的，因此必须通过一个标志来选择每次迭代中哪个是输入/输出
    volatile bool dstOut = true;
    for (int i = 0; i < 90; i++) {
        float *in, *out;
        if (dstOut) {
            in  = d->dev_inSrc;
            out = d->dev_outSrc;
        } else {
            out = d->dev_inSrc;
            in  = d->dev_outSrc;
        }
        // 将固定热源值拷贝至输入网格
        copy_const_kernel<<<blocks, threads>>>(in);
        // 执行热扩散更新，使用纹理读取输入数据
        blend_kernel<<<blocks, threads>>>(out, dstOut);
        // 切换输入输出标志
        dstOut = !dstOut;
    }

    //  ... 后面与 7.3.3.中代码相同
}

对热传导函数的最后一个修改就是清理工作，不仅要释放全局缓冲区，还需要清除与纹理的绑定：

void anim_exit( DataBlock *d ) {
    cudaUnbindTexture( texIn );
    cudaUnbindTexture( texOut );
    cudaUnbindTexture( texConstSrc );

    cudaFree( d->dev_inSrc );
    cudaFree( d->dev_outSrc );
    cudaFree( d->dev_constSrc );

    HANDLE_ERROR( cudaEventDestroy( d->start ) );
    HANDLE_ERROR( cudaEventDestroy( d->stop ) );
}

7.3.5. 使用二维纹理内存

修改纹理引用的声明，如未声明，默认的纹理引用都是一维的，因此增加了代表位数的参数 2，表示声明的是一个二维纹理引用。

texture<float,2> texConstSrc;
texture<float,2> texIn;
texture<float,2> texOut;

二维纹理将简化 blend_kernel() 方法的实现，虽然需要将 tex1Dfetch() 调用修改为 tex2D() 调用，但不需要线性化 offset 以计算偏移，可以直接通过 x 和 y 来访问纹理。不需要担心发生溢出，tex2D() 会自动处理边界值情况（x 或 y 小于 0，则返回 0 处的值；如果值大于宽度，则返回位于宽度处的值）。

__global__ void blend_kernel(float *outSrc, const float *inSrc) {
    // ... 将线程和块索引映射为图像坐标
    // 无需计算上下左右邻居的索引和边界条件处理

    // 更新当前像素的值
    float t, l, c, r, b;
    if (dstOut) {
        t = tex2D(texIn,x,y-1);
        l = tex2D(texIn,x-1,y);
        c = tex2D(texIn,x,y);
        r = tex2D(texIn,x+1,y);
        b = tex2D(texIn,x,y+1);
    } else {
        t = tex2D(texOut,x,y-1);
        l = tex2D(texOut,x-1,y);
        c = tex2D(texOut,x,y);
        r = tex2D(texOut,x+1,y);
        b = tex2D(texOut,x,y+1);
    }
    dst[offset] = c + SPEED * (t + b + r + l - 4 * c);
}