最近新学完cuda编程基础，有些知识点比较晦涩，理解起来不太容易。有必要在学完cuda基础后做一次梳理，加深对cuda编程基础知识的理解。涉及的知识点包括cuda模型、cuda内存、事件与流、cuda库几个部分，全文有几万字（长文预警！！）。

1. cuda模型简介

编程模型是为了方便程序员使用更高级的语言进行开发，学习cuda编程首先要对cuda的2套模型有基本的了解，分别是cuda内存模型和线程模型。

内存模型：用于在异构环境下对主机与设备内存进行管理，对设备内存按访问等级和性能高低抽象出了全局内存、共享内存、寄存器内存等概念； 提供了分配和管理内存的关键字和API。

线程模型：设计了线程的层次结构，便于线程资源的申请与管理；通过一系列的内置参数和API提供了线程资源申请、线程内存访问以及线程执行控制的方法。

1.1 内存模型

cuda编程模型将异构环境分为主机端(host)和设备端(device)，典型的cuda程序（单个cuda流）执行分为5个步骤：

1. 分配设备端内存;
2. 拷贝数据到设备端内存;
3. 调用kernel在设备上进行计算;
4. 将数据从设备端内存拷回到主机端CPU内存;
5. 释放设备端内存。

内存模型中设备全局内存的管理和端到端数据传输可以通过在host端调用内置的cuda API来实现，cuda内置的内存管理API与标准的C函数类似，常用cuda内存申请和管理的API如下：

c函数	cuda函数	作用
malloc	cudaMalloc	分配内存
memcpy	cudaMemcpy	数据传输
memset	cudaMemset	初始化内存
free	cudaFree	释放内存

在host端调用这些API可以对device端的全局内存进行管理，这些API配合使用的方法如下：

// 定义数据大小和占用字节数
int nElem = 1<<20;
size_t nBytes = nElem*sizeof(float)

// float *h_A = new float[nElem];
float *h_A = (float *)malloc(nBytes); //分配 host内存
initData(h_A, nElem); //host内存数据初始化，initData为自定义的初始化函数

float *d_A
cudaMalloc((float **)&d_A, nBytes); //分配 device端GPU全局内存

cudaMemcpy(d_A, h_A, nBytes, cudaMemcpyHostToDevice); //传输数据到device端

// delete[] h_A;
free(h_A); // 释放host内存
cudaFree(d_A); // 释放device内存

这段代码中有两个值得注意的细节：

1. cudaMalloc函数第1个参数类型为二重指针，说明cudaMalloc函数要在指针d_A中存入一个地址，这个地址就是我们在device端所分配到的全局内存的首地址。
2. cudaMemcpy函数第4个参数cudaMemcpyHostToDevice，它是枚举类型cudaMemcpyKind的一个枚举值，代表内存拷贝数据传输的方向为host到device，此外通过cudaMemcpyDeviceToHost、cudaMemcpyDeviceToDevice可以分别实现device到host和device到device的传输。

上面介绍的都是一些经典款的内存管理API，高版本的cuda中还有其他更高级内存管理方法和新API，我们将会在后续cuda内存章节中进行详细整理。

1.2 线程模型

线程的层次结构

在cuda编程模型中线程具有层次结构，由内到外依次为：线程(thread)、块(block)、网格(grid)；定义的kernel是在grid级别进行启动和执行的。线程、块、网格这些逻辑上的模型与GPU现实物理层有一定的对应关系:

1. 模型中的最小执行单元thread对应物理层的Core，也叫SP（streaming processor）；
2. 物理层面没有像模型层面grid、block这样二级的阵列结构，只有一级SM（Stream Multi-processor）阵列，每个SM上有32的倍数个Core，单个block只能分布在单个SM上，cuda一个block最多能申请线程个数有限制，如1024。在为kernel申请计算资源时，应考虑所使用型号GPU中SM个数以及每个SM的线程数量，一个基本原则是让一个SM恰好能被多个block的线程占满。

线程识别

给每个thread编上号可以便捷的进行定位和识别，cuda线程模型中具体的实现方式是通过threadIdx，blockIdx等内置变量计算出一个全局的线程ID(gid)。

采用内置的dim3类型定义grid和block结构时，第1个维度具有最高优先级。假如我们定义的是一个2维的结构，维度1对应行，维度2对应列，第1个维度具有高优先级也就是通俗说的行优先，这与我们平时二维数组列优先的习惯相悖，为什么会这样设计呢？这是由于dim3类型构造器支持参数缺省造成的：dim3最多可以定义3个维度参数，靠前的参数优先赋值，缺省的参数默认值为1。因此在计算线程全局ID时，越靠前的参数也具有更高的优先级。

无论是1维、2维还是3维的阵列结构，计算gid的要旨如下：

1. 计算线程在对应block中的ID(tid);
2. 计算对应block中线程0在网格中的ID(gid_0)，计算公式为gid_0 = bid*b_size, 其中bid为block在网格中的全局ID，b_size为网格大小（单个网格线程总数）;
3. tid + gid_0 = gid;

具体的，对于1维结构：

tid = threadIdx.x, bid = blockIdx.x, b_size=blockDim.x

gid_0 = bid * b_size = blockIdx.x * blockDim.x

gid = tid+gid_0 = threadIdx.x + blockIdx.x * blockDim.x

2维结构:

tid = threadIdx.x + threadIdx.y * blockDim.x;

bid = blockIdx.x + blockIdx.y * gridDim.x;

b_size = blockDim.x * blockDim.y;

gid = tid + gid_0 = threadIdx.x + threadIdx.y * blockDim.x + (blockIdx.x + blockIdx.y * gridDim.x) * blockDim.x * blockDim.y;

对于2维网格和块，也可以在2个维度上分别计算gid，计算时将每个维度看作1维的：

gidx = threadIdx.x + blockIdx.x * blockDim.x;

gidy = threadIdx.y + blockIdx.y * blockDim.y;

通过gidx和gidy这一对ID同样也可以定位thread。

在给thread编上号实现了定位后，下一步就是分配各个thread去读取设备内存中的数据进行计算。每个thread一次可以读取设备内存中单个地址的数据。thread的gid与设备内存中的数据地址偏移（索引idx）都具有连续性，不同的gid，idx映射关系对应着不同的内存加载、写入方式，具体的访问方式后续章节再介绍。下面我们先介绍一个在2维块、网络上对2维数组进行求和的例子，用于了解cuda线程ID的识别。

template <typename T>
__global__ void sumArrayOnGPU2D(T *A, T *B, T *C, unsigned int nx, unsigned int ny){
    unsigned int ix = threadIdx.x + blockIdx.x * blockDim.x;
    unsigned int iy = threadIdx.y + blockIdx.y * blockDim.y;
    unsigned int idx = iy*nx + ix;

    if(ix<nx && iy<ny){
        C[idx] = A[idx] + B[idx];
    }
}

上述sumArrayOnGPU2D核函数中先定义了ix, iy用于定位网格中的线程，然后使用idx定义线程在内存地址上的分配逻辑。如果我们定义的是2x2的grid与block，对应的线程编号gid和内存偏移idx分别为：

图1.1 二维网格ix, iy, gid, idx对应关系

由于内存地址不像thread id一样具有2级层次结构，在这种分配方式下，同一个block的线程不是分配在一块连续的内存上。block在行上(x)尺寸不一样，内存加载的性能也会不一样。

线程执行

cuda中线程是采用单指令多线程（SIMD）的方式组织工作，kernel执行模型提出了线程束(warp)的概念，一个线程束由32个连续的线程构成，是block中单条指令执行的基本单位。因此sumArrayOnGPU2D例子中block行的尺寸应尽可能的设为32的倍数。可以在kernel中通过条件分支指定一个线程束中哪些线程执行指令、哪些线程不执行指令，这种方式称之为线程束分化，避免非必要的线程束分化可以有效提高kernel性能。

在device端的kernel中调用__syncthreads方法可以实现块级别线程同步，在host端调用cudaDeviceSynchronize可以实现系统级别的线程同步，即对所有线程进行同步。跨越块的线程不允许进行块间同步。

并行规约——线程并行与同步

下面介绍一个求一维数组元素之和的规约例子，充分的理解线程执行过程。在CPU编程中采用分治的思想求一维数组的元素之和是一种时间复杂度较低(logN)的算法，实现方式如下：

int reduceSumOnCPU(int *data, int const size){
    if(size==1){
        return data[0];
    }
    int const stride = size/2;

    for(int i=0; i<stride; i++){
        data[i] += data[i+stride];
    }
    return reduceSumOnCPU(data, stride);
}

上述reduceSumOnCPU先将数组分成2个等长部分（忽略奇偶），然后2部分对应元素相加，重复这2步，直到所有元素都聚合到索引为0的位置，最后返回。同理，在kernel中我们可以用一个循环，每一步循环将数据分成两部分，然后用1/2、1/4、1/8...的thread参与运算来实现规约，在每次循环需要考虑以下3点：

1. 让哪部分线程参与运算？
2. 加载全局内存中的哪部分数据？
3. 参与运算的线程怎么保证同步性？

对于问题1，我们希望每次循环所参与运算的线程是连续的，如果参与运算的线程不连续，这样每32个线程所构成的线程束中会有一部分参与运算，一部分不参与运算，造成线程束分化。同样问题2中加载的内存当然也希望是连续的，这样内存操作的效率更高，问题3由于只对部分参与运算的线程进行同步，在前面提到的两种线程同步方式中只有块级别同步满足对部分线程同步的要求，因此在每次循环中我们使用__syncthreads内置方法对参与运算的线程进行同步，这样就导致kernel中并行规约计算是block级别的，并不是全量数据级别的，所以我们在device端对每个块各求一个和，再在host端再对每个块的和进行相加，得到最终结果。其中使用每个块并行规约的kernel如下：

__global__ void reduceSumOnGPU(int *g_idata, int *g_odata, unsigned int n){
    unsigned int tid = threadIdx.x;
    unsigned int idx = blockIdx.x*blockDim.x + threadIdx.x;
    int *idata = g_idata + blockDim.x*blockIdx.x;

    if(idx>=n){
        return;
    }

    for(int stride=blockDim.x/2; stride > 0; stride >>= 1){
        if(tid<stride){
            idata[tid] += idata[tid+stride];
        }
        __syncthreads();
    }
    if(tid==0){
        g_odata[blockIdx.x] = idata[0];
    }
}

在上述reduceSumOnGpu的kernel传入参数中，指针g_idata指向输入数据的全局内存首地址，g_odata指向用于存放输出数据的内存首地址，n为元素的个数。在kernel中定义了一个idata的指针，这个指针使用的是线程寄存器，通俗点讲在每个线程的寄存器中存放了一个地址idata，idata指向的是这个线程所在block所负责的数据块的首地址。然后在循环体中通过变量stride控制每次让block中1/2、1/4、1/8的连续线程参与运算，对idata处2个长度分别为1/2blcokSize, 1/4blockSize的连续内存进行操作，每次运算完成后使用__syncthreads()进行同步。循环终止后将结果写入g_odata所指向的内存中。

对于以上reduceSumOnGPU的kernel，我们还可以进行2点优化：

1. reduceSumOnGPU中每个block只处理内存中的单个数据块，在这样的方式下数据吞吐量不高，可以让单个block处理多块数据，提升单个block内存操作的效率；
2. 每个block中可以申请的thread上限是固定的1024，因此循环迭代的上限是已知的log1024 = 10，具体的迭代次数根据blockSize而定，展开循环可以减少至少每轮一次的边界判定和其他循环维护的指令，优化kernel执行性能。

进行以上2点优化后代码如下：

template <unsigned int iBlockSize>
__global__ void reduceSumUnroll(int *g_idata, int *g_odata, unsigned int n){
    unsigned int tid = threadIdx.x;
    unsigned int idx = blockDim.x*blockIdx.x*4 + threadIdx.x;

    int *idata = g_idata + blockIdx.x*blockDim.x*4;

    if(idx+3*blockDim.x<n){
        int i0 = g_idata[idx];
        int i1 = g_idata[idx+blockDim.x];
        int i2 = g_idata[idx+blockDim.x*2];
        int i3 = g_idata[idx+blockDim.x*3];
        g_idata[idx] = i0+i1+i2+i3;
    }

    __syncthreads();

    if(iBlockSize>=1024&&tid<512){
        idata[tid] += idata[tid+512];
    }
    __syncthreads();

    if(iBlockSize>=512&&tid<256){
        idata[tid] += idata[tid+256];
    }
    __syncthreads();

    if(iBlockSize>=256&&tid<128){
        idata[tid] += idata[tid+128];
    }
    __syncthreads();

    if(iBlockSize>=128&&tid<64){
        idata[tid] += idata[tid+64];
    }
    __syncthreads();

    if(tid<32){
        volatile int *vmem = idata;
        vmem[tid] += vmem[tid+32];
        vmem[tid] += vmem[tid+16];
        vmem[tid] += vmem[tid+8];
        vmem[tid] += vmem[tid+4];
        vmem[tid] += vmem[tid+2];
        vmem[tid] += vmem[tid+1];
    }

    if(tid==0){
        g_odata[blockIdx.x] = idata[0];
    }
}

上述代码中模板参数iBlockSize用于指定块的大小，3~16行实现了每个block先将g_idata中4个连续数据块的数据聚合。定义的变量idx记录每个block负责的第一个数据块在g_idata中的起始索引，idata记录每个block负责数据块的首地址，idx与idata有异曲同工之妙，但分别用于不同的用途，假设blockSize = 4,对应原理如1.2所示：

图1.2 单个block处理多块内存的tid, block id, idx对应关系图

此时每个block负责数据的总长度为blockDim.x * 4，编号x的block负责起始地址：idata = g_idata + blockIdx.x * blockDim.x * 4，长度为blockDim.x * 4的数据块，block中编号x的thread需要负责聚合以下4个位置上的数据:

idata+threadIdx.x = g_idata[idx]

idata+threadIdx.x+blockDim.x = g_idata[idx+blockDim.x]

idata+threadIdx.x+2 * blockDim.x = g_idata[idx+2 * blockDim.x]

idata+threadIdx.x+3 * blockDim.x = g_idata[idx+3 * blockDim.x]

经过16行的__syncthreads()之后，长度为4 * blockSize的数据被聚合到了长度为blockSize的数据块中，聚合后数据块的首地址为idata，离得最近的两个数据块的idata相隔4 * blockSize，然后18~36行4个并行分支根据blockSize的不同对聚合后的数据块进行规约，将每个数据块的数据一步步规约到长度为64的更小数据块中。由于每个线程束的大小为32，这是指令执行的最小单位，操作单个线程束内的线程自动进行隐式同步，无需再调用__syncthreads()。在38~46行中接着通过6次规约将数据聚合到idata[0]上，值得注意的是在39行使用volatile关键字定义了一个指针，volatile关键字的作用是告诉编译器：通过vmem指针偏移进行读取或写入的数据是会经常变动的，不要在编译阶段进行优化。举个例子：

vmem[tid] += vmem[tid+32]
vmem[tid] += vmem[tid+16];

这两行代码如果不加volatile关键字，编译器可能会进行这样的优化：先从vmem+tid+32和vmem+tid+16处进行取值，然后将两个值相加之后再写入vmem+tid处，这样可以减少一次内存写入，显然在thread并行执行时会造成计算结果错误。在此处我们更希望每一次内存的读写都不要省略。

2. cuda内存

cuda内存同样也具有层次结构，根据访问层级、物理上存储位置、使用功能可以分为：全局内存、共享内存、寄存器内存、本地内存、常量内存和纹理内存。

全局内存：位于DRAM存储器上，最常用的内存，GPU显存大部分都会用作全局内存，可以被所有的线程访问，延迟高、带宽低、容量大。

共享内存：位于SRAM存储器上，每个block有自己的共享内存，能被block自己的线程访问，在物理上与L1缓存属于同一个级别，延迟低，带宽高，容量小，不同的设备48~100多KB不等，使用__shared__关键字进行定义。

寄存器内存：每个thread独享，运行速度最快。在Kernel中定义的没有关键字修辞的变量所使用的内存都属于寄存器内存。寄存器属于稀缺资源，不同GPU框架下的GPU个数从几十个到几百个不等。

本地内存：线程中定义的变量个数过多，寄存器内存不够用时会溢出到本地内存中进行定义，计算能力2.0以下的GPU架构中，本地内存在物理上与全局内存享用同一块存储区域，位于DRAM上，属于高延迟，低带宽的内存空间，在计算能力2.0以上GPU中，本地内存位于SM的一级缓存与设备的二级缓存中，位于SRAM上。

常量内存：常量内存在物理上与全局内存一样位于DRAM存储器中，区别在于它有一个专用的片上缓存，也就是说对常量内存在进行第1次读取时性能与全局内存一样差，但是在第2次及之后的读取由于有专属片上缓存，性能很快，所以常量内存适合存放需要进行反复读取的数据。常量内存以全局变量的方式进行声明，使用__constant__关键字进行修饰。所有的kernel都可以访问，由于需要用到专属缓存，所以容量小，所有设备都只有64KB。

纹理内存：是一种只读访问的全局内存，与常量内存一样存在片上缓存，在一些特定的2维数据场景下具有性能优势。

以上6种内存中全局内存、共享内存、寄存器内存和常量内存使用频率最高，其中全局内存、共享内存以及常量内存的使用存在一定的技巧，下面将详细介绍这3种内存的使用。

2.1 全局内存

全局内存的管理

在1.1节中我们已经介绍了部分内存管理和数据传输的内容，在此，我们继续补充几种内存管理和数据传输的手段。

固定主机内存： 采用标准的C函数如malloc分配的host内存是可分页的，没法在deviced端进行访问，使用cuda内置的cudaMallocHost、cudaHostAlloc函数可以分配固定的（锁页内存）host内存，在device进行访问。用cudaFreeHost进行释放，值得注意的是分配过多固定的host内存会影响host的性能。cudaHostAlloc多一个枚举类型的参数flags，用来设置内存分配的类型：

cudaHostAllocDefalt -- 等同于cudaMallocHost，默认值；

cudaHostAllocPortable -- 分配的内存对所有的GPU都可以使用；

cudaHostAllocWriteCombined -- 分配不使用L1和L2 cache的内存；主机端读会很慢，适合主机端写入，设备端读取的场景；

cudaHostAllocMapped -- 分配可以实现零拷贝功能的内存，主机端和设备端各维护一个地址，通过地址直接访问该块内存，无需传输。

零拷贝内存：无需显式拷贝就可以在device端使用，当device内存不足时利用主机内存，避免了cudaMemcpy这种显式的数据传输，PCIe传输率得到了提升，而且性能不受影响，但是数据传输（隐式拷贝）还是存在。零拷贝内存使用cudaHostAlloc指定flags为cudaHostAllocMapped进行申请，然后通过调用cudaHostGetDevicePointer将host指针与device指针关联，方法如下：

cudaHostAlloc(&h_A,nBytes, cudaHostAllocMapped);
initialData(h_A,nElem);
cudaHostGetDevicePointer(&d_A,h_A,0);
kernelFunc<<<grid, block>>>(&d_A, ...);

由于零拷贝内存每一次d_A到h_A的映射都需要PCIe进行数据传输，频繁读写会导致性能降低。

统一内存

为了简化全局内存的管理，后续版本GPU或CUDA可以对主机和设备内存进行统一托管。实现统一内存管理依赖两个特性：统一虚拟寻址和统一内存寻址。其中统一虚拟寻址支持计算能力2.0以上的GPU，统一内存寻址为cuda6.0以上版本新加入的特性。

统一虚拟寻址(UVA)：计算能力2.0及以上GPU支持，主机内存与设备内存可以共用一个虚拟的地址空间，可以隐去设备指针d_A的声明和定义，也省去了cudaHostGetDevicePointer函数将主机和设备的内存进行管理，直接可以将cudaHostAlloc获取的h_A传给核函数：

cudaHostAlloc(&h_A,nBytes, cudaHostAllocMapped);
initialData(h_A,nElem);
kernelFunc<<<grid, block>>>(&h_A, ...);

统一内存寻址：cuda6.0中加入的新特性，同样可以实现一次申请，随处调用。相比于统一虚拟寻址：统一内存寻址可以自动的对当个页面的数据进行移动，而采用统一虚拟寻址的方式管理的内存（零拷贝内存）不会自动的对内存进行移动，隐式依赖PCIe传输，性能会受影响。

托管内存采用统一内存寻址的方式进行管理，通过cudaMallocManaged可以分配托管内存，相比于采用统一虚拟寻址管理的零拷贝内存，可以自动的进行数据传输和重复指针消除。

访问模式

全局内存的访问要通过缓存，其中L1缓存不一定通过，但L2缓存一定会通过。可以通过编译指令：-Xptxas -dlcm=ca启用或 -Xptxas -dlcm=cg 禁用一级缓存。

一行一级缓存是128字节，线程对全局内存的访问以线程束为单位，如果每个线程请求4个字节的数据，如果启用了一级缓存，那么线程束中32个线程的一次访问刚好与一行一级缓存的大小能对上。如果禁用一级缓存，仅使用一级缓存，此时一行缓存的大小为32字节。

如果内存事务访问（128一行或32一行）的首地址是128或32的偶数倍时（为什么是偶数倍，cuda基本数据类型至少需要双字节的存储空间？？），可以实现对齐访问。如果一个线程束的32个线程访问一段连续的内存，可以实现合并访问。对齐与合并访问可以提升单个内存操作事务的效率，一次性尽可能多的加载或写入数据。

全局内存实现矩阵转置

与合并访问对应的是交叉访问，交叉访问是一种性能极差的访问方式。下面通过一个矩阵转置的例子，理解合并访问与交叉访问：

1. 行主序读取，列主序写入

__global__ void transposeNaiveReadRow(float *in, float *out, const int nx, const int ny){
    unsigned int ix = threadIdx.x + blockDim.x*blockIdx.x;
    unsigned int iy = threadIdx.y + blockDim.y*blockIdx.y;

    if(ix<nx&&iy<ny){
        out[ny*ix + iy] = in[nx*iy + ix];
    }
}

2. 列主序读取，行主序写入

__global__ void transposeNaiveReadCol(float *in, float *out, const int nx, const int ny){
    unsigned int ix = threadIdx.x + blockDim.x*blockIdx.x;
    unsigned int iy = threadIdx.y + blockDim.y*blockIdx.y;

    if(ix<nx&&iy<ny){
        out[nx*iy + ix] = in[ny*ix + iy];
    }
}

行主序读取、列主序写入在读取时是合并访问模式，列主序读取、行主序写入在写入时是合并访问模式。通常由于缓存的存在，在直接对全局内存使用矩阵转置的例子中列读取、行写入相比于行读取、列写入性能要更好。

2.2 共享内存

共享内存属于可编程的一级缓存，前面已经提到共享内存的部分特性：块内共享、L1缓存级别的高性能。在cudaz中共享内存有如下几个用途：

1. 作为可编程的全局内存的缓存，优化全局内存的访问模式；
2. 在块内线程间共享，作为块内线程间数据交互的通道。

共享内存的分配

共享内存在kernel中申请，通过__shared__关键字来分配，可以通过指定size静态分配1维、2维（最高支持3维）的共享内存：

__shared__ float arr[size_y][size_x];

也可以不指定size，使用extern关键字延后到运行阶段再动态分配内存：

extern __shared__ float arr[];

...
kernelFunc<<<grid, block, nByte>>>();

动态分配的共享内存在声明时只能声明成1维的。

共享内存的访问

共享内存共有32个存储空间，每个存储空间称为一个bank，32个bank恰好能被一次内存事务线程束中的32个线程同时访问。如果在1个内存事务中，1个线程束的多个线程访问了同1个bank会需要排队，发生bank冲突导致1个事务被划分为多个，造成访问效率的降低。为了避免bank冲突，cuda内部有广播机制，即当发生bank冲突时只让冲突的线程中的1个进行访问，然后广播到bank下其他的冲突线程中，这种方式虽然可以消除bank冲突，但多个线程中只有1个进行了数据的访问，带宽利用率还是较差。

下面共享内存列主序写入、行主序读取的例子在将数据写入共享内存时存在存储体冲突，在这个例子中BDIMX=32, BDIMY=32, 申请的二维共享内存为方形。：

__global__ void setColReadRow(int *out){
    __shared__ int tile[BDIMY][BDIMX];
    unsigned int idx = threadIdx.y * blockDim.x + threadIdx.x;

    tile[threadIdx.x][threadIdx.y] = idx;
    __syncthreads();

    out[idx] = tile[threadIdx.y][threadIdx.x];
}

setColReadRow第5行将idx写入共享内存tile时采用了列主序访问，存在bank冲突，加入块的行数BDIMX=32，块内thread id 0, 1, 2, 3对应的threadIdx.x, threadIdx.y分别为：

(0, 0), (1, 0), (2, 0), (3, 0)

访问共享内存为：

tile[0][0], tile[1][0], tile[2][0], tile[3][0]

由于BDIMX=32，对应全局偏移idx为：

0, 32, 64, 96

如下图2.1所示，thread 0, 1, 2, 3显然在同一个warp内，而且都访问了 bank 0，存在bank冲突。

图2.1 setColReadRow中BDIMX=32列主序访问tid, bank id, idx对应关系图

setColReadRow第5行读取tile值然后写入全局内存时采用行主序访问，不存在bank冲突，thread id 0, 1, 2, 3访问内存为：

tile[0][0], tile[0][1], tile[0][2], tile[0][3]

对应全局偏移idx为：

0, 1, 2, 3

如下图2.2所示，不存在bank冲突。

图2.2 setColReadRow中BDIMX=32行主序访问tid, bank id, idx对应关系图

假如改变方形数组的size，定义BDIMX = 16或BDIMX = 64会发生什么？

BDIMX = 16: thread 0, 1, 2, 3访问的idx分别为 0, 16, 32, 48, 其中thread 0和2落在一个bank 0内，thread 1和3落在了bank 16内，同样也存在bank冲突。

BDIMX = 64: thread 0, 1, 2, 3访问idx分别为0, 64, 128, 192，同样也存在bank冲突。

一种主动解决bank冲突的方式是对共享内存进行填充：申请共享内存时特意留出一些行，这样共享内存的size比block的大，在访问共享内存时线程id和实际访问的地址索引错位，从而避免多个thread落在同一个bank中，具体实现方式如下：

__global__ void setColReadRowPad(int *out){
    __shared__ int tile[BDIMY][BDIMX+IPAD];
    unsigned int idx = threadIdx.y*blockDim.x + threadIdx.x;
    
    tile[threadIdx.x][threadIdx.y] = idx;
    __syncthreads();
    out[idx] = tile[threadIdx.y][threadIdx.x];
}

setColReadRowPad中申请共享内存tile时列数增加了IPAD填充，定义BDIMX = 32, IPAD = 1，此时代码第5行列主序访问共享内存时thread 0, 1, 2, 3访问idx分别为

0, 33, 66, 99

如图2.3所示，idx 0, 33, 66, 99分别落在不同bank内，不存在bank冲突。

第6行进行行主序访问时，thread 0, 1, 2, 3对应idx仍为 0, 1, 2, 3。但线程 31, 32, 33对应的threadIdx.x, threadIdx.y为： (31, 0), (0, 1), (1, 1)

对应idx为31, 33和34，与没有填充存在差异。

上面介绍的setColReadRow的例子是针对方形2维共享内存，即BDIMX=BDIMY的情形，实际情况中使用的共享内存往往是矩形(BDIMX!=BDIMY)，或是动态声明的1维共享内存。针对矩形共享内存，setColReadRowPad中行主序tile[threadIdx.y][threadIdx.x]的方式可以正常访问，但是列主序tile[threadIdx.x][threadIdx.y]的方式就会出现下标越界，针对列主序，需要手动的计算行和列的索引：

unsigned int irow = idx/blockDim.y;
unsigned int icol = idx%blockDim.y;

此时矩形共享内存带填充的列主序写、行主序读recSetColReadRowPad变为：

__global__ void recSetColReadRowPad(int *out){
    __shared__ int tile[BDIMY][BDIMX+IPAD];
    unsigned int idx = threadIdx.y*blockDim.x + threadIdx.x;
    
    unsigned int irow = idx/blockDim.y;
    unsigned int icol = idx%blockDim.y;
    
    // transpose irow, icol
    tile[icol][irow] = idx;
    __syncthreads();
    
    out[idx] = tile[threadIdx.y][threadIdx.x];
}

如果使用的带填充的动态声明的一维共享内存，需要手动的计算带动态填充的1维索引： 列主序写的1维索引：

unsigned int col_idx = icol*(blockDim.x+IPAD) + irow;

行主序读的1维索引：

unsigned int row_idx = threadIdx.x + threadIdx.y*(blockDim.x+IPAD);

列主序写、行主序读recSetColReadRowDynPad为：

__global__ void recSetColReadRowDynPad(int *out){
    extern __shared__ int tile[];
    unsigned int idx = threadIdx.y*blockDim.x + threadIdx.x;
    
    unsigned int irow = idx/blockDim.y;
    unsigned int icol = idx%blockDim.y;
    
    unsigned int col_idx = irow + icol*(blockDim.x+IPAD);
    unsigned int row_idx = threadIdx.x + threadIdx.y*(blockDim.x+IPAD);
    
    tile[col_idx] = idx;
    __syncthreads();
    
    out[idx] = tile[row_idx];
}

调用时，延后传入的共享内存size = BDIMY * (BDIMX+IPAD) * sizeof(int)。

共享内存的宽度与容量

共享内存bank宽度默认为4字节，计算能力3.x及以上的GPU支持8字节宽度的bank，允许两个线程同时访问同一个bank的不同子字。因此，在相同的访问模式下，8字节64位宽度的存储器会产生更少的bank冲突。

使用cudaDeviceSetSharedMemConfig(cudaSharedMemConfig config)可以设置存储器宽度，cudaSharedMemConfig枚举值的选项有：cudaSharedMemBankSizeDefault, cudaSharedMemBankSizeFourByte, cudaSharedMemBankSizeEightByte；更改bank size不会影响共享内存的使用量，也不影响核函数的占用率。

共享内存与一级缓存共用硬件资源，区别在于共享内存以32个bank进行一次访问，而L1是以缓存的形式进行访问，CUDA为配置共享内存和L1缓存的大小提供了两种方法：

1. device级别配置：使用cudaDeviceSetCacheConfig(cudaFuncCache cacheConfig), 参数支持cudaFuncCachePreferNone-默认，cudaFuncCachePreferShared-划分更多的共享内存，cudaFuncCachePreferL1-划分更多的L1缓存，cudaFuncCachePreferEqual-等量划分共享内存和L1缓存。
2. kernel级别配置：使用cudaFuncSetCacheConfig(const void *func, enum cudaFuncCache cacheConfig)，可以对单一的kernel进行共享内存和L1缓存的配置。

使用共享的并行规约

前面章节中介绍了完全使用全局内存进行并行规约的例子，如果在并行规约中使用到共享内存，可以明显的减少对全局内存的访问，使性能得到进一步提升。使用共享内存、展开的并行规约kernel如下：

__global__ void reduceSmem(int *g_idata, int *g_odata, unsigned int n){
    extern __shared__ int smem[];

    unsigned int tid = threadIdx.x;
    unsigned int idx = threadIdx.x + blockIdx.x*blockDim.x*4;

    int tempSum = 0;
    if(idx+3*blockDim.x<=n){
        int a1 = g_idata[idx];
        int a2 = g_idata[idx+blockDim.x];
        int a3 = g_idata[idx+2*blockDim.x];
        int a4 = g_idata[idx+3*blockDim.x];
        tempSum = a1+a2+a3+a4;
    }

    smem[tid] = tempSum;
    __syncthreads();

    if(blockDim.x>=1024&&tid<512){
        smem[tid] += smem[tid+512];
    }
    __syncthreads();

    if(blockDim.x>=512&&tid<256){
        smem[tid] += smem[tid+256];
    }
    __syncthreads();

    if(blockDim.x>=256&&tid<128){
        smem[tid] += smem[tid+128];
    }
    __syncthreads();

    if(blockDim.x>=128&&tid<64){
        smem[tid] += smem[tid+64];
    }
    __syncthreads();

    if(tid<32){
        volatile int *vsmem = smem;
        vsmem[tid] += vsmem[tid+32];
        vsmem[tid] += vsmem[tid+16];
        vsmem[tid] += vsmem[tid+8];
        vsmem[tid] += vsmem[tid+4];
        vsmem[tid] += vsmem[tid+2];
        vsmem[tid] += vsmem[tid+1];
    }

    if(tid==0){
        g_odata[blockIdx.x] = smem[0];
    }
}

在reduceSmem中使用了动态声明了共享内存smem，相比于前面章节的核函数reduceSumUnroll，使用smem代替了指向全局内存的指针idata，减少了对全局内存的访问。19~47行直接在smem中执行块内原地规约，速度显然更快。

使用共享内存的矩阵转置

前面章节中介绍的矩阵转置核函数transposeNaiveReadRow, transposeNaiveReadCol在写或读的时候存在对全局内存的交叉访问，使用2维的共享内存可以有效的避免这种交叉访问。前面介绍的recSetColReadRow蕴含共享内存的转置的思想，但这种转置是在单个块内进行的，现实中的需要转置的矩阵往往非常大，在单个块内的共享内存中肯定放不下，需要借用多个块构成的grid实现，逻辑方式会稍微复杂一点：

__global__ void transposeSmem(int *in, int *out, int nx, int ny){
    __shared__ int tile[BDIMY][BDIMX];

    unsigned int ix, iy, ti, to;
    ix = blockIdx.x*blockDim.x + threadIdx.x;
    iy = blockIdx.y*blockDim.y + threadIdx.y;
    ti = iy*nx + ix;

    unsigned int bidx, brow, bcol, zix, ziy;
    bidx = threadIdx.y * blockDim.x + threadIdx.x;
    brow = bidx / blockDim.y;
    bcol = bidx % blockDim.y;
    zix = blockDim.y*blockIdx.y + bcol;
    ziy = blockDim.x*blockIdx.x + brow;

    to = ziy*ny + zix;
    if(ix < nx && iy < ny){
        tile[threadIdx.y][threadIdx.x] = in[ti];
        __syncthreads();
        out[to] = tile[bcol][brow];
    }
}

transposeSmem核函数前面定义了一大串的索引，我们不妨从后往前看，在17行到20行有2个过程，1. 将全局内存in的数据读入到共享内存tile，2.同步之后读取tile的数据到全局内存out。

第1步是按照常规的索引进行的，全局内存指针in的索引ti是行优先，共享内存索引tile[threadIdx.y][threadIdx.x]也是常规的行优先，假如grid是2x4的，block是4x8的，读取完成后ti, ix, iy对应关系如图2.3 所示。

图2.3 2x4 grid, 4x8 block常规索引对应图

写入tile阶段全局索引ti=106处数据负责的线程为：

blockIdx.x, blockIdx.y: (0,1) block中的

threadIdx.x, threadIdx.y: (2,5) 线程

这个线程在读取阶段计算出的索引brow,bcol为(2,6)对应的ti为114：

// block 4x8, blockDim.x = 4, blockDim = 8
// 二维索引(2,5)在块内对应全局索引bidx = 5*4+2 = 22;
bidx = threadIdx.y * blockDim.x + threadIdx.x; // 22


// 写入时读取的位置
brow = bidx / blockDim.y; // 2
bcol = bidx % blockDim.y; // 6

块内的线程在读和写阶段分别负责ti如图2.4所示：

图2.4 同一个线程写和读对应的ti

此时新的ix, iy分别为:

zix = blockDim.y*blockIdx.y + bcol; // 14
ziy = blockDim.x*blockIdx.x + brow; // 2

to = 2 x 32 + 14 = 78

此线程写入out时应全局索引to=78，即out中偏移量78数据的写入与in中偏移量为106数据的读取使用的是同一个线程，各索引对应关系如下图2.5所示

图2.5 2x4 grid, 4x8 block转置过程中索引关系对应图

由图2.4可知，ti=106的数据在第2个阶段将由以下线程负责

blockIdx.x, blockIdx.y: (0,1) block中的

threadIdx.x, threadIdx.y: (1,5) 线程

计算bidx = 5x4+1 = 21, brow = 2, bcol = 5, tile[5][2] = 106;

zix = 8x1+5 = 13, ziy = 4x0 + 2 = 2, to = 2 x 32 + 14 = 77

transposeSmem没有使用填充，在第20行会存在bank冲突，带有填充的共享内存转置如下：

__global__ void transposeSmemPad(int *in, int *out, int nx, int ny){
    __shared__ float tile[BDIMY][BDIMX + IPAD];

    unsigned int ix, iy, ti, to;
    ix = blockIdx.x*blockDim.x + threadIdx.x;
    iy = blockIdx.y*blockDim.y + threadIdx.y;
    ti = iy*nx + ix;

    unsigned int bidx, brow, bcol, zix, ziy;
    bidx = threadIdx.y * blockDim.x + threadIdx.x;
    brow = bidx / blockDim.y;
    bcol = bidx % blockDim.y;

    zix = blockDim.y*blockIdx.y + bcol;
    ziy = blockDim.x*blockIdx.x + brow;
    to = ziy*ny + zix;

    if(ix < nx && iy < ny){
        tile[threadIdx.y][threadIdx.x] = in[ti];
        __syncthreads();
        out[to] = tile[bcol][brow];
    }
}

transposeSmemPad相比transposeSmem，只在共享内存定义时加了IPAD, 因为IPAD加在BDIMX上，BDIMY的size不变，不影响brow, bcol, zix, ziy的计算规则。

使用动态共享内存矩阵转置的实现如下：

__global__ void transposeSmemPadDyn(float *out, float *in, int nx, int ny){
    extern __shared__ float tile[];

    unsigned int ix, iy, ti, to, row_idx;
    ix = blockIdx.x*blockDim.x + threadIdx.x;
    iy = blockIdx.y*blockDim.y + threadIdx.y;
    ti = iy*nx + ix;
    row_idx = threadIdx.y * (blockDim.x + IPAD) + threadIdx.x;

    unsigned int bidx, brow, bcol, zix, ziy, col_idx;
    bidx = threadIdx.y * blockDim.x + threadIdx.x;

    brow = bidx / blockDim.y;
    bcol = bidx % blockDim.y;
    col_idx = bcol*(blockDim.x+IPAD) + brow;

    zix = blockDim.y*blockIdx.y + bcol;
    ziy = blockDim.x*blockIdx.x + brow;
    to = ziy*ny + zix;

    if (ix < nx && iy < ny){
        tile[row_idx] = in[ti];
        __syncthreads();

        out[to] = tile[col_idx];
    }
}

区别在于需要将2维的tile索引转为1维：

row_idx = threadIdx.y * (blockDim.x + IPAD) + threadIdx.x;

col_idx = bcol*(blockDim.x+IPAD) + brow;

此外，还可以将block的size减倍，通过展开在全局内存操作时一次性加载或写入更多的数据，提升带宽的利用率，具体的实现方式不在此详述了。

2.3 常量内存

前面已经介绍了常量内存和全局内存一样存储在设备的DRAM中，区别在于常量内存有一个片上的缓存。相较于其他内存，常量内存有一个最优的访问模式：一个线程束中的全部线程访问同一个位置。

常量内存变量通过__constant__关键字定义，与c/c++中全局变量一样，在全局定义，对所有源文件可见，声明周期涵盖整个进程，对device只读，对host可读可写，在host上通过cudaMemcpyToSymbol函数进行初始化：

#define RADIUS 4

__constant__ float coef[RADIUS+1];
    
void setup_coef_constant(){
    const float h_coef[] = {0.01f, 0.05f, 0.01f, 0.04f, -0.02f};
    cudaMemcpyToSymbol(coef, h_coef, (RADIUS+1)*sizeof(float));
}

常量内存适合用来存放一些需要反复读取的常量系数，下面通过一个1维的9点模板的例子来介绍常量内存的使用。模板计算是求解很多偏微分方程的基础，这在个例子中我们要实现一个如图2.6所示，

图2.6 9点模板计算示意图

其中out的每一个元素为in数组对应位置前后各4个元素的值计算得到，定义计算公式如下：

void cpu_stencil_1d(float *in, float *out, int n){
    for(int i = 0; i<n; i++){
        float temp = a1 * (in[i + 1] - in[i - 1])
                    + a2 * (in[i + 2] - in[i - 2])
                    + a3 * (in[i + 3] - in[i - 3])
                    + a4 * (in[i + 4] - in[i - 4]);
        out[i] = temp;
    }
}

这个函数在调用时h_in需要加上RADIUS偏置：

float *h_in = (float*)malloc(sizeof(float)*(n+2*RADIUS);
float *h_out = (float*)malloc(sizeof(float)*n);

cpu_stencil_1d(h_in+RADIUS, h_out, n);

相同9点模板计算逻辑在device端实现kernel如下：

__global__ void stencil_1d(float *in, float *out, int n){
    __shared__ float smem[BDIM + 2*RADIUS]; // reverse 2*radius

    int idx = threadIdx.x + blockIdx.x * blockDim.x;
    int sidx = threadIdx.x + RADIUS; // offset radius

    // read radius ~ radius+BDIM
    if(idx<n){
        smem[sidx] = in[idx];
    }

    // read 2 heads
    if(threadIdx.x<RADIUS){
        smem[sidx-RADIUS] = in[idx-RADIUS]; //0 ~ RADIUS
        smem[sidx+BDIM] = in[idx+BDIM]; //RADIUS+BDIM ~ RADIUS+BDIM+RADIUS
    }
    __syncthreads();

    float temp = 0.0f;
#pragma unroll
    for(int i=1; i<=RADIUS;i++){
        temp += coef[i]*(smem[sidx+i]-smem[sidx-i]);
    }

    out[idx] = temp;
}

在stencil_1d中定义了一个长度为BDIM + 2 * RADIUS的共享内存，8 ~ 16行进行了2次从全局内存in到共享内存smem数据的传输：

1. 8 ~ 10行为第1次传输，传输的是中间部分的数据，由块内所有线程参与，将in中长度为BDIM的数据写入smem中的RADIUS ~ RADIUS+BDIM处；
2. 第13 ~ 16行为第2次传输，传输了2段数据，只有thread id<RADIUS的线程参与，将in中长度为RADIUS的2段数据写入smem中对应的0 ~ RADIUS和 RADIUS+BDIM ~ RADIUS+BDIM+RADIUS处。

第13行if分支带来线程束分化，整个数据传输过程中，块内threadID 0~RADIUS-1的线程一共接受了3次指令，第1次为代码第9行参与传输长度为BDIM的中间段数据，第2次和第3次在代码第14、15行，分别独自进行两头长度为RADIUS数据的传输。3次指令是独立的，不需要进行同步；在传输完成后，进入计算之前进行1次线程同步（第17行）。

第19~25行为计算部分，20行定义的宏在预处理阶段提示编译器将循环展开，stencil_1d在调用的时候也需要对in加偏置：

stencil_1d<<<grid, block>>>(d_in + RADIUS, d_out, n);

stencil_1d 8~16行数据移动如图2.7所示

图2.7 stencil_1d数据移动简图，假定RADIUS=4, BDIM=8, n=16

2.4 寄存器内存

前面介绍了寄存器的简单用法，在kernel中不加任何修饰符定义的变量都使用的寄存器内存，寄存器为每个线程私有，生命周期与kernel相同。计算能力3.0以上的GPU支持同一个线程束的两个线程彼此去访问对方寄存器的值，这种相互访问通过线程束洗牌指令实现。

线程在block内通过threadIdx编号，在warp内通束内索引编号(laneID), 一维block下的warpID和laneID计算公式如下：

warpID = threadIdx.x/32

laneID = threadIdx.x%32

线程束洗牌指令

洗牌指令是一种内置的cuda函数，可以对束内寄存器进行操作，下面是一个常用的洗牌指令:

int __shfl_sync(unsigned int mask, int val, int srcLane, int width=warpSize)

__shfl_sync将线程束内ID为srcLane对应线程的val值从寄存器中取出来，然后赋值给束内其他线程，width参数可以控制在多大的范围内进行洗牌操作，默认值为warpSize，可以的取值还有16, 8, 4, 2。mask参数用来控制线程束分化，设置不同的mask值可实现线程束中不同的thread参与洗牌，在本文中我们不去深究mask的用法，各个例子都使用full mask，即让线程束中32个线程全都参与洗牌，这是为了和旧版本(cuda9.0以前的版本)中不带_sync后缀的洗牌指令兼容。如果width = 16, srcLane=2, 执行__shfl_sync后会将laneID=2线程的val值从寄存器中取出，然后赋给第0~15号线程，laneID=18线程的val值赋给第16~31号线程。使用__shfl指令可以实现线程束级的广播操作：

#define FULL_MASK 0xffffffff

__global__ void shfl_broadcast(int *d_in, int *d_out, int const srcLane){
    int value = d_in[threadIdx.x];
    value = __shfl_sync(FULL_MASK, value, srcLane);
    d_out[threadIdx.x] = value;
}

如果调用shfl_broadcast(d_in, d_out, 1)， shfl_broadcast在各线程中先定义value变量，前面介绍过这种在kernel中不加任何修饰符定义的变量是存放在各线程的寄存器中的。然后__shfl_sync洗牌指令在1号线程寄存器中将值取出，并重新赋给线程束内所有线程，此时所有线程的value变量值都与1号线程原先存放在寄存器中值相同，实现内线程束内的广播功能。

其他的线程束洗牌指令原理与__shfl_sync类型，只是功能上有些不同，常用的如下：

__shfl_up_sync(unsigned int mask, int val, unsigned int offset, int width = warpSize): 每个线程存放val变量的寄存器值复制到当前线程id+offset的线程中, 运算后id<delta的线程寄存器值不变；

__shfl_down_sync(unsigned int mask, int val, unsigned int offset, int width = warpSize): 每个线程存放val变量的寄存器值复制到当前线程id-offset的线程中, 运算后id>=width-offset的线程寄存器值不变；

__shfl_xor_sync(unsigned int mask, int val, unsigned int offset, int width = warpSize): 束内每个线程与id+offset线程交换存放val变量寄存器的值。

图2.8 洗牌指令执行示意图

以上每个洗牌指令都根据输入的val数据类型（float, double）不同进行了重载。

线程束洗牌指令实例

下面通过几个例子加深对洗牌指令的理解。

1. 环绕移动

下面test_shfl_wrap的例子将在线程束内实现以下功能

d_in：0, 1, 2, 3, 4 ...29, 30, 31

d_out：2, 3, 4, 5, 6 ...31, 0, 1

__global__ void test_shfl_wrap(int *d_in, int *d_out, int const offset){
    int val = d_in[threadIdx.x];
    val = __shfl_sync(FULL_MASK, val, threadIdx.x+offset);
    d_out[threadIdx.x] = val;
}

2. 交换线程中数组对应索引处的值

下面test_shfl_swap核函数在每个线程下声明了数组value, 然后在内联函数swap中交换了各线程数组value索引secondIdx对应值。定义blockSize = 4, SEGM = 4, 元素个数n = 16, 调用

block.x /= SEGM;

test_shfl_swap<<<grid,block>>>(d_in,d_out,1,0,3);

实现如下交换功能：

d_in: 0, 1, 2, 3, 4, 5, 6, 7, || 8, 9, 10, 11, 12, 13, 14, 15

d_out: 7, 1, 2, 3, 4, 5, 6, 0,|| 15, 9, 10, 11, 12, 13, 14, 8

#define SEGM 4

__inline__ __device__ void swap(int *value, int threadIdx, int offset, int firstIdx, int secondIdx){
    bool isEvenThread = ((threadIdx & 1) == 0);
    
    // 对偶数线程中、数组value、索引为firstIdx和secondIdx的值进行交换
    if(isEvenThread){
        int temp = value[firstIdx];
        value[firstIdx] = value[secondIdx];
        value[secondIdx] = temp;
    }
    
    // 用洗牌指令交换两个线程中、数组value对应索引为secondIdx的值
    value[secondIdx] = __shfl_xor_sync(FULL_MASK, value[secondIdx], offset);

    // 第二次交换偶数线程、数组value、firstIdx和secondIdx处的值
    if(isEvenThread){
        int temp = value[firstIdx];
        value[firstIdx] = value[secondIdx];
        value[secondIdx] = temp;
    }
}

__global__ void test_shfl_swap(int *d_in, int *d_out, const int offset, int firstIdx, int secondIdx){
    int idx = threadIdx.x * SEGM;
    int value[SEGM];
    for(int i=0; i<SEGM; i++){
        value[i] = d_in[idx+i];
    }
    swap(value, threadIdx.x, offset, firstIdx, secondIdx);
    for(int i=0; i<SEGM; i++){
        d_out[idx+i] = value[i];
    }
}

首先在每个线程中声明了一个长度为4的数组value，从d_in中加载16个数据到4个线程的value中。在swap中先判断threadIdx是否为偶数，如果是偶数，则交换数组value索引firstIdx、secondIdx的值，经过第6行到第10行后:

3, 1, 2, 0,|| 4, 5, 6, 7,|| 11, 9, 10, 8,|| 12, 13, 14, 15

第12行__shfl_xor_sync，交换了相邻两个线程中value数组索引为secondIdx值：

3, 1, 2, 7,|| 4, 5, 6, 0,|| 11, 9, 10, 15,|| 12, 13, 14, 8

最后再次对索引为偶数的线程，交换数组value索引firstIdx、secondIdx的值:

7, 1, 2, 3,|| 4, 5, 6, 0,|| 15, 9, 10, 11,|| 12, 13, 14, 8

3. 用洗牌指令实现并行规约

回顾前面使用共享内存进行规约，核函数reduceSmem运行分为2步：先将元素之和汇聚到32个点上，然后在线程束内对32个点进行规约最终汇聚到block内的1个点。第2步线程束内规约也可以使用洗牌指令__shfl_xor_sync实现。

int localSum = smem[tid];
if(tid<32){
    localSum += __shfl_xor_sync(FULL_MASK, localSum, 16);
    localSum += __shfl_xor_sync(FULL_MASK, localSum, 8);
    localSum += __shfl_xor_sync(FULL_MASK, localSum, 4);
    localSum += __shfl_xor_sync(FULL_MASK, localSum, 2);
    localSum += __shfl_xor_sync(FULL_MASK, localSum, 1);

}

if(tid==0){
    g_odata[blockIdx.x] = localSum;
}

此外，并行规约也可以使用2次束内规约进行实现：

1. 在warp级别进行第1次规约，每个warp汇聚一个值并将值存到共享内存中;
2. 将共享内存中的值读入到第一个warp中，在第一个warp中进行第2次规约，得到block内最终规约值。

这个方法有个限制条件就是每个block中warp的数量(warpNum)不大于warpsize 32，不然在第2次规约时共享内存的值在一个warp中放不下。由于单个block中线程数最多1024个，因此warpNum不会超过32。下面的核函数reduceShfl进行了实现：

#define WARP_NUM 16 // warpNum不能大于32

__inline__ __device__ int warpReduce(int mySum){
    mySum += __shfl_xor_sync(FULL_MASK,mySum, 16);
    mySum += __shfl_xor_sync(FULL_MASK,mySum, 8);
    mySum += __shfl_xor_sync(FULL_MASK,mySum, 4);
    mySum += __shfl_xor_sync(FULL_MASK,mySum, 2);
    mySum += __shfl_xor_sync(FULL_MASK,mySum, 1);
    return mySum;
}

__global__ void reduceShfl(int *g_idata, int *g_odata, unsigned int n){
    __shared__ int smem[WARP_NUM];

    unsigned int idx = threadIdx.x + blockIdx.x*blockDim.x;

    if(idx>=n){
        return;
    }

    int mySum = g_idata[idx];

    int laneIdx = threadIdx.x % warpSize;
    int warpIdx = threadIdx.x / warpSize;

    mySum = warpReduce(mySum);

    if(laneIdx==0){
        smem[warpIdx] = mySum;
    }

    __syncthreads();

    // WARP_NUM <= 32
    mySum = threadIdx.x < WARP_NUM ? smem[laneIdx] : 0;

    if(warpIdx==0){
        mySum = warpReduce(mySum);
    }

    if(threadIdx.x == 0){
        g_odata[blockIdx.x] = mySum;
    }
}

warpReduce使用洗牌指令实现了warp内规约。reduceShfl中，第21 ~ 30行进行了第1次规约，在块同步后进行了第2次规约。

3. cuda流与事件

cuda中kernel和API的执行是在网格或设备级别的，这些操作有些需要同步，有些可以异步进行，对这些cuda操作按同步情况和执行顺序进行调度管理可以有效的提高程序的执行效率。cuda中使用流和事件来对这些操作进行调度和监测。

3.1 流与事件的基本概念

cuda流介绍

cuda流是一系列需要前后相互依赖执行的cuda操作，属于一种更粗粒度的并发单位，流与流之间是可以完全异步的。cuda流分为两种类型：

空流: 默认的流，不需要显式的声明，前面接触到的在host端执行的cuda操作如数据传输、kernel执行都是在默认的空流下执行的；

非空流: 需要显式声明的流，将不需要相互依赖、可以异步执行的cuda操作放在不同的流中，不同流中的cuda操作可以在时间上重叠。

声明和创建流

非空流的声明、创建和销毁方式如下：

cudaStream_t stream; \\ 声明非空流
cudaStreamCreate(&stream); \\ 创建非空流
cudaStreamDestroy(stream); \\ 销毁(释放)非空流

执行cudaStreamDestory会立即返回，如果stream中还有未完成的操作，流中正在使用的资源并不会立即释放，只有当所有操作完成后才会释放。使用cudaStreamQuery(stream)可以检查流中操作是否已经完成，如果全部完成会返回cudaSuccess，如果未全部完成会返回cudaErrorNotReady。使用cudaStreamSyncronize(stream)可以进行流同步，cudaStreamSyncronize会阻塞主机，直到流中所有的操作完成。

非空流还可以指定优先级声明：

cudaStreamCreateWithPriority(cudaStream_t* pStream, unsigned int flags, int priority)

cudaStreamCreateWithPrioriy中第3个参数priority用来控制流的优先级，数值越小优先级越高。流的优先级只对kernel有效，高优先级流中的kernel可以在低优先级流中插队执行，但是主机与设备间的内存操作不受流的优先级影响。priority的取值有范围限制，使用cudaDeviceGetStreamPriorityRange(int *leastPriority, int *greatestPriority)可以获取priority的最小和最大取值，如果priority值超出了范围，则会被设置为最大或最小值。

第2个参数flags用来指定流的阻塞行为：指定显式声明的非空流会不会被默认的空流中的操作所阻塞。可选枚举值有

cudaStreamDefault： 默认值，显式定义的非空流会被空流所阻塞；

cudaStreamNonBlocking：显式定义的非空流不会被空流所阻塞，可以与空流隐式同步执行。

为流分配cuda操作

创建完流之后，kernel可以使用运算符<<<>>>中的第4个参数分配流，其他支持异步执行的cuda自带或第三方API中通过指定stream参数进行分配：

kernel_1<<<grid, block, 0, stream>>>(); // 分配kernel到stream中

cudaMallocAsync(&d_A, nBytes, stream); // 分配异步执行的内存操作到stream中

cuda事件介绍

cuda事件用来对cuda流中的各个执行点打上一个标记，用来监测流的执行，可以当作是加入到流中的一个cuda操作，但不会影响其他操作的执行。只有当流运行到事件标记节点后，事件才开始进行。cuda事件声明、创建和销毁方式如下

cudaEvent_t event; //声明事件
cudaEventCreate(&event) //创建事件
cudaEventDestory(evnet) //销毁事件

与流的销毁一样，只有当事件完成后对应的资源才会被释放。

使用cudaEventRecord(event, stream)向流中插入事件，stream可以缺省，默认为空流；

使用cudaEventSynchronize(event)可以阻塞主机等待流执行完事件；

使用cudaEventQuery(event)可以查询事件执行的情况，如果事件执行完成会返回cudaSuccess。

cuda事件一个常见的功能就是计算流中一段cuda操作的耗时：

cudaEvent_t start, stop;
CHECK(cudaEventCreate(&start));
CHECK(cudaEventCreate(&stop));

CHECK(cudaEventRecord(start, 0));
kernel_1<<<grid, block>>>();
CHECK(cudaEventRecord(stop,0));

CHECK(cudaEventSynchronize(stop));
float kernel_time;
CHECK(cudaEventElapsedTime(&kernel_time, start, stop));

在这段主机端代码中创建了2个事件start, end，并用cudaEventRecord将它们放置到默认的空流中kernel_1执行前后的2个节点上，然后同步到stop节点，使用cudaEventElapsedTime计算两个事件运行时间间隔的毫秒数，放到kernel_time中。

还可以通过设置阻塞行为来创建事件：

cudaError_t cudaEventCreateWithFlags(cudaEvent_t* event, unsigned int flags);

第2个参数flags控制事件的阻塞行为，枚举值可选有4种：

cudaEventDefault --默认，不阻塞
cudaEventBlockingSync --可以阻塞host
cudaEventDisableTiming --可以流间同步
cudaEventInterprocess --事件可以跨进程

3.2 使用流执行cuda操作

流中异步执行的cuda操作有2类：kernel函数和cuda API，其中cuda API以异步执行的内存操作为主。kernel和cuda API在执行时所处的硬件环境有很大的区别。因此流的异步执行可以分为：kernel的并发执行、kernel和cuda API的重叠执行。

kernel的并发执行

一般在GPU计算资源充足时，分配在不同流中的kernel执行是可以完全并发的。但是在某些老版的GPU架构中流的启动需要排队执行，不同流的kernel间看上去似乎存在先后依赖的关系，需要使用特定的方法消除这种虚假的依赖。

openMP

OpenMP是一种CPU并行编程模型，相关API位于<omp.h>头文件中，使用OpenMP可以在主机端使用多个主机线程启动kernel，让每个主机线程单独管理一个流。具体实现方法如下：

omp_set_num_threads(4);
#pragma omp parallel
{
    int i = omp_get_thread_num();
    kernel_1<<<grid,block, 0, streams[i]>>>();
    kernel_2<<<grid,block, 0, streams[i]>>>();
    kernel_3<<<grid,block, 0, streams[i]>>>();
    kernel_4<<<grid,block, 0, streams[i]>>>();
};

第一行使用omp_set_num_threads设置在主机端需要启动4个线程，编译指令#pragma omp parallel指定花括号间代码为并行代码块。在并行代码块中通过omp_get_thread_num()获取每个主机线程的id，然后用id将主机线程与流进行关联，每个流都将kernel_1~kernel_4执行一遍。

流间依赖

在比较复杂的cuda程序中，需要在一个流中阻塞cuda操作直到另一个流中的操作完成，这种流间的依赖关系需要借助cuda事件来实现，方式如下：

n_streams = 4;

cudaEvent_t *kernelEvent;
kernelEvent = (cudaEvent_t *) malloc(n_streams * sizeof(cudaEvent_t));

for (int i = 0; i < n_streams; i++){
    CHECK(cudaEventCreateWithFlags(&(kernelEvent[i]),cudaEventDisableTiming));
}

for (int i = 0; i < n_streams; i++){
    kernel_1<<<grid, block, 0, streams[i]>>>();
    kernel_2<<<grid, block, 0, streams[i]>>>();
    kernel_3<<<grid, block, 0, streams[i]>>>();
    kernel_4<<<grid, block, 0, streams[i]>>>();

    CHECK(cudaEventRecord(kernelEvent[i], streams[i]));
    CHECK(cudaStreamWaitEvent(streams[n_streams - 1], kernelEvent[i], 0));
}

代码3~8行给每个流创建了一个同步事件，然后执行4个kernel在每个流中同步执行，第16行将同步事件分配到流中，第17行使用cudaStreamWaitEvent控制最后一个流(streams[n_streams-1])需要等待其他流执行完后才能开始执行。

GPU的硬件资源是有上限的，当各kernel并行执行所需的资源总数超出GPU硬件资源上限时，流与流之间也会根据等待执行的kernel所需资源进行自动调配，出现不完全的并发。

kernel与cuda异步API重叠执行

与kernel重叠执行的异步API以数据传输的API为主，kernel与数据传输关系可以分为两种：

1. 传输的数据要在kernel上使用，那么数据传输需要与kernel位于同一个流中，且排在kernel前；
2. 传输的数据无需在kernel上使用，数据传输和kernel可以位于不同流中。

下面这段代码使用多个流对数组分段进行求和，位于不同流中kernel与数据传输是重叠执行的：

#define NSTREAM 4
int iElem = nElem / NSTREAM
int iBytes = iElem*sizeof(int)

for (int i = 0; i < NSTREAM; i++){
    int ioffset = i * iElem;
    CHECK(cudaMemcpyAsync(&d_A[ioffset], &h_A[ioffset], iBytes,
                          cudaMemcpyHostToDevice, stream[i]));
    CHECK(cudaMemcpyAsync(&d_B[ioffset], &h_B[ioffset], iBytes,
                          cudaMemcpyHostToDevice, stream[i]));
    sumArrays<<<grid, block, 0, stream[i]>>>(&d_A[ioffset], &d_B[ioffset],
            &d_C[ioffset], iElem);
    CHECK(cudaMemcpyAsync(&gpuRef[ioffset], &d_C[ioffset], iBytes,
                          cudaMemcpyDeviceToHost, stream[i]));
}

代码中将待求和的数组分成了4段，用4个流各负责一段数据的求和，每次循环中向单个流中分配：主机到设备数据传输--求和--设备到主机数据传输的操作，这种写法称之为深度优先分配。还有一种通过3个循环，每个循环分别给流分配：主机到设备数据传输--求和--设备到主机数据传输3个操作中的一个，这种写法称为广度优先分配：

for (int i = 0; i < NSTREAM; i++){
    int ioffset = i * iElem;
    CHECK(cudaMemcpyAsync(&d_A[ioffset], &h_A[ioffset], iBytes,
                          cudaMemcpyHostToDevice, stream[i]));
    CHECK(cudaMemcpyAsync(&d_B[ioffset], &h_B[ioffset], iBytes,
                          cudaMemcpyHostToDevice, stream[i]));
}

for (int i = 0; i < NSTREAM; i++){
    int ioffset = i * iElem;
    sumArrays<<<grid, block, 0, stream[i]>>>(&d_A[ioffset], &d_B[ioffset],
            &d_C[ioffset], iElem);
}

for (int i = 0; i < NSTREAM; i++){
    int ioffset = i * iElem;
    CHECK(cudaMemcpyAsync(&gpuRef[ioffset], &d_C[ioffset], iBytes,
                          cudaMemcpyDeviceToHost, stream[i]));
}

3.2 流回调

cuda流中可以插入在主机端运行的函数，这种运行方式称为流回调，和流中的其他kernel或异步API一样，回调函数也是按插入的顺序执行的。回调函数是在主机端执行的，执行时会阻塞流，流回调是一种CPU与GPU同步的方式，功能十分强大，在流中插入（注册）回调函数的方法如下：

cudaError_t cudaStreamAddCallBack(cudaStream_t stream, 
       cudaStreamCallback_t callback, 
       void *userData, 
       unsigned int flags
);

其中callback为主机端函数名，userData为callback中使用的参数，flags是为未来新版本cuda预留的参数，目前还没有意义，默认值为0，不用设置。callback函数签名需要是以下格式：

void CUDART_CB my_callback(cudaStream_t stream, cudaError_t status, void *data);

通过以下的方式可以在流中添加（注册）回调函数：

for (int i = 0; i < n_streams; i++){
    stream_ids[i] = i;
    kernel_1<<<grid, block, 0, streams[i]>>>();
    kernel_2<<<grid, block, 0, streams[i]>>>();
    kernel_3<<<grid, block, 0, streams[i]>>>();
    kernel_4<<<grid, block, 0, streams[i]>>>();
    CHECK(cudaStreamAddCallback(streams[i], my_callback,
                (void *)(stream_ids + i), 0));
}

循环中为每个流配置了4个kernel和1个回调函数。

4. cuda库

对于绝大多数cuda开发者来说，手写kernel功底再深，所实现的通用算子性能也只能无限逼近于cuda官方库函数，官方库如cuBLAS实现的算子已经将一些通用矩阵-向量运算性能优化到了极致。在日常开发中多使用cuda标准库不仅可以减少开发量，还可以明显的提升程序的性能。常用的cuda官方库有几大类：偏底层的用于数学计算的Math Libraries、并行算法的Parallel Algorithm Libraries、深度学习的Deep Learning Libraries、图像视频的Image and Video Libraries等。本文重在cuda基础知识的梳理，不会对所有库的原理和功能全部进行梳理，重点介绍2个常用的Math Libraries： cuBLAS, cuTensor。 未来会开新的专题对深度学习库如：cuDNN, TensorRT等进行详细整理，如果想研究其他cuda库的原理和功能，可以去nvidia官网了解。

Math Libraries几个库的用途如下:
cuBLAS: 稠密矩阵运算的线性代数库，BLAS库的cuda实现；
cuTensor: 张量运算线性代数库。

和其他c/c++第三方库一样，cuda库在使用前需要进行链接，可以使用nvcc -l命令(如-lcublas)进行链接。如果使用cmake构建工具，在CMakeLists.txt文件add_executable命令下方加入target_link_libraries命令，如：

target_link_libraries(cuDemos cublas.lib)

其中cuDemos为项目名。

4.1 cuda库的工作流程

cuda库函数在主机端进行调用，涉及内存分配、数据传输、数据计算几个方面，通用的cuda函数的工作流如下：

1. 创建一个句柄管理运行环境，所谓句柄有点类似于python或java框架下中的Session，如SparkSession、tensorflow Session；
2. 分配设备内存；
3. 如果输入数据类型不是cuda库支持的数据类型，需要进行数据类型转化；
4. 数据传入设备内存；
5. 配置运算条件，并执行计算;
6. 计算结果传回主机内存；
7. 释放cuda资源；

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

#include <cublas_v2.h>
#include "../include/utils.h"

#define M 6
#define N 5
#define IDX2C(i,j,ld) (((j)*(ld))+(i))

#define CHECK_CUBLAS(call){                                                    \
    cublasStatus_t err=call;                                                   \
    if (err != CUBLAS_STATUS_SUCCESS){                                         \
        fprintf(stderr, "Got CUBLAS error %d at %s:%d\n", err, __FILE__,       \
                __LINE__);                                                     \
        exit(1);                                                               \
    }                                                                          \
}

static __inline__ void modify (cublasHandle_t handle, float *m, int ldm, int n, int p, int q, float alpha, float beta){
    cublasSscal (handle, n-q, &alpha, &m[IDX2C(p,q,ldm)], ldm);
    cublasSscal (handle, ldm-p, &beta, &m[IDX2C(p,q,ldm)], 1);
}

void cuBlasProcess(){
    cublasHandle_t handle;
    float* devPtrA;
    float* a = 0;

    a = (float *)malloc (M * N * sizeof (*a));
    CHECK(cudaMalloc ((void**)&devPtrA, M*N*sizeof(*a)));
//    CHECK(cudaMallocManaged((void**)&devPtrA, M*N*sizeof(*a)));

    for (int j = 0; j < N; j++) {
        for (int i = 0; i < M; i++) {
            a[IDX2C(i,j,M)] = (float)(i * N + j + 1);
            devPtrA[IDX2C(i,j,M)] = (float)(i * N + j + 1);
        }
    }

    CHECK_CUBLAS(cublasCreate(&handle));
    CHECK_CUBLAS(cublasSetMatrix (M, N, sizeof(*a), a, M, devPtrA, M));

    modify (handle, devPtrA, M, N, 1, 2, 16.0f, 12.0f);

    CHECK_CUBLAS(cublasGetMatrix (M, N, sizeof(*a), devPtrA, M, a, M));
    cudaFree (devPtrA);

    CHECK_CUBLAS(cublasDestroy(handle));
    for (int j = 0; j < N; j++) {
        for (int i = 0; i < M; i++) {
            printf ("%7.0f", a[IDX2C(i,j,M)]);
        }
        printf ("\n");
    }
    free(a);
}

上面主机端代码演示了通用库cublas的工作流，内联函数modify中使用cublasSscal对矩阵进行了操作，具体的运算逻辑我们先暂时不用管。相比于手写kernel主机端调用代码，上述代码有两处值得注意的地方：

1. 数据传输部分使用cublasSetVector、cublasSetMatrix代替了cudaMemcpy方法，也可以使用第32行注释的cudaMallocManaged方法进行统一内存管理，省去显式的cublasSetVector、cublasSetMatrix；
2. 运算部分没有显式的定义grid, block等信息，这也正是cuda库优势之一，cuda库通常会在内部隐式的定义最优的grid和block，省去人工grid, block调优过程。

4.2 常用cuda库API介绍

cublas库是传统线性代数BLAS库的cuda实现，cublas目前有3套API：cuBLAS,cuBLASXt和cuBLASLt。其中cuBLAS最原始，需要手动的在主机代码中实现host、device端的内存管理和数据传输；cuBLASXt使用内统一内存管理，省去了主机端显式的内存管理API的调用；cuBLASLt是一款更灵活的轻量级库，3套API的头文件如下：

#include <cublas.h> // 旧版 cuBLAS
#include <cublas_v2.h> //新版 cuBLAS
//#include <cublasXt.h>
//#include <cublasLt.h>

新版本cuBLAS库使用cublas_v2.h头文件，在本文中将主要介绍cuBLAS新版本的API，如果想研究另外2套API可以去官网了解。

由于cuBLAS库最早是Fortran语言编写，Fortran语言下数组内存布局采用的是列主序，索引起始值为1，而C/C++数组布局采用的是行主序，索引起始值为0。M行N列的数组，行列索引(i, j)对应全局索引：

行主序：N*i+j;
列主序：M*j+i;

在对主机内存中数组的元素进行写入和读取时，本文将采用以下IDX2C实现列主序、0初始索引：

#define IDX2C(i,j,ld) (((j)*(ld))+(i))

cuBLAS

cuBLAS API返回值

所有的cuBLAS API都返回cublasStatus_t，官网介绍了cuBLAS API在返回不同的cublasStatus_t枚举值时代表的运行状况和解决办法：

表4.1 Nvidia官网上不同cublasStatus_t枚举值介绍

cuBLAS 运行参数配置

cublasSetPointerMode

使用cublasSetPointerMode(cublasHandle_t handle, cublasPointerMode_t mode)函数用来设置标量形指针变量（通常为记录cublas函数运算的结果的指针）传入cuBLAS API的形式，cublasPointerMode_t有两个枚举值：

CUBLAS_POINTER_MODE_HOST: 以主机端指针的形式传递(默认模式);
CUBLAS_POINTER_MODE_DEVICE: 以设备端指针的形式传递;

下面的例子将演示cublasSetPointerMode的使用。

#define M   1024
#define N   1024
#define P   1024

CHECK_CUBLAS(cublasSetPointerMode(cublas_handle, CUBLAS_POINTER_MODE_DEVICE));
for (i = 0; i < M; i++){
    CHECK_CUBLAS(cublasSasum(cublas_handle, P, d_C + (i * P), 1, d_row_sums + i));
}

/*
 * Set cuBLAS back to host pointer mode, indicating that all scalars are
 * passed as host pointers.
 */
CHECK_CUBLAS(cublasSetPointerMode(cublas_handle, CUBLAS_POINTER_MODE_HOST));
/*
 * Do the final sum of the sum of all rows to produce a total for the whole
 * output matrix.
 */
CHECK_CUBLAS(cublasSasum(cublas_handle, M, d_row_sums, 1, &total_sum));

上述代码片段中使用了2次cublasSasum函数，cublasSasum函数签名如下：

cublasStatus_t  cublasSasum(cublasHandle_t handle, int n,
                            const float           *x, int incx, float  *result)

最后一个参数result会记录函数的输出，既可以是device指针，也可以是host指针：

表4.2 官网上cublasSasum函数参数说明

所以需要与cublasSetPointerMode函数联用，16行第1次调用时在循环中对数组d_c的每行进行求和，将每行之和写入到d_row_sums中，传入result处参数为d_row_sums，是一个device端指针，在调用之前的12行需要将pointer mode设置为CUBLAS_POINTER_MODE_DEVICE。28行第2次调用时对每行之和进行第二次求和，result处传入的指针变量为&total_sum，是一个host端指针，需要在调用之前将pointer mode设置为CUBLAS_POINTER_MODE_HOST。
使用cublasGetPointerMode(cublasHandle_t handle, cublasPointerMode_t *mode)可以获取cublas当前的pointer mode并将值写入*mode中。

cublasSetAtomicsMode

使用cublasStatus_t cublasSetAtomicsMode(cublasHandlet handle, cublasAtomicsMode_t mode)设置cuBLAS函数在计算过程中是否可以使用原子操作，原子操作会导致线程一个个排队去执行，降低kernel运行性能，非万不得已尽量不要使用。cublasAtomicsMode_t有两个枚举值：

CUBLAS_ATOMICS_NOT_ALLOWED: cuBLAS函数中禁用原子操作，默认模式
CUBLAS_ATOMICS_ALLOWED: cuBLAS函数中启用原子操作

cublasSetAtomicsMode同样存在镜像方法cublasGetAtomicsMode，用来获取当前的atomics mode。

cublasSetMathMode

使用cublasStatus_t cublasSetMathMode(cublasHandle_t handle, cublasMath_t mode)函数可以设置cuBLAS函数进行数学运算的计算精度。cublasMath_t有以下枚举值：

CUBLAS_DEFAULT_MATH：精度与输入数据一致，这种mode性能最好，默认模式；
CUBLAS_PEDANTIC_MATH：使用规定的精度和标准化算法；
CUBLAS_TF32_TENSOR_OP_MATH: 使用TF32 tensor cores实现单精度加速；
CUBLAS_MATH_DISALLOW_REDUCED_PRECISION_REDUCTION: 计算过程中向最高精度对齐；
CUBLAS_TENSOR_OP_MATH: 已弃用的模式。

cublasSetWorkspace

使用cublasStatus_t cublasSetWorkspace(cublasHandle_t handle, void *workspace, size_t workspaceSizeInBytes)为cuBLAS库函数开辟单独的工作空间，工作空间的开辟是以流为单位的，如果不单独的进行设置，handle下的cuBLAS库函数和进程中的其他的核函数共享同一个默认的工作空间；workspaceSizeInBytes最少要设置256 bytes，推荐的设置大小如下：

表4.3 官网上推荐的workspaceSize

workSpace如果设置不成功，cublasSetStream函数会对cublas handle工作空间重置回默认工作空间。

cublasSetStream

使用cublasStatus_t cublasSetStream(cublasHandle_t handle, cudaStream_t streamId)为cublas handle分配cuda流；如果不进行设置，handle对应cublas函数将在默认的空流下运行

cublasLoggerConfigure

cublasStatus_t cublasLoggerConfigure(
    int             logIsOn,
    int             logToStdOut,
    int             logToStdErr,
    const char*     logFileName)

用于配置cuBLAS库运行logger, 参数logIsOn控制日志是否打开的总开关，默认关闭（0）,logToStdOut为标准输出I/O流的开关，默认关闭；logToStdErr为标准error I/0控制开关，默认关闭；logFileName设置打印日志文件名。

cublasGetLoggerCallback

cublasStatus_t cublasGetLoggerCallback(cublasLogCallback* userCallback)用来设置日志回调udf。

cuBLAS 标量-向量-矩阵运算

cuBLAS运算API一共分为3级，level-1 function用于标量-向量运算，level-2 function用于向量-矩阵运算，level-3 function用于矩阵-矩阵运算。无论哪一级运算API，都至少有4套函数分别用于float, double, float complex(复数), double complex 4种不同数据类型。通过在函数名中间设置4种不同字母实现模板函数的功能：cublas<t>asum()中，函数名中的"<t>"被下表中不同字母（S、D、C、Z）代替表示用于不同数据类型输入值。

表4.4 官网上不同type字母与数据类型对照表

cuBlas function 具有如下的通用规则：

1. 函数命名：cublas前缀+数据类型缩写（S, D, C, Z）+功能(gemm, asum, imax, gemv...)；
2. 参数列表中靠前位置参数为cublasHandle_t类型句柄、cublasOperation_t类型的运算前矩阵操作，cublasOperation_t共三个枚举值：CUBLAS_OP_N表示无操作直接用原始矩阵进行运算，CUBLAS_OP_T表示运算前先对矩阵进行转置，CUBLAS_OP_C表示运算前先求原始矩阵的共轭矩阵；
3. 中间位置参数通常为向量、矩阵数据的长度，通常使用n、m、k表示；
4. 尾部的参数为标量、向量、矩阵指针以及主维度大小（行），其中标量用小写希腊字母α、β表示，向量用小写字母x、y、z表示，大写字母A、B、C表示矩阵，主维度大小用ld表示，如lda表示A矩阵主维度，ldb表示B矩阵主维度；

Level-1 function

下面将整理一部分cublas L1级别的运算API，包括cublasI<t>amax()、cublas<t>asum()、cublas<t>dot()、cublas<t>nrm2()、cublas<t>scal()、cublas<t>swap()；
其他Level-1 function可以去Nvidia官网查看。

cublasI<t>amax()

cublasStatus_t cublasIsamax(cublasHandle_t handle, int n,
                            const float *x, int incx, int *result)
cublasStatus_t cublasIdamax(cublasHandle_t handle, int n,
                            const double *x, int incx, int *result)
cublasStatus_t cublasIcamax(cublasHandle_t handle, int n,
                            const cuComplex *x, int incx, int *result)
cublasStatus_t cublasIzamax(cublasHandle_t handle, int n,
                            const cuDoubleComplex *x, int incx, int *result)

在x中每incx个元素取一个元素，从取出的元素中找出最大元素对应的索引（1-based），最大元素不唯一取第一个，将结果返回result中：
x: [2.0,4,5,1,2,4,1,5], incx: 3
每3个元素取出1个：[2,1,1]
[2,1,1]中第一个元素最大，result: 1
如果incx设为1，则返回x中最大值元素对应的第一个索引。

cublas<t>asum()

cublasStatus_t  cublasSasum(cublasHandle_t handle, int n,
                            const float           *x, int incx, float  *result)
cublasStatus_t  cublasDasum(cublasHandle_t handle, int n,
                            const double          *x, int incx, double *result)
cublasStatus_t cublasScasum(cublasHandle_t handle, int n,
                            const cuComplex       *x, int incx, float  *result)
cublasStatus_t cublasDzasum(cublasHandle_t handle, int n,
                            const cuDoubleComplex *x, int incx, double *result)

在x中每隔incx个元素取一个元素，求和，并将结果返回result中。

cublas<t>dot()

cublasStatus_t cublasSdot (cublasHandle_t handle, int n,
                           const float           *x, int incx,
                           const float           *y, int incy,
                           float           *result)
cublasStatus_t cublasDdot (cublasHandle_t handle, int n,
                           const double          *x, int incx,
                           const double          *y, int incy,
                           double          *result)
cublasStatus_t cublasCdotu(cublasHandle_t handle, int n,
                           const cuComplex       *x, int incx,
                           const cuComplex       *y, int incy,
                           cuComplex       *result)
cublasStatus_t cublasCdotc(cublasHandle_t handle, int n,
                           const cuComplex       *x, int incx,
                           const cuComplex       *y, int incy,
                           cuComplex       *result)
cublasStatus_t cublasZdotu(cublasHandle_t handle, int n,
                           const cuDoubleComplex *x, int incx,
                           const cuDoubleComplex *y, int incy,
                           cuDoubleComplex *result)
cublasStatus_t cublasZdotc(cublasHandle_t handle, int n,
                           const cuDoubleComplex *x, int incx,
                           const cuDoubleComplex *y, int incy,
                           cuDoubleComplex       *result)

在x中每incx个元素取一个元素，在y中每隔incy个元素取一个元素:
x = {1.0,2,3,4,5,6,7,8}, y = {9.0, 10,11,12,13,14,15,16}, incx=2, incy=4;
取出:{1.0, 3, 5, 7}, {9, 13}
运算后：
result = 48 = (1*9+3*13)
此组API中复数部分又分出了以c和u为后缀的2类，计算复数内积时先要在x或y向量中中选一个，对其中元素取共轭，然后再执行乘累加，以c结尾的API表示乘累加前对x向量元素取共轭。

cublas<t>nrm2()

cublasStatus_t  cublasSnrm2(cublasHandle_t handle, int n,
                            const float           *x, int incx, float  *result)
cublasStatus_t  cublasDnrm2(cublasHandle_t handle, int n,
                            const double          *x, int incx, double *result)
cublasStatus_t cublasScnrm2(cublasHandle_t handle, int n,
                            const cuComplex       *x, int incx, float  *result)
cublasStatus_t cublasDznrm2(cublasHandle_t handle, int n,
                            const cuDoubleComplex *x, int incx, double *result)

从x中每incx个元素取出一个，计算l2范数，结果存入result中。

cublas<t>scal()

cublasStatus_t  cublasSscal(cublasHandle_t handle, int n,
                            const float           *alpha,
                            float           *x, int incx)
cublasStatus_t  cublasDscal(cublasHandle_t handle, int n,
                            const double          *alpha,
                            double          *x, int incx)
cublasStatus_t  cublasCscal(cublasHandle_t handle, int n,
                            const cuComplex       *alpha,
                            cuComplex       *x, int incx)
cublasStatus_t cublasCsscal(cublasHandle_t handle, int n,
                            const float           *alpha,
                            cuComplex       *x, int incx)
cublasStatus_t  cublasZscal(cublasHandle_t handle, int n,
                            const cuDoubleComplex *alpha,
                            cuDoubleComplex *x, int incx)
cublasStatus_t cublasZdscal(cublasHandle_t handle, int n,
                            const double          *alpha,
                            cuDoubleComplex *x, int incx)

将x中每incx个元素乘上alpha，复数API中多分出了2个，分别用于alpha是float、double复数的情形。

cublas<t>swap()

cublasStatus_t cublasSswap(cublasHandle_t handle, int n, float           *x,
                           int incx, float           *y, int incy)
cublasStatus_t cublasDswap(cublasHandle_t handle, int n, double          *x,
                           int incx, double          *y, int incy)
cublasStatus_t cublasCswap(cublasHandle_t handle, int n, cuComplex       *x,
                           int incx, cuComplex       *y, int incy)
cublasStatus_t cublasZswap(cublasHandle_t handle, int n, cuDoubleComplex *x,
                           int incx, cuDoubleComplex *y, int incy)

从x中每incx个元素取出一个，从y中每incy个元素取出一个，进行互换。

Level-2 function

L2级别用于m-v运算cublas函数将介绍cublas<t>gbmv(), cublas<t>gbmv(), cublas<t>trsv，其他函数用法都大同小异。

cublas<t>gemv()

cublasStatus_t cublasSgemv(cublasHandle_t handle, cublasOperation_t trans,
                           int m, int n,
                           const float           *alpha,
                           const float           *A, int lda,
                           const float           *x, int incx,
                           const float           *beta,
                           float           *y, int incy)
cublasStatus_t cublasDgemv(cublasHandle_t handle, cublasOperation_t trans,
                           int m, int n,
                           const double          *alpha,
                           const double          *A, int lda,
                           const double          *x, int incx,
                           const double          *beta,
                           double          *y, int incy)
cublasStatus_t cublasCgemv(cublasHandle_t handle, cublasOperation_t trans,
                           int m, int n,
                           const cuComplex       *alpha,
                           const cuComplex       *A, int lda,
                           const cuComplex       *x, int incx,
                           const cuComplex       *beta,
                           cuComplex       *y, int incy)
cublasStatus_t cublasZgemv(cublasHandle_t handle, cublasOperation_t trans,
                           int m, int n,
                           const cuDoubleComplex *alpha,
                           const cuDoubleComplex *A, int lda,
                           const cuDoubleComplex *x, int incx,
                           const cuDoubleComplex *beta,
                           cuDoubleComplex *y, int incy)

用于通用矩阵(m)-向量(v)乘法运算，通用矩阵-向量乘法运算格式： $y = αop(A)x + βy$ 其中op(A)通过参数cublasOperation_t trans控制，不同枚举值对应操作如下：

API要点如下：

1. 如果trans为转置或共轭，x,y向量的维度要随发生互换；
2. 如果incx, incy > 1, 最好保证x,y的元素够多，保证采样后的长度仍能满足计算要求，如果采样后长度不够，API或自动用0对齐之后再进行运算。 如矩阵A、x、y取值如下：

x = {1.0, 2, 3, 4};
y = {5.0, 6 ,7 ,8};
A = {1, 3, 0, 0,
     2, 1, 3, 0,
     0, 2, 1, 3,
     0, 0, 2, 1};

如果alpha = 1, beta = 1, incx = 1, incy = 1，运算后y = 1*{7, 13, 19, 10} + 1*{5, 6, 7, 8} = {12.0, 19, 26, 18)
如果incx = 2, incy = 1, 运算后y = 1*A*{1,3,0,0} + 1*{5, 6, 7, 8} = {15,11,13,8};
如果incx = 2, incy = 2, 运算后y = 1*A*{1, 3} + 1*{5, 7} = {15, 12}, 然后用15, 12 替换y中的{5, 7}, 最终y = {15, 6, 12, 8}

cublas<t>gbmv()

cublasStatus_t cublasSgbmv(cublasHandle_t handle, cublasOperation_t trans,
                           int m, int n, int kl, int ku,
                           const float           *alpha,
                           const float           *A, int lda,
                           const float           *x, int incx,
                           const float           *beta,
                           float           *y, int incy)
cublasStatus_t cublasDgbmv(cublasHandle_t handle, cublasOperation_t trans,
                           int m, int n, int kl, int ku,
                           const double          *alpha,
                           const double          *A, int lda,
                           const double          *x, int incx,
                           const double          *beta,
                           double          *y, int incy)
cublasStatus_t cublasCgbmv(cublasHandle_t handle, cublasOperation_t trans,
                           int m, int n, int kl, int ku,
                           const cuComplex       *alpha,
                           const cuComplex       *A, int lda,
                           const cuComplex       *x, int incx,
                           const cuComplex       *beta,
                           cuComplex       *y, int incy)
cublasStatus_t cublasZgbmv(cublasHandle_t handle, cublasOperation_t trans,
                           int m, int n, int kl, int ku,
                           const cuDoubleComplex *alpha,
                           const cuDoubleComplex *A, int lda,
                           const cuDoubleComplex *x, int incx,
                           const cuDoubleComplex *beta,
                           cuDoubleComplex *y, int incy)

用于通用的带状矩阵(m)-向量(v)乘法运算，带状矩阵是一种特殊的稀疏矩阵，非零元素分布在对角线。 参数kl、ku表示下带宽和上带宽，下面矩阵A的kl = 1, ku = 2；注意此API中输入矩阵A为被压缩后带状矩阵。压缩方式如下：
A(i, j) 元素放入 A(ku+1+i-j,j)中，其中j=1, ...n; i∈[max(1, j-ku), min(m, j+kl)] (1-based)

// 压缩前
A = {1, 3, 0, 0,
     2, 1, 3, 0, 
     0, 2, 1, 3,
     0, 0, 2, 1}
     
// 压缩后
Ab = {0, 3, 3, 3
      1, 1, 1, 1
      2, 2, 2, 0}

输入API的矩阵为列优先存储Ab, m,n不变, 此时lda为3，使用示例如下：

void cuBlasGemvProcess(){
    cublasHandle_t handle;
    int m = 4, n = 4;
    int kl = 1, ku = 1;

    float *x ,*y, *A, *Ab;


    float hx[] = {1,2,3,4};
    float hy[] = {5, 6,7,8};

    // 主机端显式定义列优先存储的数组
    
    float hA[] = {1,2,0,0,
                3,1,2,0,
                0,3,1,2,
                0,0,3,1};

    float hAb[] = {0,1,2,
                  3,1,2,
                  3,1,2,
                  3,1,0};

    int nAb = sizeof(hAb)/sizeof(hAb[0]);

    CHECK_CUBLAS(cublasCreate(&handle));
    CHECK(cudaMallocManaged((void**)&x, n*sizeof(float)));
    CHECK(cudaMallocManaged((void**)&y, m*sizeof(float)));
    CHECK(cudaMallocManaged((void**)&A, m*n*sizeof(float)));
    CHECK(cudaMallocManaged((void**)&Ab, nAb*sizeof(float)));

    // 数据拷贝到设备端
    for(int i=0; i<std::max(m,n); i++){
        if(i<n){
            x[i] = hx[i];
        }
        if(i<m){
            y[i] = hy[i];
        }
    }

    for(int j=0; j<n; j++){
        for(int i=0; i<m; i++){
            A[IDX2C(i,j,m)] = hA[IDX2C(i,j,m)];
        }
    }

    for(int i=0; i<nAb; i++){
        Ab[i] = hAb[i];
    }

    int incx = 1, incy = 1;
    float alpha = 1;
    float beta = 0;
    
    // 运算
    CHECK_CUBLAS(cublasSgbmv(handle,CUBLAS_OP_N,m,n,kl,ku,&alpha,Ab,3,x,incx,&beta,y,incy));
//    CHECK_CUBLAS(cublasSgemv(handle,CUBLAS_OP_N,m,n,&alpha,A,m,x,incx,&beta,y,incy));
    
    // 结果写回主机内存
    float *ry = (float*)malloc(m*sizeof(float));
    CHECK_CUBLAS(cublasGetVector(m,sizeof(float), y, 1, ry, 1));
    
    // 打印结果
    printf("ry:\n");
    for(int i=0; i<m; i++){
        printf("ry[%d]:%f\n", i,ry[i]);
    }
    free(ry);

//    printf("result: %.3f",result);

    cudaFree(x);
    cudaFree(y);
    cudaFree(A);
//
    CHECK_CUBLAS(cublasDestroy(handle));
}

cublas<t>trsv

cublasStatus_t cublasStrsv(cublasHandle_t handle, cublasFillMode_t uplo,
                           cublasOperation_t trans, cublasDiagType_t diag,
                           int n, const float           *A, int lda,
                           float           *x, int incx)
cublasStatus_t cublasDtrsv(cublasHandle_t handle, cublasFillMode_t uplo,
                           cublasOperation_t trans, cublasDiagType_t diag,
                           int n, const double          *A, int lda,
                           double          *x, int incx)
cublasStatus_t cublasCtrsv(cublasHandle_t handle, cublasFillMode_t uplo,
                           cublasOperation_t trans, cublasDiagType_t diag,
                           int n, const cuComplex       *A, int lda,
                           cuComplex       *x, int incx)
cublasStatus_t cublasZtrsv(cublasHandle_t handle, cublasFillMode_t uplo,
                           cublasOperation_t trans, cublasDiagType_t diag,
                           int n, const cuDoubleComplex *A, int lda,
                           cuDoubleComplex *x, int incx)

用于求矩阵线性代数方程$op(A)x = b$，其中A为三角矩阵，参数uplo控制A为上三角还是下三角，输入时x为b，输出时用得到的解覆盖x，diag控制是否将A的对角线元素视为1，API不会校验A是否为奇异矩阵。

Level-3 function

L3级别cublas函数用于m-m运算，将介绍cublas<t>gemm(),

cublas<t>gemm()

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)
cublasStatus_t cublasDgemm(cublasHandle_t handle,
                           cublasOperation_t transa, cublasOperation_t transb,
                           int m, int n, int k,
                           const double          *alpha,
                           const double          *A, int lda,
                           const double          *B, int ldb,
                           const double          *beta,
                           double          *C, int ldc)
cublasStatus_t cublasCgemm(cublasHandle_t handle,
                           cublasOperation_t transa, cublasOperation_t transb,
                           int m, int n, int k,
                           const cuComplex       *alpha,
                           const cuComplex       *A, int lda,
                           const cuComplex       *B, int ldb,
                           const cuComplex       *beta,
                           cuComplex       *C, int ldc)
cublasStatus_t cublasZgemm(cublasHandle_t handle,
                           cublasOperation_t transa, cublasOperation_t transb,
                           int m, int n, int k,
                           const cuDoubleComplex *alpha,
                           const cuDoubleComplex *A, int lda,
                           const cuDoubleComplex *B, int ldb,
                           const cuDoubleComplex *beta,
                           cuDoubleComplex *C, int ldc)
cublasStatus_t cublasHgemm(cublasHandle_t handle,
                           cublasOperation_t transa, cublasOperation_t transb,
                           int m, int n, int k,
                           const __half *alpha,
                           const __half *A, int lda,
                           const __half *B, int ldb,
                           const __half *beta,
                           __half *C, int ldc)

用于通用矩阵乘法运算$C=αop(A)op(B)+βC$，其中多出的一个函数cublasHgemm表示输入矩阵数据类型为半精(fp 16)。API中有m, n, k三个参数，m: A和C的行数; n: B和C的列数; k :A的列数和B的行数。运算完的结果写入C中，对C进行覆盖。使用示例如下

void cuBlasGemmProcess(){
    cublasHandle_t handle;
    int m = 5, n = 4, k = 3;

    float *A, *B, *C;

    // 主机端显式定义列优先存储的数组
    // A: 5x3 rows = 5, cols = 3
    float hA[] = {7,4,2,8, 2,
                  1,3,3,6,-4,
                  2,-2,-5,3,1};
    // B: 3x4 rows = 4, cols = 4
    float hB[] = {2,0,1,
                  5,1,0,
                  -3,-2,2,
                  1,2,-1};
    // C: 5x4 rows = 5, cols =4
    float hC[] = {2,0,4,2,1,
                  5,1,0,0,-2,
                  -3,-2,2,0,5,
                  1,2,-1,2,-7};

    CHECK_CUBLAS(cublasCreate(&handle));
    CHECK(cudaMallocManaged((void**)&A, m*k*sizeof(float)));
    CHECK(cudaMallocManaged((void**)&B, k*n*sizeof(float)));
    CHECK(cudaMallocManaged((void**)&C, m*n*sizeof(float)));

    // 数据拷贝到设备端
    for(int j=0; j<k; j++){
        for(int i=0; i<m; i++){
            A[IDX2C(i,j,m)] = hA[IDX2C(i,j,m)];
        }
    }

    for(int j=0; j<n; j++){
        for(int i=0; i<k; i++){
            B[IDX2C(i,j,k)] = hB[IDX2C(i,j,k)];
        }
    }

    for(int j=0; j<n; j++){
        for(int i=0; i<m; i++){
            C[IDX2C(i,j,m)] = hC[IDX2C(i,j,m)];
        }
    }

    float alpha = 1, beta = 1;

    // 运算
    CHECK_CUBLAS(cublasSgemm(handle, CUBLAS_OP_N,CUBLAS_OP_N,m,n,k,&alpha,A,m,B,k,&beta,C,m));

    // 数据从设备端拷贝回主机端
    float *rC = (float*)malloc(m*n*sizeof(float));
    CHECK_CUBLAS(cublasGetMatrix(m,n,sizeof(float ),C,m,rC,m));

    // 列优先存储的数组转置成行优先方式打印出来
    printf("rC:\n");
    for (int i = 0; i < m; i++) {
        for(int j=0; j<n; j++) {
            int idx = j*m+i;
            printf("rx[%d][%d]:%f, ", i,j, rC[idx]);
        }
        printf("\n");
    }
    free(rC);

    // 释放内存

    cudaFree(A);
    cudaFree(B);
    cudaFree(C);
//
    CHECK_CUBLAS(cublasDestroy(handle));
}

cublas<t>trsm()

cublasStatus_t cublasStrsm(cublasHandle_t handle,
                           cublasSideMode_t side, cublasFillMode_t uplo,
                           cublasOperation_t trans, cublasDiagType_t diag,
                           int m, int n,
                           const float           *alpha,
                           const float           *A, int lda,
                           float           *B, int ldb)
cublasStatus_t cublasDtrsm(cublasHandle_t handle,
                           cublasSideMode_t side, cublasFillMode_t uplo,
                           cublasOperation_t trans, cublasDiagType_t diag,
                           int m, int n,
                           const double          *alpha,
                           const double          *A, int lda,
                           double          *B, int ldb)
cublasStatus_t cublasCtrsm(cublasHandle_t handle,
                           cublasSideMode_t side, cublasFillMode_t uplo,
                           cublasOperation_t trans, cublasDiagType_t diag,
                           int m, int n,
                           const cuComplex       *alpha,
                           const cuComplex       *A, int lda,
                           cuComplex       *B, int ldb)
cublasStatus_t cublasZtrsm(cublasHandle_t handle,
                           cublasSideMode_t side, cublasFillMode_t uplo,
                           cublasOperation_t trans, cublasDiagType_t diag,
                           int m, int n,
                           const cuDoubleComplex *alpha,
                           const cuDoubleComplex *A, int lda,
                           cuDoubleComplex *B, int ldb)

用于求解矩阵线性代数方程组 $op(A)X = B$ (side=CUBLAS_SIDE_LEFT)或 $Xop(A) = B$ (side=CUBLAS_SIDE_LEFT)，其中A为三角矩阵，参数uplo控制A为上三角还是下三角，输入时x为b，输出时用得到的解覆盖x，diag控制是否将A的对角线元素视为1，API不会校验A是否为奇异矩阵。

cuTensor

cuTensor是nvidia官方开发的用于张量运算的cuda库，API相比于cuBLAS更高一级，主要特性有： 支持维度高（上限64维）、支持混精、任意的数据布局、支持计划缓存。cuTensor的张量和部分运算API在设计上使用了描述符（Descriptor），所谓的描述符其实是一个类，用于被描述对象的属性托管，这样大量属性相关操作就可以用Descriptor进行统一管理。与cuBLAS库一样，cuTensor API同样也需要传入用于管理上下文信息的句柄(cutensorHandle_t)。返回值为API执行状态(cutensorStatus_t)。不同点在于cuTensor中张量运算函数抽象程度更高，不像cuBLAS中一样，各level-function对应到每一类具体的向量、矩阵操作。cuTensor中张量运算函数有3类：元素级别运算（element-wise）、张量间收缩运算、张量规约运算，每种运算由少量抽象程度更高的API实现，在3类API中通过一个枚举类型的参数来控制具体的运算规则，具体方式如下表所示：

表4.5 cuTensor不同运算类型的函数介绍

cuTensor配置

cuTensor不是cuda自带的库，需要单独的进行安装和配置。本文使用的环境为windows 11 + clion专业版 + VS 2022内置的MSVC 1936 + CUDA 11.8 + cutensor 1.7.0 + CMake 3.24。首先需要下载并安装cutensor，安装好之后的目录下有inculde、lib两个文件夹：

将include下的内容copy到cuda下的include中：

lib目录下有4个版本，选择其中一个（本文选择的是11）：

将lib/11目录下的.dll文件copy到cuda的bin目录下：

lib/11目录下的.lib文件copy到cuda的lib/x64目录下： 最后，在项目CMakeLists.txt文件中添加link：

target_link_libraries(项目名 cublas.lib cutensor.lib)

至此完成了cutensor环境配置，windows环境下其他三方库都可以采用这种办法安装。

描述符(Descriptor)

cuTensor中有2种描述符：cutensorTensorDescriptor_t, cutensorConstractionDescriptor_t。其中cutensorTensorDescriptor_t用于描述一个张量的属性，在学习API之前我们需要对cuTensor中张量的属性和命名方式有个了解。

cuTensor中抽象出了mode、extent、stride概念，一个n阶（n维）张量有n个mode，每个mode拥有1个extent（单维长度），stride指同一维度上索引相同的两个元素在连续的物理内存上相距多远。如1 个4x8x16的张量，mode-a, mode-b, mode-c对应的extent-a, extent-b, extent-c分别为4，8，16；stride-a, stride-b, stride-c分别为1, 4, 4x8。这种布局是数据恰好在内存中连续存储的情形，通用情况下：stride(1)>=1, stride(i)>=stride(i-1)*extend(i-1)。

使用cutensorInitTensorDescriptor()可以对张量的cutensorTensorDescriptor_t进行初始化，完整签名如下：

cutensorStatus_t cutensorInitTensorDescriptor(
    const cutensorHandle_t *handle, 
    cutensorTensorDescriptor_t *desc, 
    const uint32_t numModes, 
    const int64_t extent[], const int64_t stride[], 
    cudaDataType_t dataType, cutensorOperator_t unaryOp
)

其中：

cutensorHandle_t *handle: 句柄
cutensorTensorDescriptor_t *desc: 需要进行初始化的Descriptor
numModes: mode数量
int64_t extent[]: extent数组
int64_t stride[]: stride数组，如果设置为NULL，将根据extent按连续存储的方式自动计算
cudaDataType_t dataType: 张量元素的数据类型 
cutensorOperator_t unaryOp: 参与运算前对张量中的元素执行的操作，需要设置为支持一元操作的枚举值

具体的使用方法如下：

// 定义modeSize, extent, datatype;
int nmode = 3;
int64_t extent[] = {96, 32, 16};
cudaDataType_t type = CUDA_R_32F;

cutensorHandle_t* handle;
CHECK_CUTENSOR( cutensorCreate(&handle));

cutensorTensorDescriptor_t desc;
CHECK_CUTENSOR(cutensorInitTensorDescriptor(handle,
                                            &desc,
                                            nmode,
                                            extent,
                                            NULL,/*stride*/
                                            type, CUTENSOR_OP_IDENTITY));

cutensorConstractionDescriptor_t 用于描述参与收缩运算的张量A、B、C分别是怎样的、输出的张量D是怎样的。所以cutensorContractionDescriptor_t初始化方法要稍微复杂一些：

cutensorStatus_t cutensorInitContractionDescriptor(
    const cutensorHandle_t *handle, 
    cutensorContractionDescriptor_t *desc, 
    const cutensorTensorDescriptor_t *descA, 
    const int32_t modeA[], 
    const uint32_t alignmentRequirementA, 
    const cutensorTensorDescriptor_t *descB, 
    const int32_t modeB[], 
    const uint32_t alignmentRequirementB, 
    const cutensorTensorDescriptor_t *descC, 
    const int32_t modeC[], 
    const uint32_t alignmentRequirementC, 
    const cutensorTensorDescriptor_t *descD,
    const int32_t modeD[], 
    const uint32_t alignmentRequirementD,
    cutensorComputeType_t typeCompute
)

cutensorInitContractionDescriptor()中为每个张量设有一组参数，以张量A为例：

const cutensorTensorDescriptor_t *descA, 
const int32_t modeA[], 
const uint32_t alignmentRequirementA, 

其中modeA[ ]长度需要与descA中的设置的extent相对应，值通常为字符，遵循爱因斯坦符号。如以下格式：

int32_t modeA[] = {'m','u','n','v'};

modeA[ ]中字符出现的顺序会影响数据在内存中的布局，但在运算过程中cutensor库会自动进行隐式转置，优化运算性能。参数alignmentRequirementA指定：在张量A中根据descA，满足最小对齐要求最多需要多少个字节。可以通过以下API计算：

cutensorStatus_t cutensorGetAlignmentRequirement(
    const cutensorHandle_t *handle, 
    const void *ptr, 
    const cutensorTensorDescriptor_t *desc, 
    uint32_t *alignmentRequirement
)

参数typeCompute用于指定收缩运算的数据类型。cutensorInitContractionDescriptor()完整的使用示例如下：

void cuTensorContractionInitProcess(){
    typedef float floatTypeA;
    typedef float floatTypeB;
    typedef float floatTypeC;

    // 定义cuda数据类型
    cudaDataType_t typeA = CUDA_R_32F;
    cudaDataType_t typeB = CUDA_R_32F;
    cudaDataType_t typeC = CUDA_R_32F;
    cutensorComputeType_t typeCompute = CUTENSOR_COMPUTE_32F;

    printf("Include headers and define data types\n");


    // 定义mode和extent
    std::vector<int> modeC{'m','u','n','v'};
    std::vector<int> modeA{'m','h','k','n'};
    std::vector<int> modeB{'u','k','v','h'};
    int nmodeA = modeA.size();
    int nmodeB = modeB.size();
    int nmodeC = modeC.size();

    std::unordered_map<int, int64_t> extent;
    extent['m'] = 96;
    extent['n'] = 96;
    extent['u'] = 96;
    extent['v'] = 64;
    extent['h'] = 64;
    extent['k'] = 64;

    std::vector<int64_t> extentA;
    for(auto mode : modeA)
        extentA.push_back(extent[mode]);

    std::vector<int64_t> extentB;
    for(auto mode : modeB)
        extentB.push_back(extent[mode]);

    std::vector<int64_t> extentC;
    for(auto mode : modeC)
        extentC.push_back(extent[mode]);

    printf("Define modes and extents\n");


    // 计算张量A、B、C元素数量和存储所需字节数量
    size_t elementsA = 1;
    for(auto mode : modeA)
        elementsA *= extent[mode];
    size_t elementsB = 1;
    for(auto mode : modeB)
        elementsB *= extent[mode];
    size_t elementsC = 1;
    for(auto mode : modeC)
        elementsC *= extent[mode];

    size_t sizeA = sizeof(floatTypeA) * elementsA;
    size_t sizeB = sizeof(floatTypeB) * elementsB;
    size_t sizeC = sizeof(floatTypeC) * elementsC;

    // 分配统一内存
    float *A_d, *B_d, *C_d;
    cudaMallocManaged((void**)&A_d, sizeA);
    cudaMallocManaged((void**)&B_d, sizeB);
    cudaMallocManaged((void**)&C_d, sizeC);

    // 初始化数据
    for(int64_t i = 0; i < elementsA; i++)
        A_d[i] = (((float) rand())/RAND_MAX - 0.5)*100;
    for(int64_t i = 0; i < elementsB; i++)
        B_d[i] = (((float) rand())/RAND_MAX - 0.5)*100;
    for(int64_t i = 0; i < elementsC; i++)
        C_d[i] = (((float) rand())/RAND_MAX - 0.5)*100;

    printf("Allocate, initialize and transfer tensors\n");


    cutensorHandle_t* handle;
    CHECK_CUTENSOR( cutensorCreate(&handle));

    // 定义并初始化 cutensorTensorDescriptor
    cutensorTensorDescriptor_t descA;
    CHECK_CUTENSOR( cutensorInitTensorDescriptor( handle,
                                                  &descA,
                                                  nmodeA,
                                                  extentA.data(),
                                                  NULL,/*stride*/
                                                  typeA, CUTENSOR_OP_IDENTITY ) );

    cutensorTensorDescriptor_t descB;
    CHECK_CUTENSOR( cutensorInitTensorDescriptor( handle,
                                                  &descB,
                                                  nmodeB,
                                                  extentB.data(),
                                                  NULL,/*stride*/
                                                  typeB, CUTENSOR_OP_IDENTITY ) );

    cutensorTensorDescriptor_t descC;
    CHECK_CUTENSOR( cutensorInitTensorDescriptor( handle,
                                                  &descC,
                                                  nmodeC,
                                                  extentC.data(),
                                                  NULL,/*stride*/
                                                  typeC, CUTENSOR_OP_IDENTITY ) );

    printf("Initialize cuTENSOR and tensor descriptors\n");


    // 计算各个张量的alignmentRequirement
    uint32_t alignmentRequirementA;
    CHECK_CUTENSOR( cutensorGetAlignmentRequirement( handle,
                                                     A_d,
                                                     &descA,
                                                     &alignmentRequirementA) );

    printf("alignmentRequirementA: %d\n", alignmentRequirementA);

    uint32_t alignmentRequirementB;
    CHECK_CUTENSOR( cutensorGetAlignmentRequirement( handle,
                                                     B_d,
                                                     &descB,
                                                     &alignmentRequirementB) );
    printf("alignmentRequirementB: %d\n", alignmentRequirementB);

    uint32_t alignmentRequirementC;
    CHECK_CUTENSOR( cutensorGetAlignmentRequirement( handle,
                                                     C_d,
                                                     &descC,
                                                     &alignmentRequirementC) );

    printf("alignmentRequirementC: %d\n", alignmentRequirementC);

    printf("Query best alignment requirement for our pointers\n");

    // 定义并初始化 cutensorContractionDescriptor_t
    cutensorContractionDescriptor_t desc;
    CHECK_CUTENSOR( cutensorInitContractionDescriptor( handle,
                                                       &desc,
                                                       &descA, modeA.data(), alignmentRequirementA,
                                                       &descB, modeB.data(), alignmentRequirementB,
                                                       &descC, modeC.data(), alignmentRequirementC,
                                                       &descC, modeC.data(), alignmentRequirementC, // 覆盖C作为输出
                                                       typeCompute));
}

tensor运算

前面介绍了cutensor中张量运算分为：元素级别运算、张量间收缩运算以及张量规约运算，其中张量间的收缩运算相对来讲更复杂一些，在运算前需要定义descriptor、plan、workspace等前置操作。在上一个小节中介绍了如何定义descriptor，在此基础上我们首先介绍张量间的收缩运算，使用cutensorContraction()函数可以实现通用收缩运算：

 cutensorStatus_t cutensorContraction(
    const cutensorHandle_t *handle, 
    const cutensorContractionPlan_t *plan, 
    const void *alpha, 
    const void *A, 
    const void *B, 
    const void *beta, 
    const void *C, 
    void *D, 
    void *workspace, 
    uint64_t workspaceSize, 
    cudaStream_t stream
)

其中cutensorContractionPlan_t *plan参数为最核心参数，使用cutensorInitContractionPlan()定义：

cutensorStatus_t cutensorInitContractionPlan(
    const cutensorHandle_t *handle, 
    cutensorContractionPlan_t *plan, 
    const cutensorContractionDescriptor_t *desc, 
    const cutensorContractionFind_t *find, 
    const uint64_t workspaceSize
)

其中desc参数采用上一小节中已介绍的方法定义，find参数使用cutensorInitContractionFind定义：

cutensorStatus_t cutensorInitContractionFind(
    const cutensorHandle_t *handle, 
    cutensorContractionFind_t *find, 
    const cutensorAlgo_t algo
)

此外cutensorContraction()还有两个关键参数workspace和workspaceSize用来设置运算中优化空间（设备内存）和空间大小。通常先根据desc和find计算得到workspaceSize，然后分配workspace的设备内存：

size_t worksize = 0;
CHECK_CUTENSOR( cutensorContractionGetWorkspaceSize(handle,
                                              &desc,
                                              &find,
                                              CUTENSOR_WORKSPACE_RECOMMENDED,
                                              &worksize ))

void *work = nullptr;
if(worksize > 0){
    if( cudaSuccess != cudaMalloc(&work, worksize)){
        work = nullptr;
        worksize = 0;
    }
}

workspace必须对齐128字节。综上，使用cutensorContraction()进行张量间的收缩运算的流程为：

1. 准备数据并定义descriptor;
2. 定义find;
3. 通过find和descriptor计算workspaceSize并定义workSpace和plan;
4. 传入plan参数执行运算。

完整的使用示例如下：

void cuTensorContractionProcess(){
    typedef float floatTypeA;
    typedef float floatTypeB;
    typedef float floatTypeC;

    // 定义cuda数据类型
    cudaDataType_t typeA = CUDA_R_32F;
    cudaDataType_t typeB = CUDA_R_32F;
    cudaDataType_t typeC = CUDA_R_32F;
    cutensorComputeType_t typeCompute = CUTENSOR_COMPUTE_32F;

    printf("Include headers and define data types\n");


    // 定义mode和extent
    std::vector<int> modeC{'m','u','n','v'};
    std::vector<int> modeA{'m','h','k','n'};
    std::vector<int> modeB{'u','k','v','h'};
    int nmodeA = modeA.size();
    int nmodeB = modeB.size();
    int nmodeC = modeC.size();

    std::unordered_map<int, int64_t> extent;
    extent['m'] = 96;
    extent['n'] = 96;
    extent['u'] = 96;
    extent['v'] = 64;
    extent['h'] = 64;
    extent['k'] = 64;

    std::vector<int64_t> extentA;
    for(auto mode : modeA)
        extentA.push_back(extent[mode]);

    std::vector<int64_t> extentB;
    for(auto mode : modeB)
        extentB.push_back(extent[mode]);

    std::vector<int64_t> extentC;
    for(auto mode : modeC)
        extentC.push_back(extent[mode]);

    printf("Define modes and extents\n");


    // 计算张量A、B、C元素数量和存储所需字节数量
    size_t elementsA = 1;
    for(auto mode : modeA)
        elementsA *= extent[mode];
    size_t elementsB = 1;
    for(auto mode : modeB)
        elementsB *= extent[mode];
    size_t elementsC = 1;
    for(auto mode : modeC)
        elementsC *= extent[mode];

    size_t sizeA = sizeof(floatTypeA) * elementsA;
    size_t sizeB = sizeof(floatTypeB) * elementsB;
    size_t sizeC = sizeof(floatTypeC) * elementsC;

    // 分配统一内存
    float *A_d, *B_d, *C_d;
    cudaMallocManaged((void**)&A_d, sizeA);
    cudaMallocManaged((void**)&B_d, sizeB);
    cudaMallocManaged((void**)&C_d, sizeC);

    // 初始化数据
    for(int64_t i = 0; i < elementsA; i++)
        A_d[i] = (((float) rand())/RAND_MAX - 0.5)*100;
    for(int64_t i = 0; i < elementsB; i++)
        B_d[i] = (((float) rand())/RAND_MAX - 0.5)*100;
    for(int64_t i = 0; i < elementsC; i++)
        C_d[i] = (((float) rand())/RAND_MAX - 0.5)*100;

    printf("Allocate, initialize and transfer tensors\n");


    cutensorHandle_t* handle;
    CHECK_CUTENSOR( cutensorCreate(&handle));

    // 定义并初始化 cutensorTensorDescriptor
    cutensorTensorDescriptor_t descA;
    CHECK_CUTENSOR( cutensorInitTensorDescriptor( handle,
                                                  &descA,
                                                  nmodeA,
                                                  extentA.data(),
                                                  NULL,/*stride*/
                                                  typeA, CUTENSOR_OP_IDENTITY ) );

    cutensorTensorDescriptor_t descB;
    CHECK_CUTENSOR( cutensorInitTensorDescriptor( handle,
                                                  &descB,
                                                  nmodeB,
                                                  extentB.data(),
                                                  NULL,/*stride*/
                                                  typeB, CUTENSOR_OP_IDENTITY ) );

    cutensorTensorDescriptor_t descC;
    CHECK_CUTENSOR( cutensorInitTensorDescriptor( handle,
                                                  &descC,
                                                  nmodeC,
                                                  extentC.data(),
                                                  NULL,/*stride*/
                                                  typeC, CUTENSOR_OP_IDENTITY ) );

    printf("Initialize cuTENSOR and tensor descriptors\n");


    // 计算各个张量的alignmentRequirement
    uint32_t alignmentRequirementA;
    CHECK_CUTENSOR( cutensorGetAlignmentRequirement( handle,
                                                     A_d,
                                                     &descA,
                                                     &alignmentRequirementA) );

    printf("alignmentRequirementA: %d\n", alignmentRequirementA);

    uint32_t alignmentRequirementB;
    CHECK_CUTENSOR( cutensorGetAlignmentRequirement( handle,
                                                     B_d,
                                                     &descB,
                                                     &alignmentRequirementB) );
    printf("alignmentRequirementB: %d\n", alignmentRequirementB);

    uint32_t alignmentRequirementC;
    CHECK_CUTENSOR( cutensorGetAlignmentRequirement( handle,
                                                     C_d,
                                                     &descC,
                                                     &alignmentRequirementC) );

    printf("alignmentRequirementC: %d\n", alignmentRequirementC);

    printf("Query best alignment requirement for our pointers\n");

    // 定义并初始化 cutensorContractionDescriptor_t
    cutensorContractionDescriptor_t desc;
    CHECK_CUTENSOR( cutensorInitContractionDescriptor( handle,
                                                       &desc,
                                                       &descA, modeA.data(), alignmentRequirementA,
                                                       &descB, modeB.data(), alignmentRequirementB,
                                                       &descC, modeC.data(), alignmentRequirementC,
                                                       &descC, modeC.data(), alignmentRequirementC, // 覆盖C作为输出
                                                       typeCompute));

    // 定义find
    cutensorContractionFind_t find;
    CHECK_CUTENSOR( cutensorInitContractionFind(
            handle, &find,
            CUTENSOR_ALGO_DEFAULT) );

    printf("Initialize settings to find algorithm\n");


    // 根据desc和find计算空间大小
    size_t worksize = 0;
    CHECK_CUTENSOR( cutensorContractionGetWorkspaceSize(handle,
                                                        &desc,
                                                        &find,
                                                        CUTENSOR_WORKSPACE_RECOMMENDED,
                                                        &worksize ) );

    // 定义空间
    void *work = nullptr;
    if(worksize > 0)
    {
        if( cudaSuccess != cudaMalloc(&work, worksize) ) // This is optional!
        {
            work = nullptr;
            worksize = 0;
        }
    }

    printf("Query recommended workspace size and allocate it\n");


    // 定义plan
    cutensorContractionPlan_t plan;
    CHECK_CUTENSOR( cutensorInitContractionPlan(handle,
                                                &plan,
                                                &desc,
                                                &find,
                                                worksize) );

    printf("Create plan for contraction\n");


    // 执行运算
    typedef float floatTypeCompute;
    floatTypeCompute alpha = (floatTypeCompute)1.1f;
    floatTypeCompute beta  = (floatTypeCompute)0.9f;

    CHECK_CUTENSOR(cutensorContraction(handle,
                              &plan,
                              (void*)&alpha, A_d,
                              B_d,
                              (void*)&beta,  C_d,
                              C_d,
                              work, worksize, 0
                              ));
    cudaDeviceSynchronize();
    printf("Execute contraction from plan\n");

    // 释放资源
    CHECK_CUTENSOR( cutensorDestroy(handle) );

    if ( A_d ) cudaFree( A_d );
    if ( B_d ) cudaFree( B_d );
    if ( C_d ) cudaFree( C_d );
    if ( work ) cudaFree( work );

    printf("Successful completion\n");
}

cuTensor元素级别运算的API相对简单一些，无需像cutensorContraction()中一样的前置操作。元素级别运算函数一共有3个，前面表4.5中有过简单介绍，下面我们以其中支持2个向量输入的函数 cutensorElementwiseBinary()为例，介绍张量元素级别运算。

cutensorStatus_t cutensorElementwiseBinary(
    const cutensorHandle_t *handle, 
    const void *alpha, 
    const void *A, 
    const cutensorTensorDescriptor_t *descA, 
    const int32_t modeA[], 
    const void *gamma, 
    const void *C, 
    const cutensorTensorDescriptor_t *descC, 
    const int32_t modeC[], 
    void *D, 
    const cutensorTensorDescriptor_t *descD, 
    const int32_t modeD[], 
    cutensorOperator_t opAC, 
    cudaDataType_t typeScalar, 
    cudaStream_t stream
)

cutensorElementwiseBinary()实现了如下运算：

输入张量A、C先进行转置运算(Π)，然后再进行独立的一元运算（ψ），最后进行A、C间元素的二元运算（Φ），计算结果放入到D中后再进行一次转置运算，最后输出。其中转置运算和一元运算在A、C、D各自的descriptor中进行定义。cutensorElementwiseBinary()具体的使用样例如下：

void cuTensorElementProcess(){
    typedef float floatTypeA;
    typedef float floatTypeC;

    // 定义cuda数据类型
    cudaDataType_t typeA = CUDA_R_32F;
    cudaDataType_t typeC = CUDA_R_32F;

    printf("Include headers and define data types\n");


    // 定义mode和extent, A和C的维度要相同
    std::vector<int> modeA{'m','h','k','n'};
    std::vector<int> modeC{'m','h','k','n'};

    int nmodeA = modeA.size();
    int nmodeC = modeC.size();

    std::unordered_map<int, int64_t> extent;
    extent['m'] = 96;
    extent['n'] = 96;
    extent['u'] = 96;
    extent['v'] = 64;
    extent['h'] = 64;
    extent['k'] = 64;

    std::vector<int64_t> extentA;
    for(auto mode : modeA)
        extentA.push_back(extent[mode]);

    std::vector<int64_t> extentC;
    for(auto mode : modeC)
        extentC.push_back(extent[mode]);

    printf("Define modes and extents\n");


    // 计算张量A、C元素数量和存储所需字节数量
    size_t elementsA = 1;
    for(auto mode : modeA)
        elementsA *= extent[mode];

    size_t elementsC = 1;
    for(auto mode : modeC)
        elementsC *= extent[mode];

    size_t sizeA = sizeof(floatTypeA) * elementsA;
    size_t sizeC = sizeof(floatTypeC) * elementsC;

    // 分配统一内存
    float *A_d, *C_d;
    cudaMallocManaged((void**)&A_d, sizeA);
    cudaMallocManaged((void**)&C_d, sizeC);

    // 初始化数据
    for(int64_t i = 0; i < elementsA; i++)
        A_d[i] = (((float) rand())/RAND_MAX - 0.5)*100;
    for(int64_t i = 0; i < elementsC; i++)
        C_d[i] = (((float) rand())/RAND_MAX - 0.5)*100;

    printf("Allocate, initialize and transfer tensors\n");


    cutensorHandle_t* handle;
    CHECK_CUTENSOR( cutensorCreate(&handle));

    // 定义并初始化 cutensorTensorDescriptor
    cutensorTensorDescriptor_t descA;
    CHECK_CUTENSOR( cutensorInitTensorDescriptor( handle,
                                                  &descA,
                                                  nmodeA,
                                                  extentA.data(),
                                                  NULL,/*stride*/
                                                  typeA, CUTENSOR_OP_IDENTITY ) );


    cutensorTensorDescriptor_t descC;
    CHECK_CUTENSOR( cutensorInitTensorDescriptor( handle,
                                                  &descC,
                                                  nmodeC,
                                                  extentC.data(),
                                                  NULL,/*stride*/
                                                  typeC, CUTENSOR_OP_IDENTITY ) );

    printf("Initialize cuTENSOR and tensor descriptors\n");


    // 执行运算
    typedef float floatTypeCompute;
    floatTypeCompute alpha = (floatTypeCompute)1.1f;
    floatTypeCompute gamma  = (floatTypeCompute)0.9f;
    cudaDataType_t typeScalar = CUDA_R_32F;


    CHECK_CUTENSOR(cutensorElementwiseBinary(handle,
                                             (void*) &alpha,A_d, &descA, modeA.data(),
                                             (void*) &gamma,C_d, &descC, modeC.data(),
                                             C_d, &descC, modeC.data(),
                                             CUTENSOR_OP_ADD, typeScalar, 0));
    cudaDeviceSynchronize();
    printf("Execute elementwiseBinary\n");

    // 释放资源
    CHECK_CUTENSOR( cutensorDestroy(handle) );

    if ( A_d ) cudaFree( A_d );
    if ( C_d ) cudaFree( C_d );

    printf("Successful completion\n");
}

cutensor中通用规约运算规则如下：

API比较简单，在形式上与元素级别API有些像，也不需要前置操作，使用cutensorReduction()函数可以实现：

cutensorStatus_t cutensorReduction(
    const cutensorHandle_t *handle,
    const void *alpha, 
    const void *A,
    const cutensorTensorDescriptor_t *descA,
    const int32_t modeA[],
    const void *beta,
    const void *C,
    const cutensorTensorDescriptor_t *descC,
    const int32_t modeC[],
    void *D,
    const cutensorTensorDescriptor_t *descD,
    const int32_t modeD[],
    cutensorOperator_t opReduce,
    cutensorComputeType_t typeCompute,
    void *workspace,
    uint64_t workspaceSize,
    cudaStream_t stream
)

cutensorReduction()使用前需要定义workspace，使用cutensorReductionGetWorkspaceSize可以计算workspaceSize：

cutensorStatus_t cutensorReductionGetWorkspaceSize(
    const cutensorHandle_t *handle,
    const void *A,
    const cutensorTensorDescriptor_t *descA,
    const int32_t modeA[],
    const void *C,
    const cutensorTensorDescriptor_t *descC,
    const int32_t modeC[],
    const void *D,
    const cutensorTensorDescriptor_t *descD,
    const int32_t modeD[],
    cutensorOperator_t opReduce,
    cutensorComputeType_t typeCompute,
    uint64_t *workspaceSize
)

得到workspaceSize后申请workspace设备内存便可进行cutensorReduction()运算，具体样例如下：

void cuTensorReductionProcess(){
    typedef float floatTypeA;
    typedef float floatTypeC;

    // 定义cuda数据类型
    cudaDataType_t typeA = CUDA_R_32F;
    cudaDataType_t typeC = CUDA_R_32F;

    printf("Include headers and define data types\n");


    // 定义mode和extent
    std::vector<int> modeA{'m','h','k','n'};
    std::vector<int> modeC{'h','k','n'}; // 只在1个维度规约
    // std::vector<int> modeC{'n'}; // 在3个维度进行规约

    int nmodeA = modeA.size();
    int nmodeC = modeC.size();

    std::unordered_map<int, int64_t> extent;
    extent['m'] = 96;
    extent['n'] = 96;
    extent['u'] = 96;
    extent['v'] = 64;
    extent['h'] = 64;
    extent['k'] = 64;

    std::vector<int64_t> extentA;
    for(auto mode : modeA)
        extentA.push_back(extent[mode]);

    std::vector<int64_t> extentC;
    for(auto mode : modeC)
        extentC.push_back(extent[mode]);

    printf("Define modes and extents\n");


    // 计算张量A、C元素数量和存储所需字节数量
    size_t elementsA = 1;
    for(auto mode : modeA)
        elementsA *= extent[mode];

    size_t elementsC = 1;
    for(auto mode : modeC)
        elementsC *= extent[mode];

    size_t sizeA = sizeof(floatTypeA) * elementsA;
    size_t sizeC = sizeof(floatTypeC) * elementsC;

    // 分配统一内存
    float *A_d, *C_d;
    cudaMallocManaged((void**)&A_d, sizeA);
    cudaMallocManaged((void**)&C_d, sizeC);

    // 初始化数据
    for(int64_t i = 0; i < elementsA; i++)
        A_d[i] = (((float) rand())/RAND_MAX - 0.5)*100;
    for(int64_t i = 0; i < elementsC; i++)
        C_d[i] = (((float) rand())/RAND_MAX - 0.5)*100;

    printf("Allocate, initialize and transfer tensors\n");


    cutensorHandle_t* handle;
    CHECK_CUTENSOR( cutensorCreate(&handle));

    // 定义并初始化 cutensorTensorDescriptor
    cutensorTensorDescriptor_t descA;
    CHECK_CUTENSOR( cutensorInitTensorDescriptor( handle,
                                                  &descA,
                                                  nmodeA,
                                                  extentA.data(),
                                                  NULL,/*stride*/
                                                  typeA, CUTENSOR_OP_IDENTITY ) );


    cutensorTensorDescriptor_t descC;
    CHECK_CUTENSOR( cutensorInitTensorDescriptor( handle,
                                                  &descC,
                                                  nmodeC,
                                                  extentC.data(),
                                                  NULL,/*stride*/
                                                  typeC, CUTENSOR_OP_IDENTITY ) );

    printf("Initialize cuTENSOR and tensor descriptors\n");


    // 执行运算
    typedef float floatTypeCompute;
    floatTypeCompute alpha = (floatTypeCompute)1.1f;
    floatTypeCompute beta  = (floatTypeCompute)0.9f;
    cutensorComputeType_t typeCompute = CUTENSOR_COMPUTE_32F;
    uint64_t workSpaceSize;

    cutensorReductionGetWorkspaceSize(handle,
                                      A_d, &descA, modeA.data(),
                                      C_d, &descC, modeC.data(),
                                      C_d, &descC, modeC.data(),
                                      CUTENSOR_OP_ADD, typeCompute, &workSpaceSize);

    void *workspace = nullptr;
    if(workSpaceSize > 0)
    {
        if( cudaSuccess != cudaMalloc(&workspace, workSpaceSize) ){ // This is optional!

            workspace = nullptr;
            workSpaceSize = 0;
        }
    }
    printf("Query recommended workspace size and allocate it\n");

    CHECK_CUTENSOR(cutensorReduction(handle,
                                             (void*) &alpha,A_d, &descA, modeA.data(),
                                             (void*) &beta,C_d, &descC, modeC.data(),
                                             C_d, &descC, modeC.data(),
                                             CUTENSOR_OP_ADD, typeCompute,workspace,workSpaceSize, 0));



    cudaDeviceSynchronize();
    printf("Execute reduction\n");

    // 释放资源
    CHECK_CUTENSOR( cutensorDestroy(handle) );

    if ( A_d ) cudaFree( A_d );
    if ( C_d ) cudaFree( C_d );

    printf("Successful completion\n");
}

结语

作为一名算法工作者，掌握cuda编程一方面可以加深对常用算法框架底层工程设计的理解；另一方面在某些具体的业务场景下，当通用的框架性能到达瓶颈而又没有合适的第三方库支持的时候，可以尝试自己开发一些cuda算子，对线上推理服务的性能进行优化。本文对cuda基础知识进行了较为全面的梳理，在此基础上未来计划在新的篇章中继续对一些更高级的cuda库进行整理，同时也会对一些经典算法实现中使用到cuda的地方进行解读。