2023年的深度學(xué)習(xí)入門指南(10) - CUDA編程基礎(chǔ)

上一篇我們走馬觀花地看了下SIMD和GPGPU的編程。不過線條太粗了，在開發(fā)大模型時(shí)遇到問題了肯定還會(huì)暈。
所以我們還是需要深入到CUDA中去探險(xiǎn)一下。

獲取CUDA設(shè)備信息

在使用CUDA設(shè)備之前，首先我們得獲取是否支持CUDA，有幾個(gè)設(shè)備。這個(gè)可以通過cudaGetDeviceCount

int deviceCount;cudaError_t cudaError;cudaError = cudaGetDeviceCount(&deviceCount);if (cudaError == cudaSuccess) {cout << "There are " << deviceCount << " cuda devices." << endl;}

獲取了支持多少個(gè)設(shè)備了之后，我們就可以遍歷設(shè)備去用cudaGetDeviceProperties函數(shù)去查看設(shè)備信息了。

for (int i = 0; i < deviceCount; i++){cudaError = cudaGetDeviceProperties(&props, i);if (cudaError == cudaSuccess) {cout << "Device Name： " << props.name << endl;cout << "Compute Capability version: " << props.major << "." << props.minor << endl;}}

這是我在我的電腦上輸出的結(jié)果：

There are 1 cuda devices. Device Name： NVIDIA GeForce RTX 3060 Compute Capability version: 8.6 struct cudaDeviceProp {char name[256];cudaUUID_t uuid;size_t totalGlobalMem;size_t sharedMemPerBlock;int regsPerBlock;int warpSize;size_t memPitch;int maxThreadsPerBlock;int maxThreadsDim[3];int maxGridSize[3];int clockRate;size_t totalConstMem;int major;int minor;size_t textureAlignment;size_t texturePitchAlignment;int deviceOverlap;int multiProcessorCount;int kernelExecTimeoutEnabled;int integrated;int canMapHostMemory;int computeMode;int maxTexture1D;int maxTexture1DMipmap;int maxTexture1DLinear;int maxTexture2D[2];int maxTexture2DMipmap[2];int maxTexture2DLinear[3];int maxTexture2DGather[2];int maxTexture3D[3];int maxTexture3DAlt[3];int maxTextureCubemap;int maxTexture1DLayered[2];int maxTexture2DLayered[3];int maxTextureCubemapLayered[2];int maxSurface1D;int maxSurface2D[2];int maxSurface3D[3];int maxSurface1DLayered[2];int maxSurface2DLayered[3];int maxSurfaceCubemap;int maxSurfaceCubemapLayered[2];size_t surfaceAlignment;int concurrentKernels;int ECCEnabled;int pciBusID;int pciDeviceID;int pciDomainID;int tccDriver;int asyncEngineCount;int unifiedAddressing;int memoryClockRate;int memoryBusWidth;int l2CacheSize;int persistingL2CacheMaxSize;int maxThreadsPerMultiProcessor;int streamPrioritiesSupported;int globalL1CacheSupported;int localL1CacheSupported;size_t sharedMemPerMultiprocessor;int regsPerMultiprocessor;int managedMemory;int isMultiGpuBoard;int multiGpuBoardGroupID;int singleToDoublePrecisionPerfRatio;int pageableMemoryAccess;int concurrentManagedAccess;int computePreemptionSupported;int canUseHostPointerForRegisteredMem;int cooperativeLaunch;int cooperativeMultiDeviceLaunch;int pageableMemoryAccessUsesHostPageTables;int directManagedMemAccessFromHost;int accessPolicyMaxWindowSize;}

我們擇其要者介紹幾個(gè)吧：

• totalGlobalMem是設(shè)備上可用的全局內(nèi)存總量，以字節(jié)為單位。
• sharedMemPerBlock是一個(gè)線程塊可用的最大共享內(nèi)存量，以字節(jié)為單位。
• regsPerBlock是一個(gè)線程塊可用的最大32位寄存器數(shù)量。
• warpSize是線程束的大小，以線程為單位。
• memPitch是涉及通過cudaMallocPitch()分配的內(nèi)存區(qū)域的內(nèi)存復(fù)制函數(shù)允許的最大間距，以字節(jié)為單位。
• maxThreadsPerBlock是每個(gè)塊的最大線程數(shù)。
• maxThreadsDim[3]包含了一個(gè)塊的每個(gè)維度的最大尺寸。
• maxGridSize[3]包含了一個(gè)網(wǎng)格的每個(gè)維度的最大尺寸。
• clockRate是時(shí)鐘頻率，以千赫為單位。
• totalConstMem是設(shè)備上可用的常量內(nèi)存總量，以字節(jié)為單位。
• major, minor是定義設(shè)備計(jì)算能力的主要和次要修訂號(hào)。
• multiProcessorCount是設(shè)備上多處理器的數(shù)量。
• memoryClockRate是峰值內(nèi)存時(shí)鐘頻率，以千赫為單位。
• memoryBusWidth是內(nèi)存總線寬度，以位為單位。
• memoryPoolsSupported 是 1，如果設(shè)備支持使用 cudaMallocAsync 和 cudaMemPool 系列 API，否則為 0
• gpuDirectRDMASupported 是 1，如果設(shè)備支持 GPUDirect RDMA API，否則為 0
• gpuDirectRDMAFlushWritesOptions 是一個(gè)按照 cudaFlushGPUDirectRDMAWritesOptions 枚舉解釋的位掩碼
• gpuDirectRDMAWritesOrdering 參見 cudaGPUDirectRDMAWritesOrdering 枚舉的數(shù)值
• memoryPoolSupportedHandleTypes 是一個(gè)支持與 mempool-based IPC 的句柄類型的位掩碼
• deferredMappingCudaArraySupported 是 1，如果設(shè)備支持延遲映射 CUDA 數(shù)組和 CUDA mipmapped 數(shù)組
• ipcEventSupported 是 1，如果設(shè)備支持 IPC 事件，否則為 0
• unifiedFunctionPointers 是 1，如果設(shè)備支持統(tǒng)一指針，否則為 0

有了更多的信息，我們輸出一些看看：

for (int i = 0; i < deviceCount; i++){cudaError = cudaGetDeviceProperties(&props, i);if (cudaError == cudaSuccess) {cout << "Device Name： " << props.name << endl;cout << "Compute Capability version: " << props.major << "." << props.minor << endl;cout << "設(shè)備上可用的全局內(nèi)存總量:(G字節(jié))" << props.totalGlobalMem / 1024 / 1024 / 1024 << endl;cout << "時(shí)鐘頻率（以MHz為單位）:" << props.clockRate / 1000 << endl;cout << "設(shè)備上多處理器的數(shù)量:" << props.multiProcessorCount << endl;cout << "每個(gè)塊的最大線程數(shù):" << props.maxThreadsPerBlock <<endl;cout << "內(nèi)存總線寬度(位)" << props.memoryBusWidth << endl;cout << "一個(gè)塊的每個(gè)維度的最大尺寸:" << props.maxThreadsDim[0] << ","<< props.maxThreadsDim[1] << "," << props.maxThreadsDim[2] << endl;cout << "一個(gè)網(wǎng)格的每個(gè)維度的最大尺寸:" << props.maxGridSize[0] << "," << props.maxGridSize[1] << "," << props.maxGridSize[2] <<endl;}}

在我的3060顯卡上運(yùn)行的結(jié)果：

Device Name： NVIDIA GeForce RTX 3060 Compute Capability version: 8.6 設(shè)備上可用的全局內(nèi)存總量:(G字節(jié))11 時(shí)鐘頻率（以MHz為單位）:1777 設(shè)備上多處理器的數(shù)量:28 每個(gè)塊的最大線程數(shù):1024 內(nèi)存總線寬度(位)192 一個(gè)塊的每個(gè)維度的最大尺寸:1024,1024,64 一個(gè)網(wǎng)格的每個(gè)維度的最大尺寸:2147483647,65535,65535

線程塊和線程網(wǎng)格

在CUDA中，線程塊（block）和線程網(wǎng)格（grid）是兩個(gè)非常重要的概念，它們用于描述GPU執(zhí)行并行任務(wù)時(shí)的線程組織方式。線程塊是由若干個(gè)線程（thread）組成的，它們可以在同一個(gè)GPU多處理器（multiprocessor）上并行執(zhí)行。線程網(wǎng)格則是由若干個(gè)線程塊組成的，它們可以在整個(gè)GPU設(shè)備上并行執(zhí)行。每個(gè)線程塊和線程網(wǎng)格都有一個(gè)唯一的索引，用于在CUDA C/C++的GPU核函數(shù)中對線程進(jìn)行標(biāo)識(shí)和控制。

在CUDA中，使用dim3結(jié)構(gòu)體來表示線程塊和線程網(wǎng)格的維度。例如，dim3(2,2)表示一個(gè)2D線程網(wǎng)格，其中有2x2=4個(gè)線程塊；dim3(2,2,2)表示一個(gè)3D線程塊，其中有2x2x2=8個(gè)線程。在啟動(dòng)GPU核函數(shù)時(shí)，可以使用<<< >>>的語法來指定線程網(wǎng)格和線程塊的大小，例如：

dim3 dimGrid(2, 2); dim3 dimBlock(2, 2, 2); myKernel<<<dimGrid, dimBlock>>>(...);

這里使用dimGrid和dimBlock指定了線程網(wǎng)格和線程塊的大小，然后調(diào)用myKernel函數(shù)，并傳遞必要的參數(shù)。在執(zhí)行GPU核函數(shù)時(shí)，CUDA會(huì)按照指定的線程網(wǎng)格和線程塊的大小啟動(dòng)對應(yīng)的線程，并對它們進(jìn)行分配和協(xié)作，從而完成任務(wù)的并行執(zhí)行。線程塊和線程網(wǎng)格的組織方式和大小都可以根據(jù)具體的應(yīng)用場景和硬件環(huán)境進(jìn)行調(diào)整和優(yōu)化，以實(shí)現(xiàn)最優(yōu)的性能和效率。

我們再看下在核函數(shù)中如何使用線程網(wǎng)格和線程塊。

__global__ void testKernel(int val) {printf("[%d, %d]:\t\tValue is:%d\n", blockIdx.y * gridDim.x + blockIdx.x,threadIdx.z * blockDim.x * blockDim.y + threadIdx.y * blockDim.x +threadIdx.x,val); }

上面有幾個(gè)點(diǎn)我們需要解釋一下：

• __global__：并不是表明這是一個(gè)全局函數(shù)，而是表明這是一個(gè)GPU核函數(shù)。
• blockIdx：是一個(gè)內(nèi)置的變量，表示當(dāng)前線程所在的塊（block）的索引。它是一個(gè)結(jié)構(gòu)體類型，包含了三個(gè)成員變量，分別表示當(dāng)前塊在x、y、z三個(gè)維度上的索引值。
• threadIdx：也是一個(gè)內(nèi)置的變量，表示當(dāng)前線程在所在的塊中的索引。它也同樣是一個(gè)結(jié)構(gòu)體類型，包含了三個(gè)成員變量，分別表示當(dāng)前線程在x、y、z三個(gè)維度上的索引值。
• blockDim：同樣是一個(gè)內(nèi)置的變量，表示每個(gè)塊（block）的維度（dimension），包括x、y、z三個(gè)維度。

在CUDA中，每個(gè)核函數(shù)（kernel function）被分配到一個(gè)或多個(gè)塊（block）中執(zhí)行，每個(gè)塊包含若干個(gè)線程（thread），它們可以在GPU上并行執(zhí)行。通過訪問blockIdx的成員變量，可以確定當(dāng)前線程所在的塊在哪個(gè)位置，從而在核函數(shù)中進(jìn)行特定的計(jì)算。例如，可以使用blockIdx.x表示當(dāng)前線程所在的塊在x軸上的索引值。在CUDA編程中，通常需要使用blockIdx和threadIdx來確定每個(gè)線程在整個(gè)GPU并行執(zhí)行中的唯一標(biāo)識(shí)，以便進(jìn)行任務(wù)的分配和協(xié)作。

然后將dimGrid和dimBlock傳給testKernel.

// Kernel configuration, where a two-dimensional grid and// three-dimensional blocks are configured.dim3 dimGrid(2, 2);dim3 dimBlock(2, 2, 2);testKernel << <dimGrid, dimBlock >> > (10);

將下面的文件保存為kernel.cu，然后通過nvcc命令編譯，最后運(yùn)行生成的可執(zhí)行文件就可以了。

// System includes #include <stdio.h> #include <assert.h> #include <iostream>// CUDA runtime #include <cuda_runtime.h>using namespace std;__global__ void testKernel(int val) {printf("[%d, %d]:\t\tValue is:%d\n", blockIdx.y * gridDim.x + blockIdx.x,threadIdx.z * blockDim.x * blockDim.y + threadIdx.y * blockDim.x +threadIdx.x,val); }int main(int argc, char** argv) {int devID;cudaDeviceProp props;int deviceCount;cudaError_t cudaError;cudaError = cudaGetDeviceCount(&deviceCount);if (cudaError == cudaSuccess) {cout << "There are " << deviceCount << " cuda devices." << endl;}for (int i = 0; i < deviceCount; i++){cudaError = cudaGetDeviceProperties(&props, i);if (cudaError == cudaSuccess) {cout << "Device Name： " << props.name << endl;cout << "Compute Capability version: " << props.major << "." << props.minor << endl;cout << "設(shè)備上可用的全局內(nèi)存總量:(G字節(jié))" << props.totalGlobalMem / 1024 / 1024 / 1024 << endl;cout << "時(shí)鐘頻率（以MHz為單位）:" << props.clockRate / 1000 << endl;cout << "設(shè)備上多處理器的數(shù)量:" << props.multiProcessorCount << endl;cout << "每個(gè)塊的最大線程數(shù):" << props.maxThreadsPerBlock <<endl;cout << "內(nèi)存總線寬度(位)" << props.memoryBusWidth << endl;cout << "一個(gè)塊的每個(gè)維度的最大尺寸:" << props.maxThreadsDim[0] << ","<< props.maxThreadsDim[1] << "," << props.maxThreadsDim[2] << endl;cout << "一個(gè)網(wǎng)格的每個(gè)維度的最大尺寸:" << props.maxGridSize[0] << "," << props.maxGridSize[1] << "," << props.maxGridSize[2] <<endl;}}// Kernel configuration, where a two-dimensional grid and// three-dimensional blocks are configured.dim3 dimGrid(2, 2);dim3 dimBlock(2, 2, 2);testKernel << <dimGrid, dimBlock >> > (10);cudaDeviceSynchronize();return EXIT_SUCCESS; }

前面輸出的不管，我們只看后面32個(gè)線程的結(jié)果：

[1, 0]: Value is:10 [1, 1]: Value is:10 [1, 2]: Value is:10 [1, 3]: Value is:10 [1, 4]: Value is:10 [1, 5]: Value is:10 [1, 6]: Value is:10 [1, 7]: Value is:10 [0, 0]: Value is:10 [0, 1]: Value is:10 [0, 2]: Value is:10 [0, 3]: Value is:10 [0, 4]: Value is:10 [0, 5]: Value is:10 [0, 6]: Value is:10 [0, 7]: Value is:10 [3, 0]: Value is:10 [3, 1]: Value is:10 [3, 2]: Value is:10 [3, 3]: Value is:10 [3, 4]: Value is:10 [3, 5]: Value is:10 [3, 6]: Value is:10 [3, 7]: Value is:10 [2, 0]: Value is:10 [2, 1]: Value is:10 [2, 2]: Value is:10 [2, 3]: Value is:10 [2, 4]: Value is:10 [2, 5]: Value is:10 [2, 6]: Value is:10 [2, 7]: Value is:10

前面表示線程塊，后面表示線程。

大家第一次搞GPU編程的話很容易被繞暈。我來解釋一下這個(gè)計(jì)算方法。其實(shí)就是跟用一維數(shù)組來模擬多維數(shù)組是一個(gè)算法。

blockIdx.y * gridDim.x + blockIdx.x表示當(dāng)前線程所在的線程塊在二維線程網(wǎng)格中的唯一標(biāo)識(shí)。其中，gridDim.x表示線程網(wǎng)格在x方向上的線程塊數(shù)量，blockIdx.x表示當(dāng)前線程塊在x方向上的索引值，blockIdx.y表示當(dāng)前線程塊在y方向上的索引值。

threadIdx.z * blockDim.x * blockDim.y表示當(dāng)前線程在z方向上的偏移量，即前面所有線程所占用的空間大小。然后，threadIdx.y * blockDim.x表示當(dāng)前線程在y方向上的偏移量，即當(dāng)前線程在所在z平面上的偏移量。最后，threadIdx.x表示當(dāng)前線程在x方向上的偏移量，即當(dāng)前線程在所在z平面的某一行上的偏移量。

明白這一點(diǎn)之后，我們嘗試將每個(gè)線程塊從8個(gè)線程改成12個(gè)：

dim3 dimGrid(2, 2);dim3 dimBlock(2, 2, 3);testKernel << <dimGrid, dimBlock >> > (12);

運(yùn)行結(jié)果如下：

[0, 0]: Value is:12 [0, 1]: Value is:12 [0, 2]: Value is:12 [0, 3]: Value is:12 [0, 4]: Value is:12 [0, 5]: Value is:12 [0, 6]: Value is:12 [0, 7]: Value is:12 [0, 8]: Value is:12 [0, 9]: Value is:12 [0, 10]: Value is:12 [0, 11]: Value is:12 [1, 0]: Value is:12 [1, 1]: Value is:12 [1, 2]: Value is:12 [1, 3]: Value is:12 [1, 4]: Value is:12 [1, 5]: Value is:12 [1, 6]: Value is:12 [1, 7]: Value is:12 [1, 8]: Value is:12 [1, 9]: Value is:12 [1, 10]: Value is:12 [1, 11]: Value is:12 [3, 0]: Value is:12 [3, 1]: Value is:12 [3, 2]: Value is:12 [3, 3]: Value is:12 [3, 4]: Value is:12 [3, 5]: Value is:12 [3, 6]: Value is:12 [3, 7]: Value is:12 [3, 8]: Value is:12 [3, 9]: Value is:12 [3, 10]: Value is:12 [3, 11]: Value is:12 [2, 0]: Value is:12 [2, 1]: Value is:12 [2, 2]: Value is:12 [2, 3]: Value is:12 [2, 4]: Value is:12 [2, 5]: Value is:12 [2, 6]: Value is:12 [2, 7]: Value is:12 [2, 8]: Value is:12 [2, 9]: Value is:12 [2, 10]: Value is:12 [2, 11]: Value is:12

下面我們正式開啟真并發(fā)之旅，在上面的48個(gè)線程里同時(shí)計(jì)算正弦。
在GPU里計(jì)算，我們CPU上原來的數(shù)學(xué)庫不頂用了，我們要用GPU自己的，在CUDA中我們用__sinf:

__global__ void testKernel(float val) {printf("[%d, %d]:\t\tValue is:%f\n", blockIdx.y * gridDim.x + blockIdx.x,threadIdx.z * blockDim.x * blockDim.y + threadIdx.y * blockDim.x +threadIdx.x,__sinf(val* threadIdx.z * blockDim.x * blockDim.y + threadIdx.y * blockDim.x +threadIdx.x)); }

main函數(shù)里也隨便改一個(gè)：

dim3 dimGrid(2, 2);dim3 dimBlock(2, 2, 3);testKernel << <dimGrid, dimBlock >> > (0.5);

運(yùn)行結(jié)果如下：

[0, 0]: Value is:0.000000 [0, 1]: Value is:0.841471 [0, 2]: Value is:0.909297 [0, 3]: Value is:0.141120 [0, 4]: Value is:0.909297 [0, 5]: Value is:0.141120 [0, 6]: Value is:-0.756802 [0, 7]: Value is:-0.958924 [0, 8]: Value is:-0.756802 [0, 9]: Value is:-0.958924 [0, 10]: Value is:-0.279416 [0, 11]: Value is:0.656986 [1, 0]: Value is:0.000000 [1, 1]: Value is:0.841471 [1, 2]: Value is:0.909297 [1, 3]: Value is:0.141120 [1, 4]: Value is:0.909297 [1, 5]: Value is:0.141120 [1, 6]: Value is:-0.756802 [1, 7]: Value is:-0.958924 [1, 8]: Value is:-0.756802 [1, 9]: Value is:-0.958924 [1, 10]: Value is:-0.279416 [1, 11]: Value is:0.656986 [3, 0]: Value is:0.000000 [3, 1]: Value is:0.841471 [3, 2]: Value is:0.909297 [3, 3]: Value is:0.141120 [3, 4]: Value is:0.909297 [3, 5]: Value is:0.141120 [3, 6]: Value is:-0.756802 [3, 7]: Value is:-0.958924 [3, 8]: Value is:-0.756802 [3, 9]: Value is:-0.958924 [3, 10]: Value is:-0.279416 [3, 11]: Value is:0.656986 [2, 0]: Value is:0.000000 [2, 1]: Value is:0.841471 [2, 2]: Value is:0.909297 [2, 3]: Value is:0.141120 [2, 4]: Value is:0.909297 [2, 5]: Value is:0.141120 [2, 6]: Value is:-0.756802 [2, 7]: Value is:-0.958924 [2, 8]: Value is:-0.756802 [2, 9]: Value is:-0.958924 [2, 10]: Value is:-0.279416 [2, 11]: Value is:0.656986

內(nèi)存與顯存間的數(shù)據(jù)交換

上面我們是傳了一個(gè)立即數(shù)到GPU核函數(shù)。我們距離正式能使用GPU進(jìn)行CUDA編程，就差分配GPU顯存和在顯存和內(nèi)存之間復(fù)制了。

同malloc類似，CUDA使用cudaMalloc來分配GPU內(nèi)存，其原型為：

cudaError_t cudaMalloc(void **devPtr, size_t size);

參數(shù)解釋:

• devPtr: 返回分配的設(shè)備內(nèi)存的指針。
• size: 要分配的內(nèi)存大小,以字節(jié)為單位。

返回值:

• cudaSuccess: 分配成功。
• cudaErrorInvalidValue: size為零或devPtr為NULL。
• cudaErrorMemoryAllocation: 內(nèi)存分配失敗。

一般的用法，記得用完了用cudaFree釋放掉：

float* devPtr; cudaMalloc(&devPtr, size * sizeof(float)); ... cudaFree(devPtr);

分配完內(nèi)存了，然后就是從內(nèi)存復(fù)制到顯存了。同樣類似于memcpy，通過cudaMemcpy來完成。

cudaError_t cudaMemcpy(void* dst, const void* src, size_t count, cudaMemcpyKind kind);

參數(shù)解釋:

• dst: 目標(biāo)內(nèi)存的指針。
• src: 源內(nèi)存的指針。
• count: 要拷貝的內(nèi)存大小,以字節(jié)為單位。
• kind: 拷貝的類型,可以是:
  • cudaMemcpyHostToHost
  • cudaMemcpyHostToDevice
  • cudaMemcpyDeviceToHost
  • cudaMemcpyDeviceToDevice

返回值:

• cudaSuccess: 拷貝成功。
• cudaErrorInvalidValue: count或dst或src為NULL。
• cudaErrorMemoryAllocation: 內(nèi)存分配失敗。

下面我們來寫一個(gè)用CUDA計(jì)算平方根的例子：

const int n = 1024;size_t size = n * sizeof(float);float* h_in = (float*)malloc(size);float* h_out = (float*)malloc(size);float* d_in, * d_out;// Initialize input arrayfor (int i = 0; i < n; ++i) {h_in[i] = (float)i;}// Allocate device memorycudaMalloc(&d_in, size);cudaMalloc(&d_out, size);// Copy input data to devicecudaMemcpy(d_in, h_in, size, cudaMemcpyHostToDevice);// Launch kernelint threadsPerBlock = 256;int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;sqrtKernel << <blocksPerGrid, threadsPerBlock >> > (d_in, d_out, n);// Copy output data to hostcudaMemcpy(h_out, d_out, size, cudaMemcpyDeviceToHost);// Verify resultsfor (int i = 0; i < n; ++i) {if (fabsf(h_out[i] - sqrtf(h_in[i])) > 1e-5) {printf("Error: h_out[%d] = %f, sqrtf(h_in[%d]) = %f\n", i, h_out[i], i, sqrtf(h_in[i]));}}printf("Success!\n");// Free memoryfree(h_in);free(h_out);cudaFree(d_in);cudaFree(d_out);

大家關(guān)注線程塊數(shù)和線程數(shù)這兩個(gè)，我們這里沒有用多維，就是用兩個(gè)整數(shù)計(jì)算的：

int threadsPerBlock = 256;int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;sqrtKernel << <blocksPerGrid, threadsPerBlock >> > (d_in, d_out, n);

我們用4個(gè)塊，每個(gè)塊有256個(gè)線程。

此時(shí)，就不用計(jì)算y和z了，只計(jì)算x維度就可以：

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

但是要注意，blockIdx和threadIdx仍然是三維的，y和z維仍然是有效的，只不過它們變成0了。

我們的核函數(shù)這樣寫：

__global__ void sqrtKernel(float* in, float* out, int n) {int i = blockIdx.x * blockDim.x + threadIdx.x;if (i < n) {out[i] = sqrtf(in[i]);printf("[%d, %d]:\t\tValue is:%f\n", blockIdx.y * gridDim.x + blockIdx.x,threadIdx.z * blockDim.x * blockDim.y + threadIdx.y * blockDim.x +threadIdx.x, out[i]);} }

當(dāng)然了，因?yàn)閎lock和thread的y和z都是0，跟只寫x是沒啥區(qū)別的：

__global__ void sqrtKernel(float* in, float* out, int n) {int i = blockIdx.x * blockDim.x + threadIdx.x;if (i < n) {out[i] = sqrtf(in[i]);printf("[%d, %d]:\t\tValue is:%f\n", blockIdx.x, threadIdx.x, out[i]);} }

使用封裝好的庫

除了CUDA運(yùn)行時(shí)之外，針對主要的應(yīng)用場景，NVidia也提供了很多專門的庫。

比如針對矩陣運(yùn)算，就有cuBLAS庫。有的庫是跟隨CUDA工具包一起安裝的，比如cuBLAS, cuFFT。也有的庫需要專門下載安裝，比如cudnn庫。

這里強(qiáng)調(diào)一下，所謂的庫，不是在核函數(shù)中要調(diào)用的模塊，而是將整個(gè)需要在核函數(shù)里面要實(shí)現(xiàn)的功能全封裝好了。所以在使用封裝庫的時(shí)候，并不需要nvcc，就是引用一個(gè)庫就好了。

我們來看一個(gè)使用cuBLAS庫來計(jì)算矩陣乘法的例子。

cuBLAS庫來計(jì)算矩陣乘法要用到的主要的函數(shù)有4個(gè)：

• cublasCreate: 創(chuàng)建cublas句柄
• cublasDestroy：釋放cublas句柄
• cublasSetVector：在CPU和GPU內(nèi)存間復(fù)制數(shù)據(jù)
• cublasSgemm：矩陣乘法運(yùn)算

cublasStatus_t cublasSetVector(int n, int elemSize, const void *x, int incx, void *y, int incy)

其中:

• n 是要拷貝的元素個(gè)數(shù)
• elemSize是每個(gè)元素的大小（以字節(jié)為單位）
• x是主機(jī)端（CPU）內(nèi)存中的數(shù)據(jù)起始地址
• incx是x中相鄰元素之間的跨度
• y是GPU設(shè)備內(nèi)存中的數(shù)據(jù)起始地址
• incy是y中相鄰元素之間的跨度

cublasStatus_t cublasSgemm(cublasHandle_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)

其中:

• handle是cuBLAS句柄；
• transa是A矩陣的轉(zhuǎn)置選項(xiàng)，取值為CUBLAS_OP_N或CUBLAS_OP_T，分別表示不轉(zhuǎn)置和轉(zhuǎn)置；
• transb是B矩陣的轉(zhuǎn)置選項(xiàng)；m、n、k分別是A、B、C矩陣的維度；
• alpha是一個(gè)標(biāo)量值，用于將A和B矩陣的乘積縮放到C矩陣中；
• A是A矩陣的起始地址；
• lda是A矩陣中相鄰列之間的跨度；
• B是B矩陣的起始地址；
• ldb是B矩陣中相鄰列之間的跨度；
• beta是一個(gè)標(biāo)量值，用于將C矩陣中的值縮放；
• C是C矩陣的起始地址；
• ldc是C矩陣中相鄰列之間的跨度。

我們簡化寫一個(gè)例子，主要說明函數(shù)的用法：

#include <stdio.h> #include <cuda_runtime.h> #include <cublas_v2.h>int main() {int m = 1024, n = 1024, k = 1024;float* h_A = (float*)malloc(m * k * sizeof(float));float* h_B = (float*)malloc(k * n * sizeof(float));float* h_C = (float*)malloc(m * n * sizeof(float));for (int i = 0; i < m * k; ++i) {h_A[i] = (float)i;}for (int i = 0; i < k * n; ++i) {h_B[i] = (float)i;}float* d_A, * d_B, * d_C;cudaMalloc(&d_A, m * k * sizeof(float));cudaMalloc(&d_B, k * n * sizeof(float));cudaMalloc(&d_C, m * n * sizeof(float));// Copy data from host to devicecublasSetVector(m * k, sizeof(float), h_A, 1, d_A, 1);cublasSetVector(k * n, sizeof(float), h_B, 1, d_B, 1);// Initialize cuBLAScublasHandle_t handle;cublasCreate(&handle);// Do matrix multiplicationconst float alpha = 1.0f, beta = 0.0f;cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, m, n, k,&alpha, d_A, m, d_B, k, &beta, d_C, m);// Copy data from device to hostcublasGetVector(m * n, sizeof(float), d_C, 1, h_C, 1);// Free memoryfree(h_A);free(h_B);free(h_C);cudaFree(d_A);cudaFree(d_B);cudaFree(d_C);// Destroy cuBLAS handlecublasDestroy(handle);return 0; }

當(dāng)然，上面的只是個(gè)例子，沒有做錯(cuò)誤處理，這樣是不對的。
我們參考官方的例子：

#include <stdio.h> #include <stdlib.h> #include <string.h>/* Includes, cuda */ #include <cublas_v2.h> #include <cuda_runtime.h> #include <helper_cuda.h>/* Matrix size */ #define N (275)/* Host implementation of a simple version of sgemm */ static void simple_sgemm(int n, float alpha, const float *A, const float *B,float beta, float *C) {int i;int j;int k;for (i = 0; i < n; ++i) {for (j = 0; j < n; ++j) {float prod = 0;for (k = 0; k < n; ++k) {prod += A[k * n + i] * B[j * n + k];}C[j * n + i] = alpha * prod + beta * C[j * n + i];}} }/* Main */ int main(int argc, char **argv) {cublasStatus_t status;float *h_A;float *h_B;float *h_C;float *h_C_ref;float *d_A = 0;float *d_B = 0;float *d_C = 0;float alpha = 1.0f;float beta = 0.0f;int n2 = N * N;int i;float error_norm;float ref_norm;float diff;cublasHandle_t handle;/* Initialize CUBLAS */printf("simpleCUBLAS test running..\n");status = cublasCreate(&handle);if (status != CUBLAS_STATUS_SUCCESS) {fprintf(stderr, "!!!! CUBLAS initialization error\n");return EXIT_FAILURE;}/* Allocate host memory for the matrices */h_A = reinterpret_cast<float *>(malloc(n2 * sizeof(h_A[0])));if (h_A == 0) {fprintf(stderr, "!!!! host memory allocation error (A)\n");return EXIT_FAILURE;}h_B = reinterpret_cast<float *>(malloc(n2 * sizeof(h_B[0])));if (h_B == 0) {fprintf(stderr, "!!!! host memory allocation error (B)\n");return EXIT_FAILURE;}h_C = reinterpret_cast<float *>(malloc(n2 * sizeof(h_C[0])));if (h_C == 0) {fprintf(stderr, "!!!! host memory allocation error (C)\n");return EXIT_FAILURE;}/* Fill the matrices with test data */for (i = 0; i < n2; i++) {h_A[i] = rand() / static_cast<float>(RAND_MAX);h_B[i] = rand() / static_cast<float>(RAND_MAX);h_C[i] = rand() / static_cast<float>(RAND_MAX);}/* Allocate device memory for the matrices */if (cudaMalloc(reinterpret_cast<void **>(&d_A), n2 * sizeof(d_A[0])) !=cudaSuccess) {fprintf(stderr, "!!!! device memory allocation error (allocate A)\n");return EXIT_FAILURE;}if (cudaMalloc(reinterpret_cast<void **>(&d_B), n2 * sizeof(d_B[0])) !=cudaSuccess) {fprintf(stderr, "!!!! device memory allocation error (allocate B)\n");return EXIT_FAILURE;}if (cudaMalloc(reinterpret_cast<void **>(&d_C), n2 * sizeof(d_C[0])) !=cudaSuccess) {fprintf(stderr, "!!!! device memory allocation error (allocate C)\n");return EXIT_FAILURE;}/* Initialize the device matrices with the host matrices */status = cublasSetVector(n2, sizeof(h_A[0]), h_A, 1, d_A, 1);if (status != CUBLAS_STATUS_SUCCESS) {fprintf(stderr, "!!!! device access error (write A)\n");return EXIT_FAILURE;}status = cublasSetVector(n2, sizeof(h_B[0]), h_B, 1, d_B, 1);if (status != CUBLAS_STATUS_SUCCESS) {fprintf(stderr, "!!!! device access error (write B)\n");return EXIT_FAILURE;}status = cublasSetVector(n2, sizeof(h_C[0]), h_C, 1, d_C, 1);if (status != CUBLAS_STATUS_SUCCESS) {fprintf(stderr, "!!!! device access error (write C)\n");return EXIT_FAILURE;}/* Performs operation using plain C code */simple_sgemm(N, alpha, h_A, h_B, beta, h_C);h_C_ref = h_C;/* Performs operation using cublas */status = cublasSgemm(handle, CUBLAS_OP_N, CUBLAS_OP_N, N, N, N, &alpha, d_A,N, d_B, N, &beta, d_C, N);if (status != CUBLAS_STATUS_SUCCESS) {fprintf(stderr, "!!!! kernel execution error.\n");return EXIT_FAILURE;}/* Allocate host memory for reading back the result from device memory */h_C = reinterpret_cast<float *>(malloc(n2 * sizeof(h_C[0])));if (h_C == 0) {fprintf(stderr, "!!!! host memory allocation error (C)\n");return EXIT_FAILURE;}/* Read the result back */status = cublasGetVector(n2, sizeof(h_C[0]), d_C, 1, h_C, 1);if (status != CUBLAS_STATUS_SUCCESS) {fprintf(stderr, "!!!! device access error (read C)\n");return EXIT_FAILURE;}/* Check result against reference */error_norm = 0;ref_norm = 0;for (i = 0; i < n2; ++i) {diff = h_C_ref[i] - h_C[i];error_norm += diff * diff;ref_norm += h_C_ref[i] * h_C_ref[i];}error_norm = static_cast<float>(sqrt(static_cast<double>(error_norm)));ref_norm = static_cast<float>(sqrt(static_cast<double>(ref_norm)));if (fabs(ref_norm) < 1e-7) {fprintf(stderr, "!!!! reference norm is 0\n");return EXIT_FAILURE;}/* Memory clean up */free(h_A);free(h_B);free(h_C);free(h_C_ref);if (cudaFree(d_A) != cudaSuccess) {fprintf(stderr, "!!!! memory free error (A)\n");return EXIT_FAILURE;}if (cudaFree(d_B) != cudaSuccess) {fprintf(stderr, "!!!! memory free error (B)\n");return EXIT_FAILURE;}if (cudaFree(d_C) != cudaSuccess) {fprintf(stderr, "!!!! memory free error (C)\n");return EXIT_FAILURE;}/* Shutdown */status = cublasDestroy(handle);if (status != CUBLAS_STATUS_SUCCESS) {fprintf(stderr, "!!!! shutdown error (A)\n");return EXIT_FAILURE;}if (error_norm / ref_norm < 1e-6f) {printf("simpleCUBLAS test passed.\n");exit(EXIT_SUCCESS);} else {printf("simpleCUBLAS test failed.\n");exit(EXIT_FAILURE);} }

一些更高級(jí)的特性

有了上面的基礎(chǔ)，我們就可以寫一些可以運(yùn)行在GPU上的代碼了。

結(jié)束之前，我們再看幾個(gè)稍微高級(jí)一點(diǎn)的特性。

__device__關(guān)鍵字

之前我們學(xué)習(xí)核函數(shù)的__global__關(guān)鍵字。核函數(shù)既可以被CPU端調(diào)用，也可以被GPU調(diào)用。

如果我們想編寫只能在GPU上運(yùn)行的函數(shù)，我們就可以使用__device__.

使用__device__定義的函數(shù)或變量只能在設(shè)備代碼中使用，無法在主機(jī)端代碼中使用。在CUDA程序中，通常使用__host__和__device__關(guān)鍵字來指定函數(shù)或變量在主機(jī)端和設(shè)備端的執(zhí)行位置。使用__device__定義的函數(shù)或變量可以在設(shè)備代碼中被其他函數(shù)調(diào)用，也可以在主機(jī)端使用CUDA API將數(shù)據(jù)從主機(jī)內(nèi)存?zhèn)鬏數(shù)皆O(shè)備內(nèi)存后，由設(shè)備上的函數(shù)處理。

GPU函數(shù)的內(nèi)聯(lián)

與CPU函數(shù)一樣，GPU上的函數(shù)也可以內(nèi)聯(lián)，使用__forceinline__關(guān)鍵字。

并發(fā)的"?:"三目運(yùn)算符

在C語言中，"?:"三目運(yùn)算符只能做一次判斷。
現(xiàn)在來到了GPU的世界，并發(fā)能力變強(qiáng)了，可以做多次判斷了。

我們來看個(gè)例子：

__device__ __forceinline__ int qcompare(unsigned &val1, unsigned &val2) {return (val1 > val2) ? 1 : (val1 == val2) ? 0 : -1; }

PTX匯編

在上一篇我們學(xué)習(xí)SIMD指令的時(shí)候，我們基本都要內(nèi)聯(lián)匯編。那么在CUDA里面是不是有匯編呢？
答案是肯定的，既然要做性能優(yōu)化，那么肯定要挖掘一切潛力。
不過，為了避免跟架構(gòu)過于相關(guān)，NVidia給我們提供的是一種中間指令格式PTX（Parallel Thread Execution）。
PTX assembly是CUDA的一種中間匯編語言，它是一種與機(jī)器無關(guān)的指令集架構(gòu)（ISA），用于描述GPU上的并行線程執(zhí)行。PTX assembly可以被編譯成特定GPU家族的實(shí)際執(zhí)行的機(jī)器碼。使用PTX assembly可以實(shí)現(xiàn)跨GPU的兼容性和性能優(yōu)化。

我們來看一段內(nèi)嵌匯編：

static __device__ __forceinline__ unsigned int __qsflo(unsigned int word) {unsigned int ret;asm volatile("bfind.u32 %0, %1;" : "=r"(ret) : "r"(word));return ret; }

其中用到的bfind.u32指令用于查找一個(gè)無符號(hào)整數(shù)中最右邊的非零位（即最低有效位），并返回其位位置。該指令將無符號(hào)整數(shù)作為操作數(shù)輸入，并將最低有效位的位位置輸出到目的操作數(shù)中。
“=r”(ret)表示輸出寄存器,返回結(jié)果保存在ret中。
“r”(word)表示輸入寄存器,將參數(shù)word作為輸入。

GPU特有的算法

最后一點(diǎn)要強(qiáng)調(diào)的時(shí)，很多時(shí)候?qū)⒋a并行化，并不是簡簡單單的從CPU轉(zhuǎn)到GPU，而很有可能是要改變算法。

比如，quicksort是一個(gè)(nlog(n))的算法，而bitonic sort是個(gè)$(nlo{g}^2(n))$的算法。但是，bitonic sort更適合于在GPU加速。所以我們在CPU上的quicksort改成bitonic sort算法會(huì)更好一些。

小結(jié)

在Intel CPU還是8+4核20線程的時(shí)候，GTX 1060顯卡做到1280個(gè)CUDA核，3060是3584個(gè)CUDA核，3090是10496個(gè)CUDA核，4090有16384個(gè)CUDA核。每個(gè)CUDA核上可以起比如1024個(gè)線程。
所以，如果有大量可以并發(fā)的任務(wù)，應(yīng)該毫不猶豫地將其寫成核函數(shù)放到GPU上去運(yùn)行。

GPU編程既沒有那么復(fù)雜，完全可以快速上手像寫CPU程序一樣去寫。但也不是那么簡單，適合GPU可能需要改用特殊的算法。

而基于大量簡單Transformers組成的大模型，恰恰是適合高并發(fā)的計(jì)算。

總結(jié)

以上是生活随笔為你收集整理的2023年的深度学习入门指南(10) - CUDA编程基础的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯(cuò)，歡迎將生活随笔推薦給好友。