搭建科技网站价格,视频拍摄教程,京东商家入驻入口官网,wordpress博客示例【CUDA】CUBLAS
在深入了解之前#xff0c;提前运行预热#xff08;warmup#xff09;和基准测试#xff08;benchmark runs#xff09; 是获得准确执行时间的关键。如果不进行预热运行#xff0c;cuBLAS 的首次运行会有较大的开销#xff08;大约 45 毫秒#xff09;…【CUDA】CUBLAS
在深入了解之前提前运行预热warmup和基准测试benchmark runs 是获得准确执行时间的关键。如果不进行预热运行cuBLAS 的首次运行会有较大的开销大约 45 毫秒会导致结果偏差。基准测试能够更准确地计算出平均执行时间。 cuBLASCUDA 基本线性代数子程序库
简介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 中需要单独下载。 CUTLASSCUDA 模板线性代数子程序库
简介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 dimensionconst 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 // 提供对 halfFP16数据类型的支持
#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::vectorhalf h_A_half(M * K);std::vectorhalf 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/*矩阵乘法描述符 cublasLtMatmulDescCreateCreate 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::vectorhalf 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_castfloat(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::vectordouble cublas_times;std::vectordouble 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::durationdouble 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::durationdouble 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;
}参考https://github.com/Infatoshi/cuda-course/tree/master
