记录学习Cuda编程的项目

任务清单：配置GPU环境以及创建GPU版的Hello world!

主要参考这篇文章

首先检查N卡的驱动，并安装对应的CUDA，然后安装CUDNN

因为我的是GTX1060，所以只能支持CUDA11.4，注意这里我才过坑，安装了11.5,11.6版本在后面对cu程序进行编译的时候会报错，因此要注意显卡能支持的最高CUDA版本，然后根据CUDA版本确定 Visual Studio的版本，

因为官网默认下载最新版本，实际上是不支持古早的CUDA版本的，因此可以根据上面文章提到连接查询官方CUDA支持的VS范围，注意此处的版本十分重要，如果版本错误会导致后面编译失败

我最开始下载的是VS 2022 community 的版本，即使在里面选择安装2019,2017的编译器仍然会报错误，提示头文件找不到，或者版本的问题，改回2019版本的就没问题了。

Linux环境只需要按照官方手册依次安装适合你GPU型号的驱动和CUDA即可，这里不做过多的赘述。

对于GPU编程主要分为两部分，一部分是主机(host)，另一部分则是设备(Device)

主机代码主要在CPU上运行，设备代码则是在GPU上运行完成

我们可以在CPU上运行逻辑分支，然后调用GPU进行计算，但这里GPU的计算与主机是异步的。

调用GPU的函数之后，控制权会返回给CPU，主机代码依然会往下运行，同时CPU后面的代码可以主动调用同步操作cudaDeviceSynchronize函数来同步等待GPU完成计算返回结果。

当然也不是每次调用GPU进行计算之后都需要调用同步操作。GPU和CPU之间的同步分为显式同步和隐式同步两种。 这里调用cudaDeviceSynchronize被称为显式同步。隐式同步则是在函数调用之后会自动等待GPU完成操作。

例如:cudaMemcpy在主机和设备间拷贝数据时，这里CPU也会阻塞来隐式同步数据。

GPU编程主要语法与C语言类似，主要是多了一些特殊标识符

$$ \begin{array}{c|c|c|c} \hline 限定符 & 执行 & 调用 & 备注\\ __{global}__ & 设备端执行 & 主机调用，也可以从算力3以上的设备中调用 & 返回类型必须是void \\ __{device}__ & 设备端执行 & 仅能从设备中调用 \\ __{host}__ & 主机端执行 & 仅能从主机调用 & 可省略 \\ \hline \end{array} $$

对于GPU调用时，其线程创建模型为两个部分，第一个部分叫网格(grid),第二个叫线程块(block),线程块包含了N个线程，结构图如下：

                                            每一个Grid
----------------------------      -------------------------------
|        |        |        |      |         |         |         |
| grid1  | grid2  | grid3  |      | block1  | block2  | block3  |
|--------------------------|      |-----------------------------|
|        |        |        |      |         |         |         |
|        |        |        | -->  |         |         |         |
|--------------------------|      |-----------------------------|
|        |        |        |      |         |         |         |
|        |        |        |      |         |         |         |
----------------------------      -------------------------------

Grid和Block都由一个三维坐标标识(dim3)，每个(x,y,z) 唯一确定一个网格或者一个块。

因此我们可以创建一个维度(3,3,3)的网格，每个网格又是一个维度为(4,4,4)的线程块。因此总线程数为 $3^3(网格数) * 4^3(线程数)$

1. 网格的维度由线程块的数量来表示
2. 线程块的维度由线程数来表示

在代码MatrixSum中：

dim3 block = (16);
dim3 grid = ((matrix_len + block.x - 1) / block.x);

这里计算所需要网格数的代码((matrix_len + block.x - 1) / block.x) 是为了让所有的数据能被网格所包括，因此加上block.x-1,这也是为了规避除法的向下取整，导致无法囊括到所有数据。

------------------------------------------------------------
| data1 | data2 |       |       |       |       |       |  |
|       |       |       |       |       |       |       |  |
------------------------------------------------------------
                                                            ^
                                                            |
-----------------------------------------------------------------
| block1| block2|       |       |       |       |       |       |
|       |       |       |       |       |       |       |       |
-----------------------------------------------------------------
    1       2       3       4       5       6       7       8

GPU有多个流式多处理器(SM)，执行时每一个SM分配多个线程块，SM分配的线程块数由其中的资源来决定

CUDA中采用单指令多线程架构(SIMT)，每32个线程为一组称为线程束(warp) 线程束中的所有线程同时执行相同的指令

每一个SM都将分配给他的所有线程块的线程按照线程束进行划分，其中一个线程块被分配到一个SM后会一直存在其中，直到完成线程任务

               -------    -------------
           |——>|warp | -->| 32 thread |  --|     ---------
           |   -------    -------------    |---> | block |
-------    |   -------    -------------    |     ---------
| SM  |  --|——>|warp | -->| 32 thread |  --|
-------    |   -------    -------------
           |   -------    -------------
           |——>|warp | -->| 32 thread |
               -------    -------------

warp才是SM上的执行单位，因此不同线程束之间的进度可能会不一致，也就会导致线程块之间不同的线程以不同的速度前进。

编译指令：

nvcc -O3 按照级别3来优化主机代码
nvcc -G -g 生成Debug信息

在硬件上线程都是一维排开，尽管线程块可能是1,2,3维。同时根据threadIdx.x的连续值来划分线程束，对于多维的线程块，同样可以根据下列公式转换成一维排列。

一维线程 (x)： $$ thread_{i} = threadIdx.x $$

二维线程 (x,y)： $$ thread_{i} = threadIdx.x + threadIdx.y \times blockDim.x $$

三维线程块 (x,y,z) $$ thread_{i} = threadIdx.z \times (blockDim.x \times blockDim.y) + threadIdx.y \times blockDim.x + threadIdx.x $$

可以看到线程坐标按照 x为内部维度，y作为外部维度，z作为最外面的维度进行递增。

注意到：$一个线程块中的线程束的数量 = 向上取整({一个线程块中的线程数量 \over 线程束大小})$

因此线程束不会跨线程块进行分配，如果当一个线程块的线程数不是线程束大小的偶数倍会造成资源浪费。并会影响线程束的调用，从而影响效率

GPU对于带有逻辑判断的分支语句，例如if...else...并没有复杂的分支预测的能力。但对于同一个线程束中所有线程必须在同一周期中执行相同的指令，因此分支语句会导致一个线程束中不同的线程走到不同的分支语句，这就是线程束的分化。

线程束的分化会明显导致性能下降

对于代码：chapter3/branch.cu，关闭编译器的分支预测优化功能进行运行测试。 编译：nvcc -g -G branch.cu -o branch

运行：sudo ./ncu --metrics smsp__sass_average_branch_targets_threads_uniform.pct branch

可以看到在RTX3050上分支效率结果为：

mathKernal2:80%
mathKernal3:100%
mathKernal4:71.43%

线程束级别的分支效率几乎没有分化的分支

每个线程束的资源都保存在SM单元中，进行线程束切换的时候不会产生性能损失

同时每个SM内核只有固定数量的共享内存和寄存器组， 因此每个线程消耗的寄存器越少，SM能同时处理的线程束越多 每个线程块消耗的共享内存越多，SM能同时处理的线程块就越少

指令发出和完成这段时钟周期被称为指令延迟

在指令延迟这段时间中，每个调度器都有一个符合条件的线程束可用于执行，那么就能隐藏这段延迟（流水线）

延迟分为：算术指令延迟/内存指令延迟

隐藏延迟所需要的活跃线程束数量可以由如下公式进行计算：$线程束数量 = 线程束的延迟\times吞吐量$

吞吐量由SM中每个周期的操作数确定，SM中以线程束为执行单位，因此一个周期会同时有线程束大小的操作执行(warpSize)，上述式子换算成操作单位可以表示成：$操作数 = 指令延迟\times吞吐量$

吞吐量的单位为：操作/周期 指令延迟单位为：周期

根据操作数的需要同时很容易反推需要多少个线程，线程束以及每个SM至少多少个线程束。

同时不要忽略SM的线程束数量同时也受限于硬件资源

占用率的计算公式如下： $$ 占用率 = {活跃线程束 \over 最大线程束} $$

最大线程束数量可以在设备信息中查询，通过cudaGetDeviceProperties获取设备参数信息的结构体，maxThreadsPerMultiProcessor/warpSize 可以得到最大的线程束数量

例如本机为NVIDIA GeForce RTX 3050 Laptop GPU,其中最大的线程束数量为maxThreadsPerMultiProcessor/warpSize = 1536/32 = 48

可以对编译器添加参数--ptxas-options=-v 或在CMakeLists.txt中添加set(CUDA_NVCC_FLAGS --ptxas-options=-v) 来查看每个核函数使用了多少个寄存器以及共享内存的资源。

提高利用率我们需要设置合适的线程块配置。过大过小都会影响资源的利用率。

小线程块：会在所有资源被充分利用之前到达SM的线程束数量的限制

大线程块：会导致每个SM中每个线程可用的硬件资源过少

网格和线程块大小划分准则：

1. 保持每个块中的线程是线程束大小的倍数
2. 避免块太小，每个块至少有表128或者256个线程
3. 根据内核资源调整块的大小
4. 块的数量要远多于SM的数量，达到足够的并行

线程块之间的同步: __device void __syncthreads() 系统级同步：cudaDevicesSynchronize()

线程块是以线程束为单位执行，因此同一个线程块中的不同线程束会处于不同的程序点，所以提供__syncthreads来同步块中的所有线程

同步之前线程块中的所有线程产生的共享内存或全局内存的修改操作等，在同步后会对线程块中的所有线程可见，并且是安全的。因此可用于线程间的通讯

老版本cuda可以使用nvprof来查看程序运行的各种指标,现在逐渐替换成ncu工具来检查程序运行的性能指标。

更多参数对照可以参考官方提供的CUDA手册《NSIGHT COMPUTE COMMAND LINE INTERFACE》 Ref:

nvprof --metrics achieved_occupancy ./xxx
ncu --metrics sm__warps_active.avg.pct_of_peak_sustained_active ./xxx

该指标越高越好

核函数运行时间

nvprof ./xxx
nsys nvprof ./xxx

在矩阵求和的例子中，不同参数下SM占用率和耗时分别为：

(32,32):53.24%  22.10ms
(32,16):70.08%  21.60ms
(16,32):73.10%  21.04ms
(16,16):78.04%  21.19ms

从上面结果可见第二种比第一种有更多的块，有着更好性能，但第四中比第三种也有更多的块，但性能并没有更好的提升，虽然SM占用率较高

./ncu --metrics l1tex__t_bytes_pipe_lsu_mem_global_op_ld.sum.per_second ./xxx

可以看到结果如下：

(32,32):78.99 Gb/s 
(32,16):89.97 Gb/s
(16,32):90.65 Gb/s
(16,16):91.21 Gb/s

其中第四中具有最高的内存读取效率，但运行时间并没有第三种情况快，所以高内存读取效率并不代表着有最好的性能

./ncu --metrics smsp__sass_average_data_bytes_per_sector_mem_global_op_ld.pct ./xxx

可以看到结果如下：

(32,32):100%
(32,16):100%
(16,32):100%
(16,16):100%

我们可以结合上述三个指标来判断最佳的块组合

对于一个数组求和或者求最大值的问题而言，使用CPU顺序求和在数据量比较大的情况下是比较复杂的，时间复杂度O(n),因此可以考虑使用并行归约来进行求解。

并行归约的主要**是，对于一个线程块划分数据区域，每个数据块处理一个区域，然后逐层将数据结果汇总

参见代码:reduceInteger.cu:matrixSum

方法matrixSum主要方式是如图所示：

            -------------------------------------
data:       |  0  |  1  |  2  |  3  |  4  |  5  |
            -------------------------------------
               |     |     |     |     |     |
            -------------------------------------
threadId:   |  0  | (1) |  2  | (3) |  4  | (5) |
            -------------------------------------
               |     |     |     |     |     |   
               -------     -------     ------- 
               |           |           |
            -------------------------------------
data:       |  0  |     |  2  |     |  4  |     |
            -------------------------------------

其中偶数线程负责进行相邻步长的相加,奇数线程则不进行处理。matrixSumNeighbored 方法则是在线程的处理上进行了排序整合

方法matrixSumNeighbored主要方式是如图所示：

            -------------------------------------
data:       |  0  |  1  |  2  |  3  |  4  |  5  |
            -------------------------------------
               |     |     |     |     |     |
            -------------------------------------
threadId:   |  0  |     |  1  |     |  2  |     |
            -------------------------------------
               |     |     |     |     |     |   
               -------     -------     ------- 
               |           |           |
            -------------------------------------
data:       |  0  |     |  2  |     |  4  |     |
            -------------------------------------

使用相邻的线程分别间隔一个步长区处理不同的数据区域进行归并

这两种方法总线程数并没变，只是改变了活跃线程的排布，对于512个线程而言，第一种方法按照32划分线程束共有16个线程束. 每个线程束中有一半的线程处于非活跃状态，而第二种将活跃线程全部集中在前8个线程束中,由于线程束中每个线程执行的命令操作要保持一样，因此第一种方法会因为线程的分化影响性能。而第二种则完全避免了这种情况。

所以确保线程束中减少线程分化有利于提高性能

同第二种方法类似，matrixSumInterleaved使用交叉归约的方式进行求和。主要流程如下：

            -------------------------------------
data:       |  0  |  1  |  2  |  3  |  4  |  5  |
            -------------------------------------
               |     |     |     |     |     |
               -------------------     |     |
               |     |     |           |     |
               V     -------------------     |
                     |     |                 |
                     V     -------------------
                           |
                           V
            -------------------------------------
threadId:   |  0  |  1  |  2  |     |     |     |
            -------------------------------------
               |     |     |
            -------------------------------------
data:       |  0  |  1  |  2  |     |     |     |
            -------------------------------------

每个线程先按数据长度的一半进行求和处理，然后结果汇集到前半段，然后继续处理前半段数据进行求和。注意边界判断，每次只有小于步长的线程才进行计算。并且同方法二一样活跃线程集中在块的前半部分。

__syncthreads块内数据同步，每次加法计算完，确保每块其他线程完成计算。

这三种方法的SM利用率，内存加载速率以及全局加载率如下：

matrixSum            96.80%  91.14Gb/s 25.05%
matrixSumNeighbored  93.05% 191.70Gb/s 25.05%
matrixSumInterleaved 91.92%  59.95Gb/s 96.78%

从结果可以看到交错合并计算效率最高。

展开归约是指在一个循环里将一次执行复制成多次的计算，这样可以将所有归约的次数继续缩短。 展开的次数称为展开因子

对上一节方法三进行更改，每个线程块处理两块数据内容。在进行块内合并之前，先进行块间数据累加，然后进行块内求和。可以看到效率更快。相比较与方法三拥有更高的内存加载速率。

所以一个方法中有更多的独立内存操作，会导致内存有更高的吞吐量，从而性能提升

展开归约的吞吐量和展开因子之间成正比。

对于方法四除了在块间进行归约展开以外，还能在块内归约展开。见方法五，将最后步长为[8,4,2,1]的结果展开求和，能使得合并效率更高。

如果已知一个线程内的循环次数，可以手动将其完全展开，以扩大内存的吞吐量，提高执行速率

对于一个cuda函数，我们可以使用模板添加一个变量，这个变量会在编译器编译的时候进行编译优化。使用方式如下：

template<data_type data_name> __global__ void function(){}

编译器在进行分支检查的时候，会自动删除关于data_name别且不满足条件的分支，以达到减少分支的目的。

即在设备函数运行过程中，根据需要再次调用设备函数进行计算。典型的例子就是嵌套。

嵌套执行分为了父子线程，父线程创建并调用子线程，只有在子线程完成后父线程才会完成。

显示同步：在父线程中设立栅栏点，父线程会等待子线程完成到达栅栏点。 隐式同步：父线程结束时会等待子线程完成。

动态并行的最大嵌套深度为24

1. 线程束是SM的执行调度单位，线程束的大小由设备warpSize确定
2. 确保块中的线程数是warpSize的整数倍
3. 块中的线程数不能太大也不能太少，至少[128-256]
4. 块的数量不能太少，以遮盖内存延迟
5. 归约问题
6. 展开归约有助于减少循环，增加可并行操作

1. 寄存器：线程私有
2. 本地内存：线程私有，生命周期在于线程
3. 共享内存：线程块拥有，生命周期同线程块
4. 全局内存：全部线程拥有，生命周期同整个程序
5. 常量内存：全部线程拥有（只读），生命周期同整个程序
6. 纹理内存：全部线程拥有（只读），生命周期同整个程序

寄存器用于保存线程中自定义的变量。超出硬件上限会使用本地内存替代，并且会影响运行效率

硬件会更具SM上线程块的数量对寄存器进行平分，因此核函数需要的寄存器越少，SM上能够活跃的线程块越多

除了依赖硬件的自行资源划分，我们也可以通过自定义块的数量来改变寄存器的数量。

__global__ void __lanch_bounds__(maxThreadsPerBlock,minBlocksPerMultiprocessor)

maxThreadsPerBlock：每个块最大的线程数

minBlocksPerMultiprocessor：每个SM中最小的常驻线程块数

同样可以使用编译参数-maxrregcount=32来指定所有核函数中寄存器数量的最大值

占用空间较大的结构体或者数组会保存在本地内存，以及不满足寄存器保存条件的变量都会保存在本地内存，本地内存本质同全局内存在同一块区域

使用__shared__标识的变量，存放于共享内存中，共享内存被线程块共享，并且在片上，因此速率更高。共享内存的使用同样会影响SM上活跃的线程束

共享内存在核函数内声明和定义

这块区域可以配合__syncthreads()作线程块的同步用

共享内存同一级缓存共同使用64K的片上内存,可以通过cudaFuncSetCacheConfig来指定共享内存和一级缓存的大小。

使用__constant__标识，常量内存必须定义在核函数外以及全局空间内。大小只有固定的64K，常量内存拥有只读属性，因此要从主设备通过cudaMemcpyToSymbol往里面写数据. 这个内存区域的读写能力要优于全局内存，适合保存的例如数学系数等

纹理内存拥有硬件支持的滤波和插值

全局内存可以被所有的SM访问到，全局内存的创建分为静态创建和动态创建

静态创建：使用__device__标识的变量

动态创建：cudaMalloc()来创建全局内存

全局内存的访问通过32,64,128字节的内存事务进行访问，并且会内存对齐:这三种事务加载时首地址必须被32,64,128整除。

一般来说用于加载内存请求的事务越多，数据的使用效率越低。

使用__device__标识的变量从主机赋值时，通过cudaMemcpyToSymbol或者cudaMemcpyFromSymbol来写入/读取数据

同样可以通过cudaGetSymbolAddress获取全局变量的地址,然后配合cudaMemcpy来复制数据到全局内存变量中。

cuda程序内存分配流程如下：

--------            -----------------------------------
|主机内存| --copy--> | cuda驱动在主机内存中创建的固定内存区域|
--------            -----------------------------------
                                    |
                                  copy
                                    |
                                    V
                                --------
                                |设备内存|
                                --------

主机进行内存拷贝的同时多了一步创建固定内存区域的步骤，我们可以通过cudaMallocHost直接在主机中创建固定内存区域，供设备访问。并通过cudaFreeHost来释放，这样效率会比[cudaMalloc,cudaMemcpy]一套要高。同时这个区域主机也能直接访问。

注：固定内存分配和释放的成本更高，但是提供了更高的传输吞吐量，当数据量较大的使用固定内存分配效果会更好。

通常主机和设备间不能直接互相访问变量和数据，零拷贝内存提供了一个特例。

使用零拷贝主要优势有：

1. 设备内存不足可以利用主机内存
2. 避免主机和设备间的显示内存传输

零拷贝也是属于固定内存，通过cudaHostAlloc创建一块主机内存并映射到设备内存区域，同时使用cudaFreeHost来释放.

在设备端通过cudaHostGetDevicePointer(devicePoint,hostPoint)来获取这块映射的内存指针。

devicePoint是指向hostPoint在设备上的映射，可以被设备使用访问数据。

注：设备端每次映射的传输与访问都会经过PCIe总线，会导致效率过低，这也是我们需要避免过多的主从设备见数据传输的原因。

因此共享少量数据的时候使用零拷贝内存是个不错的选择，数据量越大传输延迟越高。

对于计算能力2.0以上的设备支持统一虚拟寻址(UVA)，设备不需要通过cudaHostGetDevicePointer来手动获取主机上的内存映射，直接使用cudaHostAlloc分配的主机指针即可，可以交由设备自行进行寻址。

在统一内存寻址中，提供了一块内存池，底层系统在统一内存空间自动在主机和设备之间进行数据传输。

通过__managed__在文件范围和全局范围内添加一个静态托管变量，或者通过主机端调用cudaMallocManaged创建一块动态分配的托管内存

GPU内存请求特性：核函数的内存请求，通常是在物理内存上以128字节或者32字节内存事务来实现。

一级缓存一行是128字节，映射到内存中就是128字节的对齐段。

二级缓存则是32字节的内存事务进行访问

可以通过编译器选择是否开启一级缓存

禁止一级缓存编译参数：-Xptxas -dlcm=cg

开启一级缓存编译参数：-Xptxas -dlcm=ca

核函数对数据的访问产生一个内存事务，如果内存事务的首地址是对一级缓存大小或者二级缓存大小的偶数倍时，就出现了访问对齐，同一个线程束中的线程访问的数据是一个连续的内存块就会出现合并内存访问。

合并和对齐会最大化全局内存吞吐量

全局内存的加载顺序，如果一级缓存没被禁用，那么全局内存首先由一级缓存(128字节)尝试加载，如果被禁用则用二级缓存(32字节)尝试加载，二级缓存也被禁用那么会直接交由DRAM(128字节)的内存事务实现。

注意到，一级缓存的加载是128字节为单位，对于一个线程束而言，每个线程请求一个4字节变量，并且分布均匀不重复的内存空间，一次加载事务刚好满足，此时利用率为100%,但如果需要的数据不关于128字节对齐，那么需要两个加载事务，利用率为${128 \over 128 \times 2} = 50%$，这个时候使用粒度更细的32字节加载,可以看到加载率为${128 \over 32 \times 5} = 80%$，所以开启一级缓存在有的情况下并一定能提高总线加载效率

二级缓存的粒度是同内存粒度的

数据写入内存前都会通过二级缓存进行操作，也就是以32字节的粒度进行写入。 内存事务可以分为一段，两段，四段事务执行。分别对应32,64,128字节

两种数据布局分别如下所示：

AoS:                    SoA:           
---------------------   ---------------------
|x|y|x|y|x|y|...    |   |x|x|x|y|y|y|...    |
---------------------   --------------------- 

AoS的每次加载都会加载一半的x一半的y值,对于只用某一项的核函数会损失一般的加载效率。

对于单一命令多数据的架构，使用SoA更合适。例如：线程束中每个线程处理的指令相同，所以使用的数据项也是相同的，利于提高数据利用率的方式就是SoA模式。

1. 对齐和合并数据加载
2. 足够多的线程块，掩盖内存延迟
3. 展开技术

其中展开技术提高IO并行操作带来的效率提高最为显著，并且同时能减少内存读写事务。

矩阵转置分为两种:按行读取转置成列，按列读取转置成行。

    --------col----------         --------row----------
    |             ny    |         |     |             |
    |              |    |         |     |             |
    ------nx-------N-----         |     nx            |
    |              |    |         |     |             |
 row|              |    |      col|     |             |
    |              |    |         |     |             |
    |              |    |         --ny--N--------------
    ---------------------         |     |             |
                                  ---------------------

在这里nx作为横轴(行方向)的变化量，对于N(nx,ny) 在转置矩阵中的关系为，N(nx,ny) --> N(ny,nx)转换成一维矩阵的表达式：$ny \times col + nx \rightarrowtail nx \times row + ny$

    --------col----------         --------row----------
    |             nx    |         |     |             |
    |              |    |         |     |             |
    ------ny-------N-----         |     ny            |
    |              |    |         |     |             |
 row|              |    |      col|     |             |
    |              |    |         |     |             |
    |              |    |         --nx--N--------------
    ---------------------         |     |             |
                                  ---------------------

在这里nx作为竖轴(列方向)的变化量，对于N(ny,nx) 在转置矩阵中的关系为，N(ny,nx) --> N(nx,ny)转换成一维矩阵的表达式：$nx \times col + ny \rightarrowtail ny \times row + nx$

单独测试数据读取和复制(copyRow,copyCol)有如下结果：

        读取效率, 写入效率, 读取带宽,  写入带宽
copyCol: 12.5%   12.5%   74.23Gb/s  74.28Gb/s
copyRow:  100%    100%   41.68Gb/s  41.74Gb/s

对于按列转置和按行转置(transposeNaiveRow,transposeNaiveCol)可以发现有如下结果：

    读取效率, 写入效率,     读取带宽,               写入带宽
Col: 12.5%   100%   194.82Gb/s(交叉读取)  24.35Gb/s(合并写入)
Row:  100%  12.5%     9.04Gb/s(合并读取)  72.32Gb/s(交叉写入)

可以看到按列转置的读取带宽最高，这是因为在一级缓存的作用下，每个线程加载128字节的数据，但是只使用其中一个，但是其相邻线程进行操作时会有更高的缓存命中率。

同时使用按行拷贝会有更短的执行时间。按行拷贝时，读取效率和写入效率都会更高，这是因为读取和写入都通过一个事务即可。

对于按列转置而言，相比较与按行具有很高读取带宽相差20倍左右，而写入带宽相差仅仅3倍。这导致按列转置具有更低的执行时间。

当使用展开优化后，展开因子为4,(transposeNaiveRow_unroll4,transposeNaiveCol_unroll4)可以发现有如下结果：

            读取效率, 写入效率,     读取带宽,               写入带宽
Col_unroll4: 12.5%   100%      281.99Gb/s(交叉读取)  35.25Gb/s(合并写入)
Row_unroll4: 100%    12.5%     10.28Gb/s(合并读取)    82.23Gb/s(交叉写入)

按行展开在内存带宽的吞吐量上有些微提升，这也transposeNaiveRow_unroll4执行时间要比transposeNaiveRow少的原因。因为有更高的可并发内存操作。

对于按列展开同样相比较与transposeNaiveCol拥有更高的读取带宽和写入带宽，同样是因为可并行的内存读写操作导致。因此拥有更快的执行效率

对于矩阵操作优先使用按列操作

另一种优化方式可以对网格的坐标系进行优化,传统我们一般使用直角坐标系如下图所示

        直角坐标系(y,x)                     对角坐标系(y,x)
-------------------------        -------------------------
|(0,0)|(0,1)|(0,2)|(0,3)|        |(0,0)|(0,1)|(0,2)|(0,3)|
-------------------------        -------------------------
|(1,0)|(1,1)|(1,2)|(1,3)|        |(1,3)|(1,0)|(1,1)|(1,2)|
-------------------------        -------------------------
|(2,0)|(2,1)|(2,2)|(2,3)|        |(2,2)|(2,3)|(2,0)|(2,1)|
-------------------------        -------------------------
|(3,0)|(3,1)|(3,2)|(3,3)|        |(3,1)|(3,2)|(3,3)|(3,0)|
-------------------------        ------------------------- 

直角坐标系和对角坐标系的转换公式：

$$ y_对 = y_直 \\ x_对 = (x_直 - y_直 + gridDim.x) % gridDim.x $$

对角坐标系转换成直角坐标系 $$ y_直 = y_对 \ x_直 = (x_对 + y_对) % gridDim.x $$

内存的加载最终由DRAM硬件完成，在直角坐标系中，将线程快映射到内存区域时，由于排列过于紧密，很容易出现分区冲突(计算机组成原理存储部分)，而使用对角坐标系能一定程度的拉开相邻线程快加载数据的距离。从而减少分区冲突现象。

但是对角坐标系带来的提升不如按列读取

使用更瘦的块能提高存储的效率

在transpose函数中 将块的大小改为<<<8,32>>>,加载效率和写入效率相比较与<<<16,32>>>提升了一倍

使用统一内存管理同样能提升处理效率

共享内存是片上内存，全局内存是版载内存，共享内存提供了更高的访问速率并有更低的延迟。

在核函数内或者文件范围内，通过__shared__来标识变量。表示该变量是共享内存。或者申明一个共享内存的多维数组。

如果对于数组长度未知，我们可以通过extern来申明一个未知长度的一维数组，这个申明可以在全局或者核函数内，可见范围同其定义位置。

然后通过在每个核函数调用时通过第三个参数来定义大小。这个只支持一维数组。

extern __shared__ int temp[];

function<<<grid,block,byteSize>>>();

同时如果作为第三个参数的方式定义全局内存，那么共享内存也只能作数据存储用，无法导入外部数据。

每个线程块创建时都会创建一定的共享内存，线程块中的所有线程共享该块内存区域，并且具有相同的生命周期，如果多个线程同时访问共享内存的统一变量，该变量会通过多播发送给其他线程。

共享内存等分为32个同样大小的内存模型，叫做存储体，这些存储体可以同时被访问，这个数量同warpSize,如果线程束中的每个线程访问不同的存储体，那么可以用一个内存事务完成。否则就是多个事务，这样会降低加载效率

动态共享内存只能被申明成一维数据

共享内存的定义可以在核函数内也可以在核函数外，在核函数内则作用于仅限于核函数，在核函数外则对所有核函数可见。

多个地址请求落在同一块内存区域时，就会发生存储体冲突。

存储体请求分为三种情况：

1. 并行访问：多个地址请求多个存储体
2. 串行访问：多个地址请求一个存储体
3. 广播访问：单个地址请求一个存储体

对于共享内存的访问模式来说，不同的计算平台有不同的存储体宽度(32位,64位) 也就是每个存储体的一个读写操作的单位。并且连续的32位字映射到连续的存储体中

因此32位的访问平台有如下的地址到存储体的计算映射公式：

$$idx = {address \over 4} % 32$$

访问地址除以4个字节，转换成4字节索引，然后对32个存储体取余

      ------------------------------------   ------
   0x | 00 | 04 | 08 | 0C | 10 | 14 | 18 |...| 80 |
      ------------------------------------   ------
         |    |    |    |    |    |    |       
      ------------------------------------   ------
   /4 | 00 | 01 | 02 | 03 | 04 | 05 | 06 |...| 32 |
      ------------------------------------   ------
         |    |    |    |    |    |    |       
      ------------------------------------   ------
  %32 | 00 | 01 | 02 | 03 | 04 | 05 | 06 |...| 00 |
      ------------------------------------   ------

邻近的字被分到不同的存储体中，所以同一线程束中的线程访问邻近的地址时，不会有存储体冲突。

对于64位的访问平台，地址到存储体的映射公式为：

$$idx = {address \over 8} % 32$$

访问地址除以8个字节，转换成8字节索引，然后对32个存储体取余

注意两个线程访问64中的子字时(每个线程只需要其中的32位数据)，并不会产生存储体冲突，因为只需要一个64字的加载事务。

内存填充是解决存储体冲突的方法之一，通过在列后面添加一列元素可以使得原来数据排序错乱开。在访问同一列元素时可以落到不同的存储体中。主要能够解决按列写存储体冲突的问题。

-----------------          -----------------
| 0 | 1 | 2 | 3 |          | 0 | 1 | 2 | 3 |
-----------------    -->   -----------------
| 4 | 5 | 6 | 7 |          | x | 4 | 5 | 6 |
-----------------          -----------------

注：填充会增加每个线程块使用共享内存的大小，从而影响每个SM中可活跃的线程束数量。

共享内存通过4字节还是8字节访问可通过:cudaDeviceSetSharedMemConfig来进行配置

在前面我们介绍了__syncthreads()用来使块内线程全部同步到指定点

这里介绍一种全局内存同步的方法：内存栅栏，这个方法可以确保栅栏前内存的所有操作，在栅栏后对所有线程是可见的。

内存栅栏分为三个粒度：

1. 线程块级别：__threadfence_block()
2. 网格级别：__threadfence();
3. 系统级别:__threadfence_system()

线程块级别的可以保证线程块对共享内存和全局内存所有的内存操作 在该线程块内是可见的。

网格级别的，会挂起调用的线程，直到在网格级别内保证全局内存的读写在栅栏后是可见的

系统级别则是可以决定主机和设备之间的数据同步，主要针对全局内存，主机的固定内页内存，设备中的内存。

该修饰符可以避免编译器优化：将数据暂存于寄存器或者本地内存中。这个修饰符修饰的变量会实时更新到全局内存中。

更具5.1.3.1中 共享内存的访问模式我们可以看到，访问相邻数据时，是不会产生存储体冲突，因此对于一个二维数组 __shared__ int temp[w][h];而言，其是按照行优先进行存储成一维数组，因此我们使用行优先访问会避免存储体冲突。使用列会产生一半的冲突

更具线程的threadIdx.x连续，推荐使用 temp[threadIdx.y][threadIdx.x]

注：按列写入会产生许多存储体冲突，可以通过添加列数来解决

见代码shared_matrix_transpose.cu,这里实现了通过共享内存和不使用共享内存来进行矩阵转置。并且我们通过全局内存读取事务来衡量效率指标

命令：

sudo ./ncu --metrics l1tex__average_t_sectors_per_request_pipe_lsu_mem_global_op_ld.ratio,l1tex__average_t_sectors_per_request_pipe_lsu_mem_global_op_st.ratio ./transpose

其中转置方法分别可以得到如下的结果：

         读取事务，写入事务
copyRow : 4     ,4
naiveGmem：4     ,32
naiveSmem：4     ,32
naiveSmemCol:4   ,4

可以看到单纯的复制拷贝函数（按行读取，按行写入），读取平均使用4个事务，写入同样也使用4个事务。因为按行读取并不会发生存储体冲突，同理按行也是。

但对于传统转置方法：按行读取,按列写入，无法避免写入时的存储体冲突问题。使用共享内存也是如此。

只有对于使用共享内存暂存元素，同时改变读取方式，读取一列并按行写入数出矩阵则可以避免这个问题。

使用带填充的共享矩阵能进一步降低执行时间。

常量内存同全局内存一样位于片外，对于设备而言只读，对于主机而言可读可写

常量内存中如果所有线程访问同一个内存元素，则是访问最优的。常量内存的访问效率同需要访问的元素个数成反比。

常量内存通过 __constant__来标识，常量内存可以跨文件进行访问，同时只允许主机端通过cudaMemcpyToSymbol进行写入

并且常量内存只能通过静态申请内存空间(数组形式)，无法通过动态申请，同时只能定义在全局文件，如果启用独立编译则可以跨文件可见

更具常量内存的访问特性，多次访问是按照串行进行访问，因此对于系数而言是比较好的存储方式，一个线程束每次读取相同的系数值在这里能达到最优方式。

#pragma unroll表示由cuda编译器将循环自动展开

除了通过一级缓存进行数据读取还能通过只读缓存来读取数据，并且是一个单独的读取通道，这个可以为达到带宽瓶颈的核函数提供另一种读取方式

只读缓存比较与一级缓存有更高的离散性读取能力，其粒度为32字节，而一级缓存为128字节。

可以在读取全局内存的时候通过以下方法指定使用只读内存进行读取：

1. 内部函数__ldg
2. 全局内存添加限定指针，const __restrict__

该访问模式更适合于常量内存

线程束(warp)是SM的执行单位，同时每32个线程划分成一个线程束。线程束洗牌指令可以让线程之间相互获取对方的寄存器中的数据(变量)

洗牌指令分为四种：

1. 广播指令，将一个线程中某个变量值广播到其他线程

通过int __shfl(int var,int srcLane,int width=warpSize)函数进行广播。 var:返回值 srcLane:目标线程在线程束中的ID width：段宽，线程束每隔width进行一次分段

线程束中线程ID根据如下公式计算： $$ laneId = threadIdx.x % 32 \ warpId = threadIdx.x \div 32 $$

每32个线程重新按照[0-31]进行编号,这里32是默认线程束大小warpSize

2. 向右平移，线程束中的线程获取往右平移指定数量的变量值

int __shfl_up(int var,unsigned int delta,int wirth=warpSize) var:返回值 delta:偏移量 width：段宽，线程束每隔width进行一次分段

其中左侧未被覆盖的线程值保持不变

3. 向左平移，线程束中的线程获取往左平移指定数量的变量值

int __shfl_down(int var,unsigned int delta,int wirth=warpSize) var:返回值 delta:偏移量 width：段宽，线程束每隔width进行一次分段

其中右侧未被覆盖的线程值保持不变

4. 异位获取，每个线程更具自身线程ID按位异或来获取变量值

int __shfl_xor(int var,int laneMask,int wirth=warpSize) var:返回值 laneMask:按位异或的值 width：段宽，线程束每隔width进行一次分段

下面介绍几种特殊用法

利用广播来实现循环移动，对每个线程设置不同的参照函数值即可

__shfl(value,threadIdx.x + offset,DIMX);

每个线程分别设置当前线程偏移offset为参考值线程。这样就可以实现循环移动，并且offeset可以为负数

蝴蝶移动利用异或来进行值交换。

int __shfl_xor(value,1,DIMX)

表示使用当前线程的在线程束中的ID与1进行异或，最后得到的值即为需要进行值传递的线程ID。

int __shfl_xor(value,2,DIMX)

表示与2进行异或，这个效果同2次展开操作能得到相同结果，最终是以2个线程为单位进行相邻的交换。但是本质上并不一样这个是因为与2的偶数作为异或因子与2次展开刚好相同，当3作为异或因子时得到的结果则与3次展开不相同。

2次展开蝴蝶交换是在以1为异或参数交换的基础上进行展开。

   threadID  0       1       2 
           -----   -----   -----
           |   |   |   |   |   |   
         -----------------------------
         | 0 | 1 | 2 | 3 | 4 | 5 | 6 |
         -----------------------------
           |       |       |
         -----   -----   -----
第一次交换 | 0 |   | 2 |   | 4 |
         -----   -----   -----
               |       |       |         
             -----   -----   -----
第一次交换     | 1 |   | 3 |   | 5 |
             -----   -----   -----         

cuda核函数之间的调用大部分是异步的，为了缩短程序的整体运行时间，因此提出了流式的概念，用来隐藏核函数之间的延迟。

同一个流中的操作有严格的执行顺序，而不同的流之间则可以相互叠加。(流水线技术)

例如：我们可以利用这一点掩盖一部分的数据传输延迟，将数据传输划分成几个子部分，将子部分在不同的流中实现，通过流的并行处理，可以掩盖部分的传输延迟

流分为两种：隐式声明的流(空流)，显式声明的流(非空流)

空流是默认流，在主函数中没有特别说明和创建非空流，都使用的是默认流

非空流是实现重叠不同操作的关键

对于完整的一套cuda程序应该包含如下步骤：

1. 复制主机数据到设备上
2. 执行核函数
3. 拷贝设备数据到主机上

而对于这三步都可以通过异步进行处理。

cudaMemcpyAsync调用后控制权会立即返回给主设备。其中第五个参数是指定流的对象

cudaStreamCreate用来创建一个内置对象cudaStream_t,这个对象将会是后续异步操作版本的额外参数

$\star$ 注：使用异步传输数据时，必须使用固定内存(cudaMallocHost)或者零拷贝内存(cudaHostAlloc)来分配主机内存。如果不使用操作系统会从物理层面异步移动数据，程序会导致未定义的行为。

创建流后，可以通过调用cudaStreamDestroy来销毁流中的资源，这个函数会自动清空流中相关的资源

1. cudaStreamSynchronize用来阻塞主机，直到指定流完成所有操作。
2. cudaStreamQuery 非阻塞的查询流是否已经完成，如果完成则会返回cudaSuccess

见代码示例:stream.cu,使用nsys进行性能分析后得到如下结果

可以看到在流14拷贝数据的时候，流13正在执行核函数，并且流14和流16的拷贝函数在同时进行只不过流14是DTH,而16是HTD

可见通过流处理能很好的隐藏内存传输延迟

注：

1. 流的数量并不能无限制的增加，受硬件的影响，有些并不能支持很多路并发
2. 逻辑中的函数虽然是异步执行，但仍然按照调用的顺序进行

for(int i=0;i<STREAM;i++){
   cudaMemcpyAsync();
   kernal<<<grid,block,0,streams[i]>>>();
   cudaMemcpyAsync();
}

3. 流同步函数应该在函数逻辑外面，如果放在主逻辑中，每个流注册完函数后，会阻塞在同步函数处，直到该流完成任务，重新开始下一个循环。最终导致仍然是串行完成。

注意核函数的调用方法，第四个参数用来标识该核函数在哪个流中运行，默认是空流。第三个参数在共享内存中已经介绍，用来标识动态分配的共享内存大小。

能并行执行的流的数量并不是无限制增加，流映射到物理硬件上能同时并行的流个数是有限的。 因此就存在一个多路复用的技术，超过硬件所能并行的数量的流会和其他流公用同一个流

因此就有一个虚假的依赖关系，Hyper-Q技术则是用来解决这种情况：创建硬件工作队列，每个流分配到一个对列里面。

可以通过函数cudaStreamCreateWithPriority来分配流的优先级，这个优先级会返回给cudaStream_t句柄中,优先级的查询可以使用cudaDeviceGetStreamPriorityRange 函数

cuda提供了一个插入事件以及查询事件完成的功能，只有当事件前面的所有操作完成则事件才会完成。因此我们可以在流中插入事件点，并检测事件的完成情况来判断流程的进程

通过cudaEventCreate来创建cudaEvent_t事件对象，并通过cudaEventDestroy销毁事件及其相关资源

函数cudaEventRecord(cudaEvent_t,cudaStream_t)用来在cuda流中打上事件点， 可以通过检查这个事件的完成情况，判断流程是否到达指定点。

cudaEventSynchronize函数来阻塞流线程的调用并等待所有流到达该事件点

cudaEventQuery非阻塞查询事件是否完成，完成返回cudaSuccess

cudaEventElapsedTime获取两个事件之间的耗时

使用方法见stream.cu

非空流可以分两种:阻塞流和非阻塞流

非空流中的操作可以被空流(主机流)中的操作阻塞，如果一个非空流是阻塞流，则空流可以阻塞非空流，如果是非阻塞流，则不会阻塞空流中的操作

我们使用cudaStreamCreate创建的流是阻塞流，也就是说主机流可以阻塞该流中的操作。

主机流是默认流也就是空流，同时也是同步流，当操作发布到空流中，该操作执行前cuda会等待所有先前的操作完成，同时非空流阻塞流中的操作会被挂起，直到空流中的操作完成。见下图，展示了空流和非空流的阻塞流的调度过程

可见三个内核是按照时序执行，但是从主机角度来看main函数的打印输出在三个核函数执行前，这三个核函数的调用仍然是异步的。总的来说就是 调用的时候是异步的，但是在核函数执行的过程中仍然是同步的

在多个核函数调度同时可以看到，非空流的阻塞流依然可以并行执行，但在执行过程中被空流断开

函数cudaStreamCreateWithFlags(cudaStream_t,flags)用于创建一个指定行为(阻塞，非阻塞)的流.其中flag支持两个参数：cudaStreamDefault(阻塞),cudaStreamNonBlocking(非阻塞)，非空流的非阻塞流能够和主机同步执行，就解决了上述的问题

因此为了避免主机流对非空流的阻塞行为，有两种解决方法：1. 非空流使用非阻塞流，2. 避免在主机流里调用核函数

例如 cudaMemcpy

设备级别：cudaDeviceSynchronize，阻塞主机线程，等待设备完成所有任务

流级别：cudaStreamSynchronize，阻塞主机线程，等待流的所有操作完成，同样可以使用非阻塞同步函数cudaStreamQuery

事件级别：cudaEventSynchronize，等待核函数到达事件点。同样可以使用非阻塞同步函数cudaEventQuery

跨流同步：cudaStreamWaitEvent，等待某一个流的某一个事件发生，该方法的调用可以和等待的流不在同一个流内，因此可以做到跨流同步。这会保证等待的那个流达到要求的程序点

cudaEventCreateWithFlags可以定义不同类型的事件

cudaEventDefault:默认

cudaEventBlockingSync:避免当调用cudaEventSynchronize时阻塞调用线程的情况发生,同时减少CPU周期的消耗

cudaEventDisableTiming:表示事件只用来同步，不需要记录时间

cudaEventInterprocess：事件被用于进程间事件

更具流的并发性可以分为两种调度执行：深度优先，广度优先

深度优先：每个流按照核函数的执行顺序进行注册

广度优先:按照核函数的顺序，每次往多个流中同时注册

广度优先在单一流架构上的内核执行效率要高于深度优先，但是在多流架构上会慢与深度优先，特别是存在数据传输的流上，数据传输会竞争GPU内部的数据传输队列，从而影响性能。

OpenMP可以利用编译指令识别并行区域，减少代码的书写,见stream_block.cu:openmp

利用编译参数-Xcompiler -fopenmp xxx -lgomp

可以通过设置环境变量export CUDA_DEVICE_MAX_CONNECTIONS=32来改变最大的硬件并发连接数，对于用不到那么多流的程序可以适当的降低连接数，减少资源消耗。

而对于实际的硬件允许的最大连接数(流数)，不仅受这个参数的影响还能被硬件资源所影响。

如果核函数所启动的线程过多，那么能启动的流数会变少

GPU硬件中有两个传输队列，分别是主机拷贝数据到设备，设备复制数据到主机。这两个队列可以并行操作。并且仅在这两个操作处于不同流中时会发生，见下图

流回调是在流的操作中插入一个回调函数，等到回调函数之前的操作完成后，返回调用主机上的一个函数。可以往任意流中插入主机逻辑。也可以进行GPU和CPU的同步

在流中利用函数cudaStreamAddCallback(cudaStream_t,cudaStreamCallback_t,void * userData,flags)在流中添加回调，同时阻塞后面的操作。直到回调完成

1. cudaStreamCallback_t：流回调函数
2. userData:传给回调函数的参数
3. flag:0 未使用

回调函数的定义：在定义方法之后返回类型的后面添加 CUDART_CB标识符，表示这是一个回调函数,同时添加三个形参：流对象，流错误，自定义数据

void CUDART_CB function(cudaStream_t stream,cudaError_t status,void * data)

回调函数的限制：

1. 回调函数中不可以调用CUDA的API函数
2. 回调函数中不可以执行同步

cuda编译器对自身的一些特殊的计算操作进行了指令集的优化，相比较与标准函数，cuda内部函数拥有更少的指令，同时精度也会有所缺失，但在执行效率上要远高于标准函数。例如exp,sqrt,sin,cos等标准函数都有对应的内部函数__fsqrt_xx,__expf,__sinf,__cosf。甚至包括加减乘除等基本运算也有对应的优化内部函数，具体函数可以参考头文件device_functions.h

cuda同时提供了单精度和双精度相关的函数以及计算操作。

对于双精度而言，拥有更高的计算精度，但是拥有两倍与单精度的带宽需求。在执行效率上要慢与单精度，单精度精度不够但是符合硬件读取效率。拥有较高的执行效率

参见头文件device_atomic_functions.h,里面提供了原子操作的加(atomicAdd)，减(atomicSub)，交换(atomicExch)以及CAS(atomicCAS)

原子操作确保这些操作是线程安全的，能够保证多个线程进行数据处理时数据是正确的

重要指令的介绍：

函数atomicExch(int *M,int v)用于无条件的使用V替换M的值

CAS(compare and swap) 比较并交换，该函数主要过程是，利用一个旧值与实际值相比较，如果相等则更新成期望值。如果不等则返回不做操作。

atomicCAS(int *address,int compare,int val)函数不管是否执行成功都会返回内存指向的值，函数执行成功返回的是新值，如果执行失败返回的是旧值。因此我们可以检测旧值与返回值，来判断是否执行成功。同时利用自旋来确保新值更新到了内存。

CAS操作确保了在修改值期间，没有别的线程操作内存的数据。确保内存的数据更新能被其他线程可见以此来保证数据安全

注：原子操作指令会严重影响性能,但是会确保线程安全的操作数据

利用编译参数--ptx可以生成指令集的代码

cuda编译器提供了两种指令集优化的类型：

1. 内部函数
2. 编译参数

乘后加($a\times b +c$)指令是比较常见的浮点计算指令，可以利用编译参数--fmad=true或者--fmad=false来开启或者关闭乘后加指令。默认情况下乘后加是开启的。同样的开启是以精度为代价

利用乘法的内部函数__fmul替代乘法操作符*，会阻止乘后加的优化。

--prec-div=[true | false] 用来提高单精度除法以及倒数的精度，默认为true,但是会降低性能

--prec-sqrt = [true | false]，提高开方函数的精度，默认为true

--use_fast_math 用内部函数替换调所有标准函数，同时设置--prec-div = false以及--prec-sqrt = false

很多内部函数都带有不同的后缀:rn,rz,ru,rd

1. rn:近似值(默认)
2. rz:向零取整
3. ru:向下取整
4. rd:向上取整

大部分原子操作都只支持int,unsigned int,unsigned long long int相关类型的，对于浮点数，我们可以使用unsigned int作为中介进行转换，例如__float2uint_rn等方法

主要介绍三方库的使用，cufft/cublas

1. 创建句柄
2. 分配数据内存，并拷贝数据
3. 更具数据排列（列优先，行优先）配置函数库
4. 执行函数
5. 取回数据
6. 释放句柄

稀疏矩阵计算库,主要用于稀疏矩阵的计算，例如稀疏矩阵与矩阵乘法/稀疏矩阵向量乘法等

主要难点在于稀疏矩阵的存储格式转换，稀疏矩阵有如下几个坐标存储方式

1. COO(坐标存储方式),按照数据的行列坐标以及值进行存储，组成一个列表
2. CSR(压缩稀疏行方式),只需要知道偏移量和行的长度就能直到一行数据有哪些。这三个单独的索引数组如下：
   1. V：按照每行非零元素压缩存储
   2. C：列的索引
   3. R：R[i] = V，C中第i行的偏移量
3. CSC(压缩稀疏列方式)，与CSR大致相同，不过输入数据是按列存储，也是按列进行压缩。比CSR更省空间

线性代数计算库,主要用于线性代数的相关计算，例如向量点乘，矩阵向量乘法，矩阵与矩阵乘法，矩阵转置以及乘加算法等，具体参考相关API

主要伪代码如下：

cublasHandle_t handle;
cublasCreateHandle()
/*****Data copy to Device ****/
...

/*****调用计算，注意参数设置****/
cublas<T>gemm()

/*****释放资源*****/
cublasDestroy()

主要用到的函数有：

矩阵乘法:cublas<T>gemm,T in {S/s,D/d,C/c,Z/z} 分别对应数据类型float,double,cuFloatComplex,cuDoubleComplex,

该函数计算表达式为：$C = \alpha op(A)op(B) + \beta C$

op()表示对矩阵的操作：不操作，转置，海森矩阵

$\alpha,\beta$均为缩放系数，$\beta$为0时C表示为无效输入

leading dimension:which in the case of column-major storage is the number of rows of the allocated matrix ,主维度,按列主序存储时主维度为行数，同样按行主序存储时为列数。

矩阵求和：cublasCgeam

该函数计算表达式子：$C = \alpha op(A) + \beta op(B)$

同样的该函数可以用来求矩阵的转置，相比较与gemm只需要两个或者一个数组参与计算

矩阵求逆：cublas<T>matinvBatched，该方法为批量矩阵求逆，可以同时进行N个批次的矩阵求逆，但该方法仅支持32x32大小以下的矩阵求逆，

除了利用该方法，还能使用LU分解求逆:cublas<T>getrfBatched/cublas<T>getriBatched.先调用cublasgetrfBatched 进行LU(L:下三角单位对角矩阵，U：上三角矩阵)因式分解，然后调用 cublas<T>getriBatched输入LU得到的P，计算得到逆矩阵。

// step 1: perform in-place LU decomposition, P*A = L*U.
//      Aarray[i] is n*n matrix A[i]
    cublasDgetrfBatched(handle, n, Aarray, lda, PivotArray, infoArray, batchSize);
//      check infoArray[i] to see if factorization of A[i] is successful or not.
//      Array[i] contains LU factorization of A[i]

// step 2: perform out-of-place inversion, Carray[i] = inv(A[i])
    cublasDgetriBatched(handle, n, Aarray, lda, PivotArray, Carray, ldc, infoArray, batchSize);
//      check infoArray[i] to see if inversion of A[i] is successful or not.

方法二不限矩阵大小。

cublasSetVector/cublasSetMatrix，向量,矩阵复制函数，将主机的数据拷贝到设备上

傅里叶变换库

主要计算流程：

createHandle
cufftPlanMany()//设置计算计划
cufftExecC2C() // 执行
DestroyHandle

傅里叶变换主要在于PlanMany的设置用于设置多个fft变换 cufftPlanMany参数说明：

1. plan:cufftHandle指针
2. int rank:几维傅里叶变换
3. int * n:与rank值同维度的数组，用于描述每个fft变换有多少个元素，n[0]最外面的维度，n[rank-1]最里面的维度
4. int * inembed:输入数据的内存排布/存储维度，按行优先存储则inembed[0]=col，inembed[1]=row，按列优先存储则inembed[0]=row,inembed[1]=col
5. int istride:同一个fft变换中元素的距离
6. int idist:相邻两个fft变换中元素的距离
7. int * onembed:输出数据的内存排列
8. int ostride:同一个fft变换中元素的距离
9. int odist:相邻两个fft变换中元素的距离
10. cufftType type:fft变换类型C2C,C2R,R2C,R2R：复数到复数，复数到实数，实数到复数，实数到实数

更具type分别调用对应的执行函数

同样cufft支持多流的处理，利用cufftSetStream设置handle是处于哪个流

注：fft的执行函数必须在所有planMany之后，否则会报错CUFFT_INTERNAL_ERROR

随机数生成库

多GPU主要包含任务分配以及多GPU之间的数据传输，为了避免使用主机内存中转，cuda提供了一套点对点传输的API，用于在多GPU之间直接通过PCIe进行数据传输。

在显示调用GPU操作前，调用cudaSetDevice可以指定设备，一旦选定设备，接下来的操作都会在这个设备上。并且，该函数并不会导致主机同步。可以很快的将控制权返回给主机。

获取设备信息：可以通过cudaGetDeviceCount来获取可用GPU的数量，同时更具cudaGetDeviceProperties来获取设备信息

首先对于任意设备都需要开启点对点的功能，因为不是所有设备都支持这个功能，还要显示的查询设备是否支持这个功能，利用cudaDeviceCanAccessPeer函数来查询设备支不支持点对点传输， 然后通过cudaDeviceEnablePeerAccess(int peerDevice)来设定当前设备对peerDevice设备可进行传输。

这个功能会显示开启直到显示禁用cudaDeviceDisablePeerAccess

数据复制函数则是 cudaMemcpyPeerAsync用来将数据复制到另一个设备上