【CUDA】CUBLAS

在深入了解之前，提前运行预热（warmup）和基准测试（benchmark runs） 是获得准确执行时间的关键。如果不进行预热运行，cuBLAS 的首次运行会有较大的开销（大约 45 毫秒），会导致结果偏差。基准测试能够更准确地计算出平均执行时间。

cuBLAS：CUDA 基本线性代数子程序库

简介：cuBLAS 是 NVIDIA 的 GPU 加速线性代数运算库，广泛应用于 人工智能（AI） 和 高性能计算（HPC） 领域。
功能：提供了行业标准的 BLAS 和 GEMM（矩阵乘法）API，并支持高度优化的融合操作（fusion），充分发挥 NVIDIA GPU 的性能。
使用提示：在使用 GEMM 操作时，需要特别注意矩阵的存储布局（行优先或列优先），例如参考这里。

cuBLASLt：轻量级扩展

简介：cuBLASLt 是 cuBLAS 的轻量级扩展，提供了更灵活的 API，专注于提升特定工作负载（如深度学习模型）的性能。
特点：
- 如果单个内核无法处理问题，cuBLASLt 会将问题分解为多个子问题，并在每个子问题上运行内核。
- 支持混合精度计算，如 fp16、fp8 和 int8，可显著提升深度学习的推理速度。

cuBLASXt：支持多 GPU 扩展

简介：cuBLASXt 是 cuBLAS 的扩展版，主要针对超大规模计算，支持多 GPU 运行。
特点：
- 多 GPU 支持：能够将 BLAS 操作分布到多块 GPU 上，适合处理需要扩展 GPU 计算的大型数据集。
- 线程安全：支持多线程并发执行，在不同的 GPU 上同时运行多个 BLAS 操作。
- 适用场景：特别适用于超出单块 GPU 显存限制的大规模线性代数问题。
缺点：由于需要在 主板 DRAM 和 GPU VRAM 之间频繁传输数据，会造成内存带宽瓶颈，导致计算速度较慢。

cuBLASDx：设备端扩展（未在课程中使用）

简介：cuBLASDx 是一个设备端 API 扩展，用于直接在 CUDA 内核中执行 BLAS 计算。
特点：
- 通过融合（fusion）数值操作，进一步降低延迟，提升应用程序性能。
- 注意：cuBLASDx 并不包含在 CUDA Toolkit 中，需要单独下载。

CUTLASS：CUDA 模板线性代数子程序库

简介：cuBLAS 及其变体主要在 主机（host）端 运行，而 CUTLASS 提供了模板库，允许开发者实现高度自定义和优化的线性代数操作。
特点：
- 支持在 CUDA 中轻松融合矩阵运算。
- 对于深度学习，矩阵乘法是最重要的操作，而 cuBLAS 无法轻松实现复杂操作的融合。
补充说明：
- CUTLASS 并未用于实现 Flash Attention，后者是通过高度优化的 CUDA 内核实现的（详见论文）。

总结与应用场景

库	特点与适用场景
cuBLAS	高性能线性代数运算，适用于 AI 和 HPC。
cuBLASLt	灵活 API 和混合精度支持，专注深度学习工作负载。
cuBLASXt	多 GPU 支持，适合超大规模计算，但受限于内存带宽瓶颈。
cuBLASDx	设备端 API，可在 CUDA 内核中实现 BLAS 操作，进一步优化延迟。
CUTLASS	提供模板库，支持复杂运算融合，深度学习中高效矩阵运算的选择。

通过根据工作负载的需求选择合适的库，可以充分利用 NVIDIA GPU 的计算能力。

cuBLAS示例

#include <cublas_v2.h>
#include <cuda_fp16.h>
#include <cuda_runtime.h>
#include <stdio.h>#define M 3
#define K 4
#define N 2#define CHECK_CUDA(call)                                                                               \{                                                                                                  \cudaError_t err = call;                                                                        \if (err != cudaSuccess) {                                                                      \fprintf(stderr, "CUDA error in %s:%d: %s\n", __FILE__, __LINE__, cudaGetErrorString(err)); \exit(EXIT_FAILURE);                                                                        \}                                                                                              \}#define CHECK_CUBLAS(call)                                                              \{                                                                                   \cublasStatus_t status = call;                                                   \if (status != CUBLAS_STATUS_SUCCESS) {                                          \fprintf(stderr, "cuBLAS error in %s:%d: %d\n", __FILE__, __LINE__, status); \exit(EXIT_FAILURE);                                                         \}                                                                               \}#undef PRINT_MATRIX
#define PRINT_MATRIX(mat, rows, cols)                                       \for (int i = 0; i < rows; i++) {                                        \for (int j = 0; j < cols; j++) printf("%8.3f ", mat[i * cols + j]); \printf("\n");                                                       \}                                                                       \printf("\n");void CpuMatmul(float* A, float* B, float* C) {for (int i = 0; i < M; ++i) {for (int j = 0; j < N; ++j) {float sum = 0;for (int k = 0; k < K; ++k) {sum += A[i * K + k] * B[k * N + j];}C[i * N + j] = sum;}}
}int main() {float A[M * K] = {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f, 9.0f, 10.0f, 11.0f, 12.0f};float B[K * N] = {1.0f, 2.0f, 3.0f, 4.0f, 5.0f, 6.0f, 7.0f, 8.0f};float C_cpu[M * N], C_cublas_s[M * N], C_cublas_h[M * N];// CPU matmulCpuMatmul(A, B, C_cpu);cublasHandle_t handle;CHECK_CUBLAS(cublasCreate(&handle));float *d_A, *d_B, *d_C;CHECK_CUDA(cudaMalloc(&d_A, M * K * sizeof(float)));CHECK_CUDA(cudaMalloc(&d_B, K * N * sizeof(float)));CHECK_CUDA(cudaMalloc(&d_C, M * N * sizeof(float)));CHECK_CUDA(cudaMemcpy(d_A, A, M * K * sizeof(float), cudaMemcpyHostToDevice));CHECK_CUDA(cudaMemcpy(d_B, B, K * N * sizeof(float), cudaMemcpyHostToDevice));// cuBLAS SGEMM(单精度)float alpha = 1.0f, beta = 0.0f;/*cublasSgemm ：矩阵乘法公式 C = alpha x op(A) @ op(B) + beta x CcublasStatus_t cublasSgemm(cublasHandle_t handle,      // cuBLAS上下文// CUBLAS_OP_N：不转置（默认）。 CUBLAS_OP_T：转置。CUBLAS_OP_C：共轭转置。cublasOperation_t transa,  //指定矩阵A是否转置,cublasOperation_t transb,  // 指定矩阵B是否转置int m,                     // 矩阵C的行数int n,                     // 矩阵C的列数int k,                     // A的列数或转置后行数；B的行数或转置后列数const float *alpha,        // 标量alphaconst float *A,            // 矩阵A// 主维度（leading dimension），表示存储矩阵时每列或每行的跨度,如果矩阵是列主序，则 lda 是矩阵 A 的行数。int lda,                   // A的主维（leading dimension）const float *B,            // 矩阵Bint ldb,                   // B的主维const float *beta,         // 标量betafloat *C,                  // 输出矩阵Cint ldc                    // C的主维);*/// cuBLAS 默认使用列主序存储矩阵。如果输入矩阵是行主序存储，则需要手动调整主维度，或使用技巧重新解释。/*通过技巧将问题转换为计算 (B^T * A^T)^T，从而直接获得期望结果：矩阵 A 和 B 在内存中以行主序存储，等价于将其转置解释为列主序存储。调用 cublasSgemm 时，将矩阵按未转置（CUBLAS_OP_N）处理，令 cuBLAS 按列主序解释输入。调整矩阵的传递顺序和参数，实际计算 (B^T * A^T)^T，最终结果矩阵在内存中直接符合行主序的期望。*/CHECK_CUBLAS(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N));CHECK_CUDA(cudaMemcpy(C_cublas_s, d_C, M * N * sizeof(float), cudaMemcpyDeviceToHost));// cuBLAS HGEMMhalf *d_A_h, *d_B_h, *d_C_h;CHECK_CUDA(cudaMalloc(&d_A_h, M * K * sizeof(half)));CHECK_CUDA(cudaMalloc(&d_B_h, K * N * sizeof(half)));CHECK_CUDA(cudaMalloc(&d_C_h, M * N * sizeof(half)));// Convert to half percision to CPUhalf A_h[M * K], B_h[K * N];for (int i = 0; i < M * K; ++i) {A_h[i] = __float2half(A[i]);}for (int i = 0; i < K * N; ++i) {B_h[i] = __float2half(B[i]);}CHECK_CUDA(cudaMemcpy(d_A_h, A_h, M * K * sizeof(half), cudaMemcpyHostToDevice));CHECK_CUDA(cudaMemcpy(d_B_h, B_h, K * N * sizeof(half), cudaMemcpyHostToDevice));half alpha_h = __float2half(1.0f), bata_h = __float2half(0.0f);CHECK_CUBLAS(cublasHgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha_h, d_B_h, N, d_A_h, K, &bata_h, d_C_h, N));// Copy result back to host and convert to floathalf C_h[M * N];CHECK_CUDA(cudaMemcpy(C_h, d_C_h, M * N * sizeof(half), cudaMemcpyDeviceToHost));for (int i = 0; i < M * N; i++) {C_cublas_h[i] = __half2float(C_h[i]);}// Print resultsprintf("Matrix A (%dx%d):\n", M, K);PRINT_MATRIX(A, M, K);printf("Matrix B (%dx%d):\n", K, N);PRINT_MATRIX(B, K, N);printf("CPU Result (%dx%d):\n", M, N);PRINT_MATRIX(C_cpu, M, N);printf("cuBLAS SGEMM Result (%dx%d):\n", M, N);PRINT_MATRIX(C_cublas_s, M, N);printf("cuBLAS HGEMM Result (%dx%d):\n", M, N);PRINT_MATRIX(C_cublas_h, M, N);// Clean upCHECK_CUDA(cudaFree(d_A));CHECK_CUDA(cudaFree(d_B));CHECK_CUDA(cudaFree(d_C));CHECK_CUDA(cudaFree(d_A_h));CHECK_CUDA(cudaFree(d_B_h));CHECK_CUDA(cudaFree(d_C_h));CHECK_CUBLAS(cublasDestroy(handle));return 0;
}

输出：

Matrix A (3x4):1.000    2.000    3.000    4.000 5.000    6.000    7.000    8.000 9.000   10.000   11.000   12.000 Matrix B (4x2):1.000    2.000 3.000    4.000 5.000    6.000 7.000    8.000 CPU Result (3x2):50.000   60.000 114.000  140.000 178.000  220.000 cuBLAS SGEMM Result (3x2):50.000   60.000 114.000  140.000 178.000  220.000 cuBLAS HGEMM Result (3x2):50.000   60.000 114.000  140.000 178.000  220.000

cuBLASLt示例

#include <cublasLt.h>  // cuBLASLt 是 NVIDIA 的 cuBLAS 库的扩展，提供了更灵活和可配置的矩阵乘法接口。
#include <cuda_fp16.h>  // 提供对 half（FP16）数据类型的支持
#include <cuda_runtime.h>#include <iomanip>
#include <iostream>
#include <vector>#define CHECK_CUDA(call)                                                                                       \do {                                                                                                       \cudaError_t status = call;                                                                             \if (status != cudaSuccess) {                                                                           \std::cerr << "CUDA error at line " << __LINE__ << ": " << cudaGetErrorString(status) << std::endl; \exit(EXIT_FAILURE);                                                                                \}                                                                                                      \} while (0)#define CHECK_CUBLAS(call)                                                                   \do {                                                                                     \cublasStatus_t status = call;                                                        \if (status != CUBLAS_STATUS_SUCCESS) {                                               \std::cerr << "cuBLAS error at line " << __LINE__ << ": " << status << std::endl; \exit(EXIT_FAILURE);                                                              \}                                                                                    \} while (0)void CpuMatmul(const float* A, const float* B, float* C, int M, int N, int K) {for (int i = 0; i < M; ++i) {for (int j = 0; j < N; ++j) {float sum = 0.0f;for (int k = 0; k < K; ++k) {sum += A[i * K + k] * B[k * N + j];}C[i * N + j] = sum;}}
}void PrintMatrix(const float* matrix, int rows, int cols, const char* name) {std::cout << "Matrix " << name << ":" << std::endl;for (int i = 0; i < rows; ++i) {for (int j = 0; j < cols; ++j) {// std::setw(8)用于设置输出字段的宽度为 8 个字符，std::fixed用于设置浮点数的输出格式为固定小数点格式// std::setprecision(2)用于设置浮点数的小数点后的精度为 2 位。std::cout << std::setw(8) << std::fixed << std::setprecision(2) << matrix[i * cols + j] << ' ';}std::cout << std::endl;}
}int main() {const int M = 4, K = 4, N = 4;float h_A[M * K] = {1.0f, 2.0f,  3.0f,  4.0f,  5.0f,  6.0f,  7.0f,  8.0f,9.0f, 10.0f, 11.0f, 12.0f, 13.0f, 14.0f, 15.0f, 16.0f};float h_B[K * N] = {1.0f, 2.0f,  4.0f,  4.0f,  5.0f,  6.0f,  7.0f,  8.0f,9.0f, 10.0f, 11.0f, 12.0f, 17.0f, 18.0f, 19.0f, 20.0f};float h_C_cpu[M * N] = {0};float h_C_gpu_fp32[M * N] = {0};float h_C_gpu_fp16[M * N] = {0};// Print input matricesPrintMatrix(h_A, M, K, "A");PrintMatrix(h_B, K, N, "B");float *d_A_fp32, *d_B_fp32, *d_C_fp32;CHECK_CUDA(cudaMalloc(&d_A_fp32, M * K * sizeof(float)));CHECK_CUDA(cudaMalloc(&d_B_fp32, K * N * sizeof(float)));CHECK_CUDA(cudaMalloc(&d_C_fp32, M * N * sizeof(float)));half *d_A_fp16, *d_B_fp16, *d_C_fp16;CHECK_CUDA(cudaMalloc(&d_A_fp16, M * K * sizeof(half)));CHECK_CUDA(cudaMalloc(&d_B_fp16, K * N * sizeof(half)));CHECK_CUDA(cudaMalloc(&d_C_fp16, M * N * sizeof(half)));CHECK_CUDA(cudaMemcpy(d_A_fp32, h_A, M * K * sizeof(float), cudaMemcpyHostToDevice));CHECK_CUDA(cudaMemcpy(d_B_fp32, h_B, K * N * sizeof(float), cudaMemcpyHostToDevice));// Convert and copy data to device(FP16)std::vector<half> h_A_half(M * K);std::vector<half> h_B_half(K * N);for (int i = 0; i < M * K; ++i) h_A_half[i] = __float2half(h_A[i]);for (int i = 0; i < K * N; ++i) h_B_half[i] = __float2half(h_B[i]);CHECK_CUDA(cudaMemcpy(d_A_fp16, h_A_half.data(), M * K * sizeof(half), cudaMemcpyHostToDevice));CHECK_CUDA(cudaMemcpy(d_B_fp16, h_B_half.data(), K * N * sizeof(half), cudaMemcpyHostToDevice));// Create cuBLAS handlecublasLtHandle_t handle;CHECK_CUBLAS(cublasLtCreate(&handle));// Set up matrix descriptors for FP32/*cublasStatus_t cublasLtMatrixLayoutCreate(cublasLtMatrixLayout_t *matLayout, cudaDataType type, uint64_t rows,uint64_t cols, int64_t ld)矩阵描述符的作用 矩阵描述符告诉 cuBLAS Lt： 矩阵的数据类型：如CUDA_R_32F（单精度浮点）、CUDA_R_16F（半精度浮点）。矩阵的维度：行数和列数。矩阵在内存中的步幅：矩阵的行或列在内存中的存储间隔（通常等于矩阵的列数或行数）。*/cublasLtMatrixLayout_t matA_fp32, matB_fp32, matC_fp32;CHECK_CUBLAS(cublasLtMatrixLayoutCreate(&matA_fp32, CUDA_R_32F, K, M, K));CHECK_CUBLAS(cublasLtMatrixLayoutCreate(&matB_fp32, CUDA_R_32F, N, K, N));CHECK_CUBLAS(cublasLtMatrixLayoutCreate(&matC_fp32, CUDA_R_32F, N, M, N));// Set up matrix descriptors for FP16cublasLtMatrixLayout_t matA_fp16, matB_fp16, matC_fp16;CHECK_CUBLAS(cublasLtMatrixLayoutCreate(&matA_fp16, CUDA_R_16F, K, M, K));  // original MKKCHECK_CUBLAS(cublasLtMatrixLayoutCreate(&matB_fp16, CUDA_R_16F, N, K, N));  // original KNNCHECK_CUBLAS(cublasLtMatrixLayoutCreate(&matC_fp16, CUDA_R_16F, N, M, N));  // original MNN// Set up matrix multiplication descriptor for FP32/*矩阵乘法描述符 cublasLtMatmulDescCreate：Create new matmul operation descriptor.cublasStatus_t cublasLtMatmulDescCreate(cublasLtMatmulDesc_t *matmulDesc, cublasComputeType_t computeType,cudaDataType_t scaleType)通过创建矩阵乘法描述符，用户能够设置：1）矩阵是否需要转置、共轭转置。2）计算的精度。3）加法、激活函数等附加操作。*/cublasLtMatmulDesc_t matmulDesc_fp32;CHECK_CUBLAS(cublasLtMatmulDescCreate(&matmulDesc_fp32, CUBLAS_COMPUTE_32F, CUDA_R_32F));// Set up matrix multiplication descriptor for FP16cublasLtMatmulDesc_t matmulDesc_fp16;CHECK_CUBLAS(cublasLtMatmulDescCreate(&matmulDesc_fp16, CUBLAS_COMPUTE_16F, CUDA_R_16F));// Set matrix operation for A and BcublasOperation_t transa = CUBLAS_OP_N;cublasOperation_t transb = CUBLAS_OP_N;/*cublasLtMatmulDescSetAttribute 用于设置矩阵乘法描述符 (cublasLtMatmulDesc_t)的属性。这个函数非常关键，因为它允许用户根据不同的需求配置矩阵乘法的行为，例如设置加法、转置、矩阵的精度等。通过这种方式，您可以定制和优化矩阵乘法的计算过程。cublasStatus_t cublasLtMatmulDescSetAttribute(cublasLtMatmulDesc_t matmul_desc,  // 矩阵乘法描述符cublasLtMatmulDescAttribute_t attr, // 需要设置的属性const void *value,                 // 属性值size_t size                       // 属性值的大小);CUBLASLT_MATMUL_DESC_TRANSA作用：指定矩阵 A 是否需要转置。值：CUBLAS_OP_N：矩阵 A 不需要转置（即使用原始矩阵 A）。CUBLAS_OP_T：矩阵 A 需要转置（即对矩阵 A 执行转置操作）。CUBLAS_OP_CONJ_T：矩阵 A 需要进行共轭转置（即对矩阵 A 执行共轭转置操作）。*/CHECK_CUBLAS(cublasLtMatmulDescSetAttribute(matmulDesc_fp32, CUBLASLT_MATMUL_DESC_TRANSA, &transa,sizeof(cublasOperation_t)));CHECK_CUBLAS(cublasLtMatmulDescSetAttribute(matmulDesc_fp32, CUBLASLT_MATMUL_DESC_TRANSB, &transb,sizeof(cublasOperation_t)));CHECK_CUBLAS(cublasLtMatmulDescSetAttribute(matmulDesc_fp16, CUBLASLT_MATMUL_DESC_TRANSA, &transa,sizeof(cublasOperation_t)));CHECK_CUBLAS(cublasLtMatmulDescSetAttribute(matmulDesc_fp16, CUBLASLT_MATMUL_DESC_TRANSB, &transb,sizeof(cublasOperation_t)));// Set up alpha and betaconst float alpha = 1.0f;const float beta = 0.0f;// Perform matrix multiplication using cublasLtMatmul(FP32)/*Execute matrix multiplication (D = alpha * op(A) * op(B) + beta * C).cublasStatus_t cublasLtMatmul(cublasLtHandle_t lightHandle,       // 句柄cublasLtMatmulDesc_t computeDesc,   // 矩阵乘法的描述符const void *alpha,const void *A, cublasLtMatrixLayout_t Adesc,    //输入矩阵A,以及矩阵 A 的布局描述符Adescconst void *B, cublasLtMatrixLayout_t Bdesc,const void *beta,const void *C, cublasLtMatrixLayout_t Cdesc,    //输入矩阵 C，以及布局描述符Cdescvoid *D, cublasLtMatrixLayout_t Ddesc,          //输出矩阵 D，以及布局描述符Ddesc// // algo 是指向 cublasLtMatmulAlgo_t结构体的指针，定义了具体的矩阵乘法算法const cublasLtMatmulAlgo_t *algo,// workspace 是指向用于矩阵乘法运算的工作空间的指针,某些算法可能需要额外的内存来存储中间数据void *workspace, size_t workspaceSizeInBytes,cudaStream_t stream //CUDA 流);*/CHECK_CUBLAS(cublasLtMatmul(handle, matmulDesc_fp32, &alpha, d_B_fp32, matB_fp32, d_A_fp32, matA_fp32, &beta,d_C_fp32, matC_fp32, d_C_fp32, matC_fp32, nullptr, nullptr, 0, 0));// half alpha and betaconst half alpha_half = __float2half(1.0f);const half beta_half = __float2half(0.0f);// Perform matrix multiplication using cublasLtMatmul (FP16)CHECK_CUBLAS(cublasLtMatmul(handle, matmulDesc_fp16, &alpha_half, d_B_fp16, matB_fp16, d_A_fp16, matA_fp16,&beta_half, d_C_fp16, matC_fp16, d_C_fp16, matC_fp16, nullptr, nullptr, 0, 0));// Copy results back to hostCHECK_CUDA(cudaMemcpy(h_C_gpu_fp32, d_C_fp32, M * N * sizeof(float), cudaMemcpyDeviceToHost));std::vector<half> h_C_gpu_fp16_half(M * N);CHECK_CUDA(cudaMemcpy(h_C_gpu_fp16_half.data(), d_C_fp16, M * N * sizeof(half), cudaMemcpyDeviceToHost));// Convert half precision results to single precisionfor (int i = 0; i < M * N; ++i) {h_C_gpu_fp16[i] = __half2float(h_C_gpu_fp16_half[i]);}// Perform CPU matrix multiplicationCpuMatmul(h_A, h_B, h_C_cpu, M, N, K);// Print resultsPrintMatrix(h_C_cpu, M, N, "C (CPU)");PrintMatrix(h_C_gpu_fp32, M, N, "C (GPU FP32)");PrintMatrix(h_C_gpu_fp16, M, N, "C (GPU FP16)");// Compare CPU and GPU resultsbool fp32_match = true;bool fp16_match = true;for (int i = 0; i < M * N; ++i) {if (std::abs(h_C_cpu[i] - h_C_gpu_fp32[i]) > 1e-5) {fp32_match = false;}if (std::abs(h_C_cpu[i] - h_C_gpu_fp16[i]) > 1e-2) {  // Increased tolerance for FP16fp16_match = false;}}std::cout << "FP32 Results " << (fp32_match ? "match" : "do not match") << std::endl;std::cout << "FP16 Results " << (fp16_match ? "match" : "do not match") << std::endl;// Clean upCHECK_CUBLAS(cublasLtMatrixLayoutDestroy(matA_fp32));CHECK_CUBLAS(cublasLtMatrixLayoutDestroy(matB_fp32));CHECK_CUBLAS(cublasLtMatrixLayoutDestroy(matC_fp32));CHECK_CUBLAS(cublasLtMatrixLayoutDestroy(matA_fp16));CHECK_CUBLAS(cublasLtMatrixLayoutDestroy(matB_fp16));CHECK_CUBLAS(cublasLtMatrixLayoutDestroy(matC_fp16));CHECK_CUBLAS(cublasLtMatmulDescDestroy(matmulDesc_fp32));CHECK_CUBLAS(cublasLtMatmulDescDestroy(matmulDesc_fp16));CHECK_CUBLAS(cublasLtDestroy(handle));CHECK_CUDA(cudaFree(d_A_fp32));CHECK_CUDA(cudaFree(d_B_fp32));CHECK_CUDA(cudaFree(d_C_fp32));CHECK_CUDA(cudaFree(d_A_fp16));CHECK_CUDA(cudaFree(d_B_fp16));CHECK_CUDA(cudaFree(d_C_fp16));return 0;
}

cuBLASXt示例

cuBLASXt使用起来和cuBLAS类似。

首先和CPU的性能比较：

#include <cublasXt.h>
#include <cublas_v2.h>
#include <cuda_runtime.h>#include <cstdlib>
#include <ctime>
#include <iostream>// 定义矩阵维度，矩阵尺寸为 256 x 256
const int M = 1024 / 4;  // Rows of matrix A and C
const int N = 1024 / 4;  // Columns of matrix B and C
const int K = 1024 / 4;  // Columns of matrix A and rows of matrix B// 定义检查 cuBLAS 调用状态的宏
#define CHECK_CUBLAS(call)                                                           \{                                                                                \cublasStatus_t err = call;                                                   \if (err != CUBLAS_STATUS_SUCCESS) {                                          \std::cerr << "Error in " << #call << ", line " << __LINE__ << std::endl; \exit(1);                                                                 \}                                                                            \}int main() {// 初始化随机数生成器srand(time(0));// 在主机（CPU）上分配内存用于存储矩阵 A, B, Cfloat* A_host = new float[M * K];float* B_host = new float[K * N];float* C_host_cpu = new float[M * N];  // 用于 CPU 计算结果float* C_host_gpu = new float[M * N];  // 用于 GPU 计算结果// 初始化矩阵 A 和 B 的值为随机浮点数for (int i = 0; i < M * K; i++) {A_host[i] = (float)rand() / RAND_MAX;}for (int i = 0; i < K * N; i++) {B_host[i] = (float)rand() / RAND_MAX;}// 在 CPU 上进行矩阵乘法 C = A * Bfloat alpha = 1.0f;  // 缩放因子 alphafloat beta = 0.0f;   // 缩放因子 betafor (int i = 0; i < M; i++) {for (int j = 0; j < N; j++) {C_host_cpu[i * N + j] = 0.0f;  // 初始化 C 矩阵的元素for (int k = 0; k < K; k++) {C_host_cpu[i * N + j] += A_host[i * K + k] * B_host[k * N + j];}}}// 创建 cuBLASXt 句柄cublasXtHandle_t handle;CHECK_CUBLAS(cublasXtCreate(&handle));// 设置 GPU 设备，选择 GPU 0int devices[1] = {0};/*cublasXtDeviceSelect 用于指定 cuBLASXt 将使用哪些 GPU 设备进行计算。cuBLASXt 支持多GPU的并行计算，通过这个函数用户可以灵活选择目标设备。cublasStatus_t cublasXtDeviceSelect(cublasXtHandle_t handle, int nbDevices, const int* deviceIds);*/CHECK_CUBLAS(cublasXtDeviceSelect(handle, 1, devices));// 进行一次 warmup 运行，调用 cuBLASXt 的矩阵乘法接口进行计算/*cublasXtSgemm 是 cuBLASXt 提供的单精度矩阵乘法 (GEMM) 函数。它支持单 GPU 和多 GPU 加速，执行计算公式：C=α⋅op(A)⋅op(B)+β⋅CcublasStatus_t cublasXtSgemm(cublasXtHandle_t handle,cublasOperation_t transA,cublasOperation_t transB,int m,int n,int k,const float* alpha,const float* A,int lda,const float* B,int ldb,const float* beta,float* C,int ldc);m:矩阵 C 的行数(矩阵 A 的行数).n:矩阵 C 的列数(矩阵 B 的列数).k:矩阵 A 的列数(矩阵 B 的行数)。*/CHECK_CUBLAS(cublasXtSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, B_host, N, A_host, K, &beta, C_host_gpu, N));// 比较 CPU 和 GPU 计算结果float max_diff = 1e-4f;  // 最大容许误差for (int i = 0; i < M * N; i++) {float diff = std::abs(C_host_cpu[i] - C_host_gpu[i]);if (diff > max_diff) {std::cout << "i: " << i << " CPU: " << C_host_cpu[i] << ", GPU: " << C_host_gpu[i] << std::endl;}}std::cout << "Maximum difference between CPU and GPU results: " << max_diff << std::endl;// 释放主机内存delete[] A_host;delete[] B_host;delete[] C_host_cpu;delete[] C_host_gpu;return 0;
}

CUBLAS和CUBLAS-XT的性能比较：

#include <cublasXt.h>
#include <cublas_v2.h>
#include <cuda_runtime.h>#include <chrono>
#include <iostream>
#include <vector>#define CHECK_CUDA(call)                                                            \{                                                                               \cudaError_t err = call;                                                     \if (err != cudaSuccess) {                                                   \printf("CUDA error: %s, line %d\n", cudaGetErrorString(err), __LINE__); \exit(1);                                                                \}                                                                           \}
#define CHECK_CUBLAS(call)                                           \{                                                                \cublasStatus_t status = call;                                \if (status != CUBLAS_STATUS_SUCCESS) {                       \printf("CUBLAS error: %d, line %d\n", status, __LINE__); \exit(1);                                                 \}                                                            \}void initMatrix(float *matrix, int rows, int cols) {for (int i = 0; i < rows * cols; ++i) {matrix[i] = static_cast<float>(rand()) / RAND_MAX;}
}bool compareResults(float *result1, float *result2, int size, float tolerance) {for (int i = 0; i < size; ++i) {float diff = std::abs(result1[i] - result2[i]);float max_val = std::max(std::abs(result1[i]), std::abs(result2[i]));if (diff / max_val > tolerance) {std::cout << "Results do not match at index " << i << std::endl;std::cout << "CUBLAS: " << result1[i] << ", CUBLAS-XT: " << result2[i] << std::endl;std::cout << "Relative difference: " << diff / max_val << std::endl;return false;}}return true;
}int main() {int M = 16384;int N = 16384;int K = 16384;size_t size_A = M * K * sizeof(float);size_t size_B = K * N * sizeof(float);size_t size_C = M * N * sizeof(float);float *h_A = (float *)malloc(size_A);float *h_B = (float *)malloc(size_B);float *h_C_cublas = (float *)malloc(size_C);float *h_C_cublasxt = (float *)malloc(size_C);initMatrix(h_A, M, K);initMatrix(h_B, K, N);const int num_runs = 5;std::vector<double> cublas_times;std::vector<double> cublasxt_times;// CUBLAS{cublasHandle_t handle;CHECK_CUBLAS(cublasCreate(&handle));float *d_A, *d_B, *d_C;CHECK_CUDA(cudaMalloc(&d_A, size_A));CHECK_CUDA(cudaMalloc(&d_B, size_B));CHECK_CUDA(cudaMalloc(&d_C, size_C));CHECK_CUDA(cudaMemcpy(d_A, h_A, size_A, cudaMemcpyHostToDevice));CHECK_CUDA(cudaMemcpy(d_B, h_B, size_B, cudaMemcpyHostToDevice));const float alpha = 1.0f;const float beta = 0.0f;// Warmup runCHECK_CUBLAS(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N));CHECK_CUDA(cudaDeviceSynchronize());// Benchmark runsfor (int i = 0; i < num_runs; ++i) {auto start = std::chrono::high_resolution_clock::now();CHECK_CUBLAS(cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, d_B, N, d_A, K, &beta, d_C, N));CHECK_CUDA(cudaDeviceSynchronize());auto end = std::chrono::high_resolution_clock::now();std::chrono::duration<double> diff = end - start;cublas_times.push_back(diff.count());std::cout << "CUBLAS run " << i + 1 << " time: " << diff.count() << " seconds" << std::endl;}CHECK_CUDA(cudaMemcpy(h_C_cublas, d_C, size_C, cudaMemcpyDeviceToHost));CHECK_CUDA(cudaFree(d_A));CHECK_CUDA(cudaFree(d_B));CHECK_CUDA(cudaFree(d_C));CHECK_CUBLAS(cublasDestroy(handle));}// CUBLAS-XT{cublasXtHandle_t handle;CHECK_CUBLAS(cublasXtCreate(&handle));int devices[1] = {0};CHECK_CUBLAS(cublasXtDeviceSelect(handle, 1, devices));const float alpha = 1.0f;const float beta = 0.0f;// Warmup runCHECK_CUBLAS(cublasXtSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, h_B, N, h_A, K, &beta, h_C_cublasxt, N));// Benchmark runsfor (int i = 0; i < num_runs; ++i) {auto start = std::chrono::high_resolution_clock::now();CHECK_CUBLAS(cublasXtSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, M, K, &alpha, h_B, N, h_A, K, &beta,h_C_cublasxt, N));auto end = std::chrono::high_resolution_clock::now();std::chrono::duration<double> diff = end - start;cublasxt_times.push_back(diff.count());std::cout << "CUBLAS-XT run " << i + 1 << " time: " << diff.count() << " seconds" << std::endl;}CHECK_CUBLAS(cublasXtDestroy(handle));}// Calculate and print average timesdouble avg_cublas = 0.0, avg_cublasxt = 0.0;for (int i = 0; i < num_runs; ++i) {avg_cublas += cublas_times[i];avg_cublasxt += cublasxt_times[i];}avg_cublas /= num_runs;avg_cublasxt /= num_runs;std::cout << "Average CUBLAS time: " << avg_cublas << " seconds" << std::endl;std::cout << "Average CUBLAS-XT time: " << avg_cublasxt << " seconds" << std::endl;// Verify resultsfloat tolerance = 1e-4f;bool results_match = compareResults(h_C_cublas, h_C_cublasxt, M * N, tolerance);if (results_match) {std::cout << "Results match within tolerance." << std::endl;} else {std::cout << "Results do not match within tolerance." << std::endl;}free(h_A);free(h_B);free(h_C_cublas);free(h_C_cublasxt);return 0;
}