第7章 CUDA 库 - 站在巨人的肩膀上（CUDA Libraries）

一、导读

1.1 本章学习目标

通过本章的学习，你将能够：

了解 NVIDIA CUDA 生态系统中的主要库及其用途

掌握 cuBLAS 库的基本使用方法进行线性代数运算

学会使用 cuFFT 库进行快速傅里叶变换

掌握 cuRAND 库生成高质量随机数

理解 Thrust 库的并行算法范式

学会使用 cuDNN 加速深度学习操作

了解 NCCL 在多 GPU 通信中的作用

能够根据应用需求选择合适的 CUDA 库

理解库函数与自定义内核的性能权衡

掌握将库集成到 CUDA 程序中的最佳实践

1.2 为什么本章如此重要

编写高效的 CUDA 内核需要深入的专业知识和大量的优化工作。幸运的是，NVIDIA 和其他开发者已经创建了许多高质量的库，涵盖了常见的计算任务。这些库经过精心优化，通常在各种硬件配置上都能提供接近峰值的性能。

使用 CUDA 库的优势包括：

性能保证：库函数经过专业优化，通常优于手写的通用实现
开发效率：无需从头实现常见算法
可维护性：库 API 稳定，减少代码复杂度
跨代兼容：库会自动利用新架构的特性

然而，库并非万能解决方案。理解库的适用场景、限制和最佳使用方法，是成为高效 CUDA 程序员的关键。本章将带你探索主要的 CUDA 库，并展示如何在实际项目中有效利用它们。

库在软件生态中的地位：

在现代软件开发中，库是不可或缺的组成部分。它们封装了复杂的算法和优化技巧，使程序员能够专注于应用逻辑而非底层实现。CUDA 库生态系统尤其丰富，涵盖了从基础线性代数到高级深度学习的广泛领域。

对于科学计算和工程应用，正确使用 CUDA 库可以带来数量级的性能提升。例如，cuBLAS 的矩阵乘法性能通常是手写实现的 2-5 倍，因为它利用了 Tensor Core 等专用硬件单元。cuFFT 的 FFT 实现经过了多年的优化，能够自动选择最优算法。

对于深度学习应用，cuDNN 提供了高度优化的卷积、池化、归一化等操作，是 TensorFlow、PyTorch 等框架的底层依赖。理解 cuDNN 的工作原理有助于更好地调优深度学习模型。

1.3 核心问题引导

在深入学习之前，请思考以下问题：

为什么应该优先使用库而不是自己编写内核？

cuBLAS 支持哪些线性代数操作？

FFT 在哪些应用中至关重要？

为什么 GPU 上的随机数生成具有挑战性？

Thrust 与 STL 有什么相似之处？

cuDNN 为深度学习提供了哪些加速？

多 GPU 通信面临哪些挑战？

如何将多个库组合使用？

二、核心概念详解

7.1 CUDA 库生态系统概览

NVIDIA 提供了丰富的 CUDA 加速库，覆盖了科学计算、信号处理、机器学习等多个领域：

库	用途	典型应用
cuBLAS	基本线性代数子程序	矩阵乘法、向量运算
cuFFT	快速傅里叶变换	信号处理、图像处理
cuRAND	随机数生成	蒙特卡洛模拟、加密
cuSOLVER	线性方程组求解	科学计算、优化
cuSPARSE	稀疏矩阵运算	图算法、机器学习
cuDNN	深度神经网络	深度学习训练和推理
Thrust	并行算法模板库	排序、归约、扫描
NCCL	多 GPU 通信	分布式训练
cuTENSOR	张量运算	高性能深度学习
cuGRAPH	图算法	社交网络分析、推荐系统

选择库的原则：

功能匹配：首先确认库是否提供所需的功能

性能要求：评估库的性能是否满足应用需求

易用性：考虑学习曲线和集成难度

许可证：检查库的许可证是否与项目兼容

社区支持：活跃的社区可以提供更好的帮助

7.2 cuBLAS - 线性代数库

概述

cuBLAS（CUDA Basic Linear Algebra Subprograms）实现了 BLAS 标准，提供了三个级别的线性代数操作：

Level 1：向量-向量操作（点积、范数等），O(n) 复杂度
Level 2：矩阵-向量操作，O(n²) 复杂度
Level 3：矩阵-矩阵操作（GEMM 等），O(n³) 复杂度

BLAS 是科学计算的基础，许多高级算法（如 LAPACK、ScaLAPACK）都依赖于 BLAS 实现。cuBLAS 提供了 GPU 加速的 BLAS 实现，性能远超 CPU 版本。

基本用法

cpp

#include <cublas_v2.h>
#include <stdio.h>

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

#define CUBLAS_CHECK(call) \
    do { \
        cublasStatus_t status = call; \
        if (status != CUBLAS_STATUS_SUCCESS) { \
            printf("cuBLAS error at %s:%d\n", __FILE__, __LINE__); \
            exit(1); \
        } \
    } while(0)

int main() {
    int n = 1024;
    size_t size = n * n * sizeof(float);
    
    // 分配主机内存
    float *h_A = (float *)malloc(size);
    float *h_B = (float *)malloc(size);
    float *h_C = (float *)malloc(size);
    
    // 初始化矩阵
    for (int i = 0; i < n * n; i++) {
        h_A[i] = 1.0f;
        h_B[i] = 2.0f;
    }
    
    // 分配设备内存
    float *d_A, *d_B, *d_C;
    CUDA_CHECK(cudaMalloc((void **)&d_A, size));
    CUDA_CHECK(cudaMalloc((void **)&d_B, size));
    CUDA_CHECK(cudaMalloc((void **)&d_C, size));
    
    CUDA_CHECK(cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice));
    
    // 创建 cuBLAS 句柄
    cublasHandle_t handle;
    CUBLAS_CHECK(cublasCreate(&handle));
    
    // 执行矩阵乘法 C = alpha * A * B + beta * C
    float alpha = 1.0f;
    float beta = 0.0f;
    
    // cuBLAS 使用列主序，注意参数顺序
    CUBLAS_CHECK(cublasSgemm(handle,
                              CUBLAS_OP_N, CUBLAS_OP_N,  // 不转置
                              n, n, n,                    // m, n, k
                              &alpha,                      // alpha
                              d_A, n,                      // A, lda
                              d_B, n,                      // B, ldb
                              &beta,                       // beta
                              d_C, n));                    // C, ldc
    
    // 复制结果回主机
    CUDA_CHECK(cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost));
    
    printf("C[0] = %.1f (expected: %.1f)\n", h_C[0], (float)n * 2.0f);
    
    // 清理
    CUBLAS_CHECK(cublasDestroy(handle));
    CUDA_CHECK(cudaFree(d_A));
    CUDA_CHECK(cudaFree(d_B));
    CUDA_CHECK(cudaFree(d_C));
    free(h_A);
    free(h_B);
    free(h_C);
    
    return 0;
}

常用函数

函数	描述
`cublasSgemm`	单精度矩阵乘法
`cublasDgemm`	双精度矩阵乘法
`cublasSgemv`	单精度矩阵-向量乘法
`cublasSnrm2`	单精度向量 L2 范数
`cublasSdot`	单精度向量点积
`cublasSaxpy`	单精度 A*X+Y 操作

列主序 vs 行主序：

cuBLAS 使用列主序（Column-Major），这与 Fortran 和 MATLAB 一致，但与 C/C++ 的行主序不同。在使用 cuBLAS 时，需要注意这一点：

如果数据是行主序存储的，可以通过交换参数来正确计算
或者在传输数据时进行转置

批处理 GEMM：

cuBLAS 提供了批处理 GEMM 函数，可以一次性执行多个小型矩阵乘法：

cpp

cublasSgemmBatched(handle, transa, transb, m, n, k, 
                   &alpha, dA_array, lda, dB_array, ldb, 
                   &beta, dC_array, ldc, batchCount);

这对于深度学习中的小批量处理非常有用。

7.3 cuFFT - 快速傅里叶变换库

概述

cuFFT 实现了快速傅里叶变换（FFT），支持：

1D、2D、3D FFT
实数到复数、复数到复数变换
批量处理
多 GPU 扩展

FFT 是信号处理、图像处理、科学计算中的核心算法。它将时域/空域信号转换为频域表示，使得滤波、卷积等操作更加高效。

基本用法

cpp

#include <cufft.h>
#include <stdio.h>
#include <complex.h>

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

#define CUFFT_CHECK(call) \
    do { \
        cufftResult result = call; \
        if (result != CUFFT_SUCCESS) { \
            printf("cuFFT error at %s:%d\n", __FILE__, __LINE__); \
            exit(1); \
        } \
    } while(0)

int main() {
    int nx = 1024;
    int ny = 1024;
    
    // 分配复数数组
    int totalSize = nx * ny;
    size_t bytes = totalSize * sizeof(cufftComplex);
    
    cufftComplex *h_data = (cufftComplex *)malloc(bytes);
    cufftComplex *d_data;
    
    // 初始化数据
    for (int i = 0; i < totalSize; i++) {
        h_data[i].x = sin(2.0 * M_PI * i / 64.0);
        h_data[i].y = 0.0f;
    }
    
    CUDA_CHECK(cudaMalloc((void **)&d_data, bytes));
    CUDA_CHECK(cudaMemcpy(d_data, h_data, bytes, cudaMemcpyHostToDevice));
    
    // 创建 cuFFT 计划
    cufftHandle plan;
    CUFFT_CHECK(cufftPlan2d(&plan, nx, ny, CUFFT_C2C));
    
    // 执行正向 FFT
    CUFFT_CHECK(cufftExecC2C(plan, d_data, d_data, CUFFT_FORWARD));
    
    // 执行逆向 FFT
    CUFFT_CHECK(cufftExecC2C(plan, d_data, d_data, CUFFT_INVERSE));
    
    // 缩放结果（cuFFT 不自动归一化）
    float scale = 1.0f / (nx * ny);
    
    // 复制回主机
    CUDA_CHECK(cudaMemcpy(h_data, d_data, bytes, cudaMemcpyDeviceToHost));
    
    printf("FFT completed. First element: (%.3f, %.3f)\n", 
           h_data[0].x * scale, h_data[0].y * scale);
    
    // 清理
    CUFFT_CHECK(cufftDestroy(plan));
    CUDA_CHECK(cudaFree(d_data));
    free(h_data);
    
    return 0;
}

应用场景

图像滤波和卷积：FFT 可以将卷积转换为逐点乘法，大幅降低计算复杂度
音频信号处理：频谱分析、均衡器、降噪
偏微分方程求解：谱方法
量子化学：电子结构计算

性能提示：

FFT 的性能对数据尺寸敏感。2 的幂次尺寸通常最快

使用 cuFFT 计划缓存避免重复创建计划的开销

对于批量 FFT，使用 cufftExecC2C 的批处理版本

7.4 cuRAND - 随机数生成库

概述

cuRAND 提供了高质量的伪随机数和准随机数生成器：

伪随机数：基于 XORWOW、MRG32k3a、Philox 等算法
准随机数：基于 Sobol、Scrambled Sobol 序列

GPU 上的随机数生成具有挑战性，因为：

需要保证统计质量

需要支持大规模并行生成

需要能够重现结果（用于调试）

基本用法

cpp

#include <curand.h>
#include <stdio.h>

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

#define CURAND_CHECK(call) \
    do { \
        curandStatus_t status = call; \
        if (status != CURAND_STATUS_SUCCESS) { \
            printf("cuRAND error at %s:%d\n", __FILE__, __LINE__); \
            exit(1); \
        } \
    } while(0)

int main() {
    int n = 1 << 20;
    size_t size = n * sizeof(float);
    
    float *d_random;
    CUDA_CHECK(cudaMalloc((void **)&d_random, size));
    
    // 创建生成器
    curandGenerator_t gen;
    CURAND_CHECK(curandCreateGenerator(&gen, CURAND_RNG_PSEUDO_DEFAULT));
    
    // 设置种子
    CURAND_CHECK(curandSetPseudoRandomGeneratorSeed(gen, 1234ULL));
    
    // 生成均匀分布的随机数 [0, 1)
    CURAND_CHECK(curandGenerateUniform(gen, d_random, n));
    
    // 复制到主机验证
    float *h_random = (float *)malloc(size);
    CUDA_CHECK(cudaMemcpy(h_random, d_random, size, cudaMemcpyDeviceToHost));
    
    printf("First 10 random numbers:\n");
    for (int i = 0; i < 10; i++) {
        printf("%.6f ", h_random[i]);
    }
    printf("\n");
    
    // 生成正态分布的随机数
    CURAND_CHECK(curandGenerateNormal(gen, d_random, n, 0.0f, 1.0f));
    CUDA_CHECK(cudaMemcpy(h_random, d_random, size, cudaMemcpyDeviceToHost));
    
    printf("First 10 normal distributed numbers:\n");
    for (int i = 0; i < 10; i++) {
        printf("%.6f ", h_random[i]);
    }
    printf("\n");
    
    // 清理
    CURAND_CHECK(curandDestroyGenerator(gen));
    CUDA_CHECK(cudaFree(d_random));
    free(h_random);
    
    return 0;
}

应用场景

蒙特卡洛模拟：金融定价、风险评估
机器学习：权重初始化、数据增强
密码学：密钥生成
游戏：随机事件生成

生成器选择：

生成器	特点	适用场景
XORWOW	快速，中等质量	一般用途
MRG32k3a	高质量，较慢	科学模拟
Philox	高质量，快速	深度学习
Sobol	准随机，低差异	数值积分

7.5 Thrust - 并行算法库

概述

Thrust 是一个类似 C++ STL 的并行算法库，提供了：

排序（sort）
归约（reduce）
扫描（scan）
变换（transform）
搜索（find）

Thrust 的设计目标是提供简洁的 API，同时保持高性能。它自动选择最优的实现策略，并根据数据类型和操作进行调整。

基本用法

cpp

#include <thrust/device_vector.h>
#include <thrust/sort.h>
#include <thrust/reduce.h>
#include <thrust/transform.h>
#include <thrust/fill.h>
#include <thrust/copy.h>
#include <iostream>
#include <functional>

int main() {
    const int N = 1000000;
    
    // 创建设备向量
    thrust::device_vector<float> data(N);
    
    // 填充数据
    thrust::sequence(data.begin(), data.end());
    
    // 归约求和
    float sum = thrust::reduce(data.begin(), data.end(), 0.0f, thrust::plus<float>());
    std::cout << "Sum: " << sum << std::endl;
    
    // 变换：每个元素乘以 2
    thrust::transform(data.begin(), data.end(), data.begin(), 
                      thrust::placeholders::_1 * 2.0f);
    
    // 排序（降序）
    thrust::sort(data.begin(), data.end(), thrust::greater<float>());
    
    // 查找最大值
    float maxVal = data[0];
    std::cout << "Max value: " << maxVal << std::endl;
    
    // 计数满足条件的元素
    int count = thrust::count_if(data.begin(), data.end(), 
                                  thrust::placeholders::_1 > N);
    std::cout << "Elements > N: " << count << std::endl;
    
    return 0;
}

优势

简洁的 API：熟悉 STL 的程序员可以快速上手
自动优化：Thrust 自动选择最优的实现策略
类型安全：编译时类型检查
可扩展：支持自定义函子和迭代器

常见算法：

算法	功能	复杂度
sort	排序	O(n log n)
reduce	归约	O(n)
scan	前缀和	O(n)
transform	元素级变换	O(n)
copy_if	条件复制	O(n)

7.6 cuDNN - 深度神经网络库

概述

cuDNN 提供了深度学习原语的高度优化实现：

卷积（Convolution）
池化（Pooling）
激活函数（Activation）
归一化（Normalization）
RNN/LSTM
Transformer 相关操作

cuDNN 是 TensorFlow、PyTorch 等深度学习框架的底层依赖。它针对不同的 GPU 架构进行了深度优化，能够自动选择最优算法。

基本用法概念

cpp

// cuDNN 的使用模式：
// 1. 创建句柄
cudnnHandle_t handle;
cudnnCreate(&handle);

// 2. 创建张量描述符
cudnnTensorDescriptor_t inputDesc;
cudnnCreateTensorDescriptor(&inputDesc);
cudnnSetTensor4dDescriptor(inputDesc, ...);

// 3. 创建操作描述符（如卷积）
cudnnConvolutionDescriptor_t convDesc;
cudnnCreateConvolutionDescriptor(&convDesc);

// 4. 查询最优算法
cudnnConvolutionFwdAlgo_t algo;
cudnnGetConvolutionForwardAlgorithm(..., &algo);

// 5. 执行操作
cudnnConvolutionForward(handle, ..., inputDesc, d_input, ..., d_output);

// 6. 清理
cudnnDestroyTensorDescriptor(inputDesc);
cudnnDestroyConvolutionDescriptor(convDesc);
cudnnDestroy(handle);

性能提示：

使用 cudnnGetConvolutionForwardAlgorithm 查询最优算法

启用 Tensor Core 加速（使用 FP16 或 TF32）

选择合适的工作区大小以平衡性能和内存

7.7 NCCL - 多 GPU 通信库

概述

NCCL（NVIDIA Collective Communications Library）优化了多 GPU 和多节点之间的通信：

AllReduce
Broadcast
Reduce
AllGather
ReduceScatter

NCCL 是分布式深度学习训练的核心组件，用于在多个 GPU 之间同步梯度。

基本用法

cpp

#include <nccl.h>

// 初始化
ncclComm_t comm;
ncclCommInitRank(&comm, numRanks, uniqueId, rank);

// 执行 AllReduce
ncclAllReduce(sendBuf, recvBuf, count, ncclFloat, ncclSum, comm, stream);

// 清理
ncclCommDestroy(comm);

应用场景

分布式深度学习训练：同步梯度
大规模科学计算：聚合结果
数据并行处理：合并分区数据

三、程序员的视角

7.1 在实际项目中的应用

科学计算流水线

数据加载 → FFT (cuFFT) → 矩阵运算 (cuBLAS) → 统计分析 (Thrust) → 结果输出

深度学习训练

数据预处理 → 前向传播 (cuDNN) → 损失计算 → 反向传播 (cuDNN) → 参数更新

金融模拟

随机数生成 (cuRAND) → 路径模拟 → 统计聚合 (Thrust) → 风险分析

7.2 常见陷阱与误区

误区一：库总是最快的

事实：对于特定问题规模或数据分布，手写内核可能更快。但库提供了良好的基线性能。

误区二：忽略库的资源需求

事实：某些库操作可能需要临时缓冲区，应合理管理内存。

陷阱：数据类型不匹配

确保使用正确的精度版本（S=单精度，D=双精度）。

陷阱：忘记清理资源

库句柄和描述符需要显式销毁。

7.3 与其他层级的联系

第4章：库函数内部也遵循内存访问优化原则
第6章：库调用可以在流中异步执行
第8章：库通常是性能优化的起点

四、代码实战

7.1 完整的 cuBLAS 矩阵乘法示例

cpp

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

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

#define CUBLAS_CHECK(call) \
    do { \
        cublasStatus_t status = call; \
        if (status != CUBLAS_STATUS_SUCCESS) { \
            printf("cuBLAS error\n"); \
            exit(1); \
        } \
    } while(0)

// CPU 参考实现
void matmulCPU(float *A, float *B, float *C, int n) {
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < n; j++) {
            float sum = 0.0f;
            for (int k = 0; k < n; k++) {
                sum += A[i * n + k] * B[k * n + j];
            }
            C[i * n + j] = sum;
        }
    }
}

int main() {
    int n = 1024;
    size_t size = n * n * sizeof(float);
    
    // 分配主机内存
    float *h_A = (float *)malloc(size);
    float *h_B = (float *)malloc(size);
    float *h_C_cublas = (float *)malloc(size);
    float *h_C_cpu = (float *)malloc(size);
    
    // 初始化
    srand(time(NULL));
    for (int i = 0; i < n * n; i++) {
        h_A[i] = (float)rand() / RAND_MAX;
        h_B[i] = (float)rand() / RAND_MAX;
    }
    
    // 分配设备内存
    float *d_A, *d_B, *d_C;
    CUDA_CHECK(cudaMalloc((void **)&d_A, size));
    CUDA_CHECK(cudaMalloc((void **)&d_B, size));
    CUDA_CHECK(cudaMalloc((void **)&d_C, size));
    
    CUDA_CHECK(cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice));
    CUDA_CHECK(cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice));
    
    // 创建 cuBLAS 句柄
    cublasHandle_t handle;
    CUBLAS_CHECK(cublasCreate(&handle));
    
    // 计时
    cudaEvent_t start, stop;
    CUDA_CHECK(cudaEventCreate(&start));
    CUDA_CHECK(cudaEventCreate(&stop));
    
    float alpha = 1.0f, beta = 0.0f;
    
    CUDA_CHECK(cudaEventRecord(start));
    CUBLAS_CHECK(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N,
                              n, n, n, &alpha, d_A, n, d_B, n, &beta, d_C, n));
    CUDA_CHECK(cudaEventRecord(stop));
    CUDA_CHECK(cudaEventSynchronize(stop));
    
    float ms = 0;
    CUDA_CHECK(cudaEventElapsedTime(&ms, start, stop));
    printf("cuBLAS GEMM time: %.3f ms\n", ms);
    
    // 获取结果
    CUDA_CHECK(cudaMemcpy(h_C_cublas, d_C, size, cudaMemcpyDeviceToHost));
    
    // CPU 参考（仅对小矩阵）
    if (n <= 256) {
        clock_t cpuStart = clock();
        matmulCPU(h_A, h_B, h_C_cpu, n);
        clock_t cpuEnd = clock();
        double cpuTime = ((double)(cpuEnd - cpuStart)) / CLOCKS_PER_SEC * 1000;
        printf("CPU GEMM time: %.3f ms\n", cpuTime);
        
        // 验证
        float maxDiff = 0.0f;
        for (int i = 0; i < n * n; i++) {
            float diff = fabs(h_C_cublas[i] - h_C_cpu[i]);
            if (diff > maxDiff) maxDiff = diff;
        }
        printf("Max difference: %.6f\n", maxDiff);
    }
    
    // 清理
    CUBLAS_CHECK(cublasDestroy(handle));
    CUDA_CHECK(cudaFree(d_A));
    CUDA_CHECK(cudaFree(d_B));
    CUDA_CHECK(cudaFree(d_C));
    CUDA_CHECK(cudaEventDestroy(start));
    CUDA_CHECK(cudaEventDestroy(stop));
    free(h_A);
    free(h_B);
    free(h_C_cublas);
    free(h_C_cpu);
    
    return 0;
}

7.2 cuFFT 信号处理示例

cpp

#include <cufft.h>
#include <stdio.h>
#include <math.h>

#define CUDA_CHECK(call) \
    do { \
        cudaError_t err = call; \
        if (err != cudaSuccess) { \
            printf("CUDA error: %s\n", cudaGetErrorString(err)); \
            exit(1); \
        } \
    } while(0)

#define CUFFT_CHECK(call) \
    do { \
        cufftResult result = call; \
        if (result != CUFFT_SUCCESS) { \
            printf("cuFFT error\n"); \
            exit(1); \
        } \
    } while(0)

int main() {
    int N = 1024;
    float sampleRate = 1000.0f; // Hz
    
    // 生成测试信号：两个正弦波的叠加
    cufftComplex *h_signal = (cufftComplex *)malloc(N * sizeof(cufftComplex));
    float freq1 = 50.0f;
    float freq2 = 120.0f;
    
    for (int i = 0; i < N; i++) {
        float t = (float)i / sampleRate;
        h_signal[i].x = sin(2 * M_PI * freq1 * t) + 0.5f * sin(2 * M_PI * freq2 * t);
        h_signal[i].y = 0.0f;
    }
    
    // 分配设备内存
    cufftComplex *d_signal;
    CUDA_CHECK(cudaMalloc((void **)&d_signal, N * sizeof(cufftComplex)));
    CUDA_CHECK(cudaMemcpy(d_signal, h_signal, N * sizeof(cufftComplex), 
                          cudaMemcpyHostToDevice));
    
    // 创建 FFT 计划
    cufftHandle plan;
    CUFFT_CHECK(cufftPlan1d(&plan, N, CUFFT_C2C, 1));
    
    // 执行 FFT
    CUFFT_CHECK(cufftExecC2C(plan, d_signal, d_signal, CUFFT_FORWARD));
    
    // 获取结果
    CUDA_CHECK(cudaMemcpy(h_signal, d_signal, N * sizeof(cufftComplex), 
                          cudaMemcpyDeviceToHost));
    
    // 计算幅度谱
    printf("Frequency spectrum (first 20 bins):\n");
    for (int i = 0; i < 20; i++) {
        float magnitude = sqrt(h_signal[i].x * h_signal[i].x + 
                               h_signal[i].y * h_signal[i].y);
        float freq = (float)i * sampleRate / N;
        printf("  %.1f Hz: %.3f\n", freq, magnitude);
    }
    
    // 清理
    CUFFT_CHECK(cufftDestroy(plan));
    CUDA_CHECK(cudaFree(d_signal));
    free(h_signal);
    
    return 0;
}

7.3 Thrust 并行算法示例

cpp

#include <thrust/device_vector.h>
#include <thrust/host_vector.h>
#include <thrust/sort.h>
#include <thrust/reduce.h>
#include <thrust/transform.h>
#include <thrust/random.h>
#include <thrust/iterator/counting_iterator.h>
#include <iostream>
#include <functional>

// 自定义函子
struct square_functor {
    __host__ __device__
    float operator()(const float& x) const {
        return x * x;
    }
};

int main() {
    const int N = 1000000;
    
    // 生成随机数据
    thrust::device_vector<float> data(N);
    thrust::default_random_engine rng(1234);
    thrust::uniform_real_distribution<float> dist(-10.0f, 10.0f);
    
    for (int i = 0; i < N; i++) {
        data[i] = dist(rng);
    }
    
    // 计算均方根
    float sumOfSquares = thrust::transform_reduce(
        data.begin(), data.end(), 
        0.0f, 
        thrust::plus<float>(), 
        square_functor());
    float rms = sqrt(sumOfSquares / N);
    std::cout << "RMS: " << rms << std::endl;
    
    // 排序
    thrust::sort(data.begin(), data.end());
    
    // 计算中位数
    float median = data[N / 2];
    std::cout << "Median: " << median << std::endl;
    
    // 统计范围内的元素
    int inRange = thrust::count_if(
        data.begin(), data.end(),
        [] __host__ __device__ (float x) {
            return x >= -1.0f && x <= 1.0f;
        });
    std::cout << "Elements in [-1, 1]: " << inRange << std::endl;
    
    return 0;
}

五、深入思考

5.1 历史背景

BLAS 库的历史可以追溯到 1970 年代，成为了科学计算的事实标准。NVIDIA 将这一标准移植到 GPU 上，创建了 cuBLAS。类似的，其他库也是将成熟的算法范式适配到 GPU 架构。

5.2 前沿方向

库的组合优化：自动融合多个库调用
自适应精度：根据需求动态调整精度
稀疏加速：针对稀疏数据的专用库

5.3 延伸阅读

NVIDIA CUDA Library Documentation
BLAS 标准规范
Thrust 用户指南

六、实验笔记

实验 1：比较 cuBLAS 与手写内核的性能

目标：评估库函数相对于自定义实现的优势。

步骤：

实现简单的矩阵乘法内核

使用相同规模的 cuBLAS GEMM

比较执行时间和准确性

实验 2：探索 Thrust 的各种算法

目标：熟悉 Thrust 提供的并行算法。

步骤：

实现排序、归约、扫描等操作

与 CPU 版本比较

分析性能差异

七、本章小结

本章介绍了 CUDA 生态系统中的主要库。我们学习了 cuBLAS、cuFFT、cuRAND、Thrust 等库的基本使用方法。

关键要点包括：

库提供了经过优化的实现，应优先使用
cuBLAS 是线性代数计算的标准
cuFFT 高效处理频域变换
cuRAND 提供高质量随机数
Thrust 简化并行算法编程

在下一章中，我们将探讨性能优化的高级技巧。

八、习题

列举至少五个 CUDA 库及其主要用途。

cuBLAS 中的 Sgemm 函数的参数含义是什么？

为什么 FFT 在信号处理中如此重要？

cuRAND 支持哪些随机数分布？

Thrust 与 C++ STL 有什么相似之处？

编写代码演示如何使用 cuBLAS 进行向量点积计算。

解释 cuFFT 中正逆变换的归一化问题。

什么场景下应该考虑手写内核而不是使用库？

NCCL 在分布式训练中扮演什么角色？

设计一个使用多个库组合的科学计算流程。