2023年的深度学习入门指南(10) – CUDA编程根底

上一篇咱们蜻蜓点水地看了下SIMD和GPGPU的编程。不过线条太粗了，在开发大模型时遇到问题了必定还会晕。 所以咱们还是需求深化到CUDA中去探险一下。

获取CUDA设备信息

在运用CUDA设备之前，首先咱们得获取是否支撑CUDA，有几个设备。这个能够经过cudaGetDeviceCount

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

获取了支撑多少个设备了之后，咱们就能够遍历设备去用cudaGetDeviceProperties函数去查看设备信息了。

    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;
        }
    }

这是我在我的电脑上输出的成果：

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;
          }

咱们择其要者介绍几个吧：

• totalGlobalMem是设备上可用的大局内存总量，以字节为单位。
• sharedMemPerBlock是一个线程块可用的最大同享内存量，以字节为单位。
• regsPerBlock是一个线程块可用的最大32位寄存器数量。
• warpSize是线程束的巨细，以线程为单位。
• memPitch是涉及经过cudaMallocPitch()分配的内存区域的内存仿制函数允许的最大间距，以字节为单位。
• maxThreadsPerBlock是每个块的最大线程数。
• maxThreadsDim[3]包括了一个块的每个维度的最大尺度。
• maxGridSize[3]包括了一个网格的每个维度的最大尺度。
• clockRate是时钟频率，以千赫为单位。
• totalConstMem是设备上可用的常量内存总量，以字节为单位。
• major, minor是界说设备核算才能的首要和次要修订号。
• multiProcessorCount是设备上多处理器的数量。
• memoryClockRate是峰值内存时钟频率，以千赫为单位。
• memoryBusWidth是内存总线宽度，以位为单位。
• memoryPoolsSupported 是 1，假如设备支撑运用 cudaMallocAsync 和 cudaMemPool 系列 API，否则为 0
• gpuDirectRDMASupported 是 1，假如设备支撑 GPUDirect RDMA API，否则为 0
• gpuDirectRDMAFlushWritesOptions 是一个按照 cudaFlushGPUDirectRDMAWritesOptions 枚举解说的位掩码
• gpuDirectRDMAWritesOrdering 参见 cudaGPUDirectRDMAWritesOrdering 枚举的数值
• memoryPoolSupportedHandleTypes 是一个支撑与 mempool-based IPC 的句柄类型的位掩码
• deferredMappingCudaArraySupported 是 1，假如设备支撑延迟映射 CUDA 数组和 CUDA mipmapped 数组
• ipcEventSupported 是 1，假如设备支撑 IPC 事件，否则为 0
• unifiedFunctionPointers 是 1，假如设备支撑统一指针，否则为 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 << "设备上可用的大局内存总量:(G字节)" << props.totalGlobalMem / 1024 / 1024 / 1024 << endl;
            cout << "时钟频率（以MHz为单位）:" << props.clockRate / 1000 << endl;
            cout << "设备上多处理器的数量:" << props.multiProcessorCount << endl;
            cout << "每个块的最大线程数:" << props.maxThreadsPerBlock <<endl;
            cout << "内存总线宽度(位)" << props.memoryBusWidth << endl;
            cout << "一个块的每个维度的最大尺度:" << props.maxThreadsDim[0] << ","<< props.maxThreadsDim[1] << "," << props.maxThreadsDim[2] << endl;
            cout << "一个网格的每个维度的最大尺度:" << props.maxGridSize[0] << "," << props.maxGridSize[1] << "," << props.maxGridSize[2] <<endl;
        }
    }

在我的3060显卡上运行的成果：

Device Name： NVIDIA GeForce RTX 3060
Compute Capability version: 8.6
设备上可用的大局内存总量:(G字节)11
时钟频率（以MHz为单位）:1777
设备上多处理器的数量:28
每个块的最大线程数:1024
内存总线宽度(位)192
一个块的每个维度的最大尺度:1024,1024,64
一个网格的每个维度的最大尺度:2147483647,65535,65535

线程块和线程网格

在CUDA中，线程块（block）和线程网格（grid）是两个非常重要的概念，它们用于描述GPU履行并行使命时的线程组织方法。线程块是由若干个线程（thread）组成的，它们能够在同一个GPU多处理器（multiprocessor）上并行履行。线程网格则是由若干个线程块组成的，它们能够在整个GPU设备上并行履行。每个线程块和线程网格都有一个仅有的索引，用于在CUDA C/C++的GPU核函数中对线程进行标识和控制。

在CUDA中，运用dim3结构体来标明线程块和线程网格的维度。例如，dim3(2,2)标明一个2D线程网格，其间有2×2=4个线程块；dim3(2,2,2)标明一个3D线程块，其间有2x2x2=8个线程。在发动GPU核函数时，能够运用<<< >>>的语法来指定线程网格和线程块的巨细，例如：

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

这里运用dimGrid和dimBlock指定了线程网格和线程块的巨细，然后调用myKernel函数，并传递必要的参数。在履行GPU核函数时，CUDA会按照指定的线程网格和线程块的巨细发动对应的线程，并对它们进行分配和协作，从而完成使命的并行履行。线程块和线程网格的组织方法和巨细都能够依据具体的运用场景和硬件环境进行调整和优化，以完成最优的性能和功率。

咱们再看下在核函数中怎么运用线程网格和线程块。

__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);
}

上面有几个点咱们需求解说一下：

• __global__：并不是标明这是一个大局函数，而是标明这是一个GPU核函数。
• blockIdx：是一个内置的变量，标明当时线程地点的块（block）的索引。它是一个结构体类型，包括了三个成员变量，分别标明当时块在x、y、z三个维度上的索引值。
• threadIdx：也是一个内置的变量，标明当时线程在地点的块中的索引。它也相同是一个结构体类型，包括了三个成员变量，分别标明当时线程在x、y、z三个维度上的索引值。
• blockDim：相同是一个内置的变量，标明每个块（block）的维度（dimension），包括x、y、z三个维度。

在CUDA中，每个核函数（kernel function）被分配到一个或多个块（block）中履行，每个块包括若干个线程（thread），它们能够在GPU上并行履行。经过拜访blockIdx的成员变量，能够确定当时线程地点的块在哪个方位，从而在核函数中进行特定的核算。例如，能够运用blockIdx.x标明当时线程地点的块在x轴上的索引值。在CUDA编程中，一般需求运用blockIdx和threadIdx来确定每个线程在整个GPU并行履行中的仅有标识，以便进行使命的分配和协作。

然后将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指令编译，最终运行生成的可履行文件就能够了。

// 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 << "设备上可用的大局内存总量:(G字节)" << props.totalGlobalMem / 1024 / 1024 / 1024 << endl;
            cout << "时钟频率（以MHz为单位）:" << props.clockRate / 1000 << endl;
            cout << "设备上多处理器的数量:" << props.multiProcessorCount << endl;
            cout << "每个块的最大线程数:" << props.maxThreadsPerBlock <<endl;
            cout << "内存总线宽度(位)" << props.memoryBusWidth << endl;
            cout << "一个块的每个维度的最大尺度:" << props.maxThreadsDim[0] << ","<< props.maxThreadsDim[1] << "," << props.maxThreadsDim[2] << endl;
            cout << "一个网格的每个维度的最大尺度:" << 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个线程的成果：

[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编程的话很简单被绕晕。我来解说一下这个核算方法。其实便是跟用一维数组来模仿多维数组是一个算法。

blockIdx.y * gridDim.x + blockIdx.x标明当时线程地点的线程块在二维线程网格中的仅有标识。其间，gridDim.x标明线程网格在x方向上的线程块数量，blockIdx.x标明当时线程块在x方向上的索引值，blockIdx.y标明当时线程块在y方向上的索引值。

threadIdx.z * blockDim.x * blockDim.y标明当时线程在z方向上的偏移量，即前面全部线程所占用的空间巨细。然后，threadIdx.y * blockDim.x标明当时线程在y方向上的偏移量，即当时线程在地点z平面上的偏移量。最终，threadIdx.x标明当时线程在x方向上的偏移量，即当时线程在地点z平面的某一行上的偏移量。

明白这一点之后，咱们尝试将每个线程块从8个线程改成12个：

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

运行成果如下：

[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

下面咱们正式开启真并发之旅，在上面的48个线程里同时核算正弦。 在GPU里核算，咱们CPU上本来的数学库不顶用了，咱们要用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函数里也随意改一个：

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

运行成果如下：

[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

内存与显存间的数据交换

上面咱们是传了一个当即数到GPU核函数。咱们距离正式能运用GPU进行CUDA编程，就差分配GPU显存和在显存和内存之间仿制了。

同malloc相似，CUDA运用cudaMalloc来分配GPU内存，其原型为：

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

参数解说:

• devPtr: 回来分配的设备内存的指针。
• size: 要分配的内存巨细,以字节为单位。

回来值:

• cudaSuccess: 分配成功。
• cudaErrorInvalidValue: size为零或devPtr为NULL。
• cudaErrorMemoryAllocation: 内存分配失利。

一般的用法，记住用完了用cudaFree开释掉：

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

分配完内存了，然后便是从内存仿制到显存了。相同相似于memcpy，经过cudaMemcpy来完成。

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

参数解说:

• dst: 目标内存的指针。
• src: 源内存的指针。
• count: 要仿制的内存巨细,以字节为单位。
• kind: 仿制的类型,能够是:
  • cudaMemcpyHostToHost
  • cudaMemcpyHostToDevice
  • cudaMemcpyDeviceToHost
  • cudaMemcpyDeviceToDevice

回来值:

• cudaSuccess: 仿制成功。
• cudaErrorInvalidValue: count或dst或src为NULL。
• cudaErrorMemoryAllocation: 内存分配失利。

下面咱们来写一个用CUDA核算平方根的比如：

    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 array
    for (int i = 0; i < n; ++i) {
        h_in[i] = (float)i;
    }
    // Allocate device memory
    cudaMalloc(&d_in, size);
    cudaMalloc(&d_out, size);
    // Copy input data to device
    cudaMemcpy(d_in, h_in, size, cudaMemcpyHostToDevice);
    // Launch kernel
    int threadsPerBlock = 256;
    int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
    sqrtKernel << <blocksPerGrid, threadsPerBlock >> > (d_in, d_out, n);
    // Copy output data to host
    cudaMemcpy(h_out, d_out, size, cudaMemcpyDeviceToHost);
    // Verify results
    for (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 memory
    free(h_in);
    free(h_out);
    cudaFree(d_in);
    cudaFree(d_out);

咱们重视线程块数和线程数这两个，咱们这里没有用多维，便是用两个整数核算的：

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

咱们用4个块，每个块有256个线程。

此时，就不用核算y和z了，只核算x维度就能够：

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

但是要注意，blockIdx和threadIdx仍然是三维的，y和z维仍然是有用的，只不过它们变成0了。

咱们的核函数这样写：

__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]);
    }
}

当然了，由于block和thread的y和z都是0，跟只写x是没啥区别的：

__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运行时之外，针对首要的运用场景，NVidia也供给了很多专门的库。

比方针对矩阵运算，就有cuBLAS库。有的库是跟从CUDA工具包一同安装的，比方cuBLAS, cuFFT。也有的库需求专门下载安装，比方cudnn库。

这里着重一下，所谓的库，不是在核函数中要调用的模块，而是将整个需求在核函数里边要完成的功用全封装好了。所以在运用封装库的时分，并不需求nvcc，便是引证一个库就好了。

咱们来看一个运用cuBLAS库来核算矩阵乘法的比如。

cuBLAS库来核算矩阵乘法要用到的首要的函数有4个：

• cublasCreate: 创立cublas句柄
• cublasDestroy：开释cublas句柄
• cublasSetVector：在CPU和GPU内存间仿制数据
• cublasSgemm：矩阵乘法运算

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

其间:

• n 是要仿制的元素个数
• elemSize是每个元素的巨细（以字节为单位）
• x是主机端（CPU）内存中的数据开始地址
• incx是x中相邻元素之间的跨度
• y是GPU设备内存中的数据开始地址
• 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矩阵的转置选项，取值为CUBLAS_OP_N或CUBLAS_OP_T，分别标明不转置和转置；
• transb是B矩阵的转置选项；m、n、k分别是A、B、C矩阵的维度；
• alpha是一个标量值，用于将A和B矩阵的乘积缩放到C矩阵中；
• A是A矩阵的开始地址；
• lda是A矩阵中相邻列之间的跨度；
• B是B矩阵的开始地址；
• ldb是B矩阵中相邻列之间的跨度；
• beta是一个标量值，用于将C矩阵中的值缩放；
• C是C矩阵的开始地址；
• ldc是C矩阵中相邻列之间的跨度。

咱们简化写一个比如，首要说明函数的用法：

#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 device
    cublasSetVector(m * k, sizeof(float), h_A, 1, d_A, 1);
    cublasSetVector(k * n, sizeof(float), h_B, 1, d_B, 1);
    // Initialize cuBLAS
    cublasHandle_t handle;
    cublasCreate(&handle);
    // Do matrix multiplication
    const 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 host
    cublasGetVector(m * n, sizeof(float), d_C, 1, h_C, 1);
    // Free memory
    free(h_A);
    free(h_B);
    free(h_C);
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    // Destroy cuBLAS handle
    cublasDestroy(handle);
    return 0;
}

当然，上面的只是个比如，没有做错误处理，这样是不对的。 咱们参阅官方的比如：

#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);
  }
}

一些更高档的特性

有了上面的根底，咱们就能够写一些能够运行在GPU上的代码了。

完毕之前，咱们再看几个稍微高档一点的特性。

__device__关键字

之前咱们学习核函数的__global__关键字。核函数既能够被CPU端调用，也能够被GPU调用。

假如咱们想编写只能在GPU上运行的函数，咱们就能够运用__device__.

运用__device__界说的函数或变量只能在设备代码中运用，无法在主机端代码中运用。在CUDA程序中，一般运用__host__和__device__关键字来指定函数或变量在主机端和设备端的履行方位。运用__device__界说的函数或变量能够在设备代码中被其他函数调用，也能够在主机端运用CUDA API将数据从主机内存传输到设备内存后，由设备上的函数处理。

GPU函数的内联

与CPU函数相同，GPU上的函数也能够内联，运用__forceinline__关键字。

并发的”?:”三目运算符

在C语言中，”?:”三目运算符只能做一次判别。 现在来到了GPU的世界，并发才能变强了，能够做屡次判别了。

咱们来看个比如：

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

PTX汇编

在上一篇咱们学习SIMD指令的时分，咱们根本都要内联汇编。那么在CUDA里边是不是有汇编呢？ 答案是必定的，已然要做性能优化，那么必定要发掘全部潜力。 不过，为了避免跟架构过于相关，NVidia给咱们供给的是一种中心指令格式PTX（Parallel Thread Execution）。 PTX assembly是CUDA的一种中心汇编语言，它是一种与机器无关的指令集架构（ISA），用于描述GPU上的并行线程履行。PTX assembly能够被编译成特定GPU宗族的实际履行的机器码。运用PTX assembly能够完成跨GPU的兼容性和性能优化。

咱们来看一段内嵌汇编：

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指令用于查找一个无符号整数中最右边的非零位（即最低有用位），并回来其位方位。该指令将无符号整数作为操作数输入，并将最低有用位的位方位输出到目的操作数中。 “=r”(ret)标明输出寄存器,回来成果保存在ret中。 “r”(word)标明输入寄存器,将参数word作为输入。

GPU特有的算法

最终一点要着重的时，很多时分将代码并行化，并不是简简单单的从CPU转到GPU，而很有或许是要改变算法。

比方，quicksort是一个(nlog(n))的算法，而bitonic sort是个(nlog2(n))(nlog^2(n))(nlog2(n))的算法。但是，bitonic sort更合适于在GPU加速。所以咱们在CPU上的quicksort改成bitonic sort算法会更好一些。

小结

在Intel CPU还是8+4核20线程的时分，GTX 1060显卡做到1280个CUDA核，3060是3584个CUDA核，3090是10496个CUDA核，4090有16384个CUDA核。每个CUDA核上能够起比方1024个线程。 所以，假如有很多能够并发的使命，应该毫不犹豫地将其写成核函数放到GPU上去运行。

GPU编程既没有那么复杂，完全能够快速上手像写CPU程序相同去写。但也不是那么简单，合适GPU或许需求改用特殊的算法。

而根据很多简单Transformers组成的大模型，恰恰是合适高并发的核算。