CUDA Libraries简单介绍
上图是CUDA 库的位置。本文简要介绍cuSPARSE、cuBLAS、cuFFT和cuRAND。之后会介绍OpenACC。
- cuSPARSE线性代数库，主要针对稀疏矩阵之类的。
- cuBLAS是CUDA标准的线代库，只是没有专门针对稀疏矩阵的操作。
- cuFFT傅里叶变换
- cuRAND随机数
CUDA库和CPU编程所用到的库没有什么区别，都是一系列接口的集合，主要优点是，仅仅须要编写host代码，调用相应API就可以，能够节约非常多开发时间。并且我们全然能够信任这些库能够达到非常好的性能，写这些库的人都是在CUDA上的大能。一般人比不了。当然。全然依赖于这些库而对CUDA性能优化一无所知也是不行的，我们依旧须要手动做一些改进来挖掘出更好的性能。
下图是《CUDA C编程》中提到的一些支持的库。详细细节能够在NVIDIA开发人员论坛查看：
假设大家的APP属于上面库的应用范围，非常建议大家使用。
A Common Library Workflow
以下是一个使用CUDA库的详细步骤，当然，各个库的使用可能不尽同样，可是不会逃脱以下的几个步骤，差异基本上就是少了哪几步而已。
- 创建一个库的句柄来管理上下文信息。
- 分配device存储空间给输入输出。
- 假设输入的格式并非库中API支持的须要做一下转换。
- 填充device Memory数据。
- 配置library computation以便运行。
- 调用库函数来让GPU工作。
- 取回device Memory中的结果。
- 假设取回的结果不是APP的原始格式。就做一次转换。
- 释放CUDA资源。
- 继续其他的工作。
以下是这几个步骤的一些细节解释：
Stage1：Creating a Library Handle
CUDA库好多都有一个handle的概念，其包括了该库的一些上下文信息。比方数据格式、device的使用等。对于使用handle的库，我们第一步就是初始化这么一个东西。一般的。我们能够觉得这是一个存放在host对程序员透明的object，这个object包括了跟这个库相关联的一些信息。比如。我们可定希望全部的库的操作运行在一个特别的CUDA stream，虽然不同的库使用不同函数名字。可是大多数都会规定全部的库操作以一定的stream发生（比方cuSPARSE使用cusparseSetSStream、cuBLAS使用cublasSetStream、cuFFT使用cufftSetStream）。
stream的信息就会保存在这个handle中。
Stage2：Allocating Device Memory
本文所讲的库。其device存储空间的分配依旧是cudaMalloc或者库自己调用cudaMalloc。
仅仅有在使用多GPU编程的库时，才会使用一些定制的API来实现内存分配。
Stage3：Converting Inputs to a Library-Supported Format
假设APP的数据格式和库要求的输入格式不同的话。就须要做一次转化。比方，我们APP存储一个row-major的2D数组，可是库却要求一个column-major，这就须要做一次转换了。为了最优性能。我们应该尽量避免这样的转化，也就是尽量和库的格式保持一致。
Stage4：Populating Device Memory with Inputs
完毕上述三步后，就是将host的数据传送到device了，也就是相似cudaMemcpy的作用。之所说相似。是引文大部分库都有自己的API来实现这个功能，而不是直接调用cudaMemcpy。
比如。当使用cuBLAS的时候，我们要将一个vector传送到device，使用的就是cubalsSetVector，当然其内部还是调用了cudaMemcpy或者其他等价函数来实现传输。
Stage5：Configuring the Library
有步骤3知道。数据格式是个明显的问题。库函数须要知道自己应该使用什么数据格式。某些情况下，相似数据维度之类的数据格式信息会直接当做函数參数配置。其他的情形下。就须要手动来配置下之前说的库的handle了。还有个别情况是，我们须要管理一些分离的元数据对象。
Stage6：Executing
运行就简单多了，做好之前的步骤。配置好參数。直接调用库API。
Stage7：Retrieving Results from Device Memory
这一步将计算结果从device送回host，当然还是须要注意数据格式，这一步就是步骤4的反过程。
Stage8：Converting Back to Native Format
假设计算结果和APP的原始数据格式不同，就须要做一次转化，这一步是步骤3的反过程。
Stage9：Releasing CUDA Resources
假设上面步骤使用的内存资源不再使用就须要释放掉，正如我们曾经介绍的那样，内存的分配和释放是非常大的负担。所以希望尽可能的资源重用。比方device Memory、handles和CUDA stream这些资源。
Stage10：Continuing with the Application
继续干别的。
再次重申，上面的步骤可能会给大家使用库是非常麻烦低效的事儿。但事实上这些步骤通常是冗余的，非常多情况下，当中的非常多步骤是不必要的，在以下的章节我们会介绍几个基本的库以及其简要使用。相信看过后。你就不会觉得使用库得不偿失了。
THE CUSPARSE LIBRARY
cuSPARSE就是一个线性代数库。对稀疏矩阵之类的操作尤其独到的使用方法。使用非常宽泛。
他当对稠密和稀疏的数据格式都支持。
下图是该库的一些函数调用。从中能够对其功能有一个大致的了解。
cuSPARSE将函数以level区分，全部level 1的function仅操作稠密和稀疏的vector。
全部level2函数操作稀疏矩阵和稠密vector。
全部level3函数操作稀疏和稠密矩阵。
cuSPARSE Data Storage Formats
稠密矩阵就是当中的值大部分非零。
稠密矩阵全部值都是存储在一个多维的数组中的。
相对而言，稀疏矩阵和vector中元素主要是零。所以其存储就能够做一些文章。比方我们能够仅仅保存非零值和其坐标。cuSPARSE支持非常多种稀疏矩阵的存储方式。本文仅仅介绍当中三种。
先看一下稠密（dens）矩阵的存储方式，图示非常明显。不多说了：
Coordinate（COO）
对于稀疏矩阵中的每一个非零值，COO方式都保存其行和列坐标，因此。当通过行列检索矩阵值的时候，假设该行列值没有在存储格式中匹配到的话。必定就是零了。
我们应该注意到了。所谓稀疏矩阵要稀疏到什么程度才干使用COO呢？这个须要详细问题详细分析了，主要跟元素数据类型和索引数据类型有关。比方，一个存储32-bit的浮点类型数据的稀疏矩阵，索引使用32-bit的整型格式。那么仅仅有当非零数据少于于矩阵的三分之中的一个的时候才会节约存储空间。
Compressed Sparse Row（CSR）
CSR和COO相似，唯一不同就是非零值的行索引。
COO模式下，全部非零值都会相应一个int的行索引，而CSR则是存储一个偏移值，这个偏移值是全部属于同一行的值拥有的属性。
例如以下图所看到的，相比COO，降低了row：
由于全部存储在同一行的数据在内存中是相邻的。要找到某一行相应的值仅仅须要一个偏移量和length。
比如，假设仅仅想知道第三行的非零值，我们能够使用偏移量为2，length为2在V中检索，例如以下图所看到的：
对图中的C使用同样的偏移和length就能定位列索引，也就能全然确定一个value在矩阵中的位置。当存储一个非常大的矩阵且相对来说每行数据都非常多的时候，使用CSR比存储每一个非零值的索引要有效得多。
如今我们要考虑这些偏移地址和length的存储了。最简单的方式是创建两个数组Ro和Rl，每一个都相应一个nRows用作length。假设矩阵有大量的行就须要分配两个非常大的数组。鉴于此，我们能够使用单独的一个length为nRows+1的数组R，第i行的偏移地址就存储在R[i]。第i行的长度能够通过比較R[I+1]和R[i]值来做出推断，还有就是R[i+1]是用来存储矩阵非零值的总数的。本例中R数组例如以下：
由上图知。0行的偏移地址是0,1行偏移地址是1,2行偏移地址是2，共同拥有4个非零元素，我们能够找矩阵行为0的值及其列索引。由于R[1]-R[0]=1-0=1，说明第一行仅有一个非零值，其列索引为0，其值为3。
这样。对于每行都有多个非零值的稀疏矩阵存储，CSR比COO要节约空间。下图是CSR的完整示意图：
使用CSR格式稀疏矩阵的function非常直观。首先，我们在host定义一个CSR格式的稀疏矩阵，其代码例如以下：
float *h_csrVals;
int *h_csrCols;
int *h_csrRows;
h_csrVals用来存储非零值个数。h_csrCols存储列索引，h_csrRows存储行偏移。接下来就是分配device内存之类的常规操作：
cudaMalloc((void **)&d_csrVals, n_vals * sizeof(float));
cudaMalloc((void **)&d_csrCols, n_vals * sizeof(int));
cudaMalloc((void **)&d_csrRows, (n_rows + 1) * sizeof(int));
cudaMemcpy(d_csrVals, h_csrVals, n_vals * sizeof(float),cudaMemcpyHostToDevice);
cudaMemcpy(d_csrCols, h_csrCols, n_vals * sizeof(int),cudaMemcpyHostToDevice)；
cudaMemcpy(d_csrRows, h_csrRows, (n_rows + 1) * sizeof(int),cudaMemcpyHostToDevice);
上述三种（包括稠密矩阵）数据格式各有各擅长的方面。下图列出了cuSPARSE支持的一些数据格式以及各自的最佳使用场景：
Formatting Conversion with cuSPARSE
从前文可知，这个过程应该尽量避免。转换不仅须要有计算的开销，还有额外存储的空间浪费。还有就是在使用cuSPARSE也应该尽量发挥其在稀疏矩阵存储上的优势。由于好多APP的latency就是仅仅简单的使用稠密矩阵存储方式。由于cuSPARSE的数据格式众多，其用来转换的API也不少。下图列出了这些转换API。
左边那一列是你要转换的目标格式。为空表示不支持两种数据格式的转换，虽然如此，你还能够通过多次转换来实现未显示支持的转换API，比方dense2bsr没有被支持，可是我们能够使用dense2csr和csr2bsr两个过程来达到目的。
Demonstrating cuSPARSE
这部分演示样例代码会涉及到矩阵向量相乘，数据格式转换，以及其他cuSPARSE的特征。
// Create the cuSPARSE handle
cusparseCreate(&handle);
// Allocate device memory for vectors and the dense form of the matrix A
...
// Construct a descriptor of the matrix A
cusparseCreateMatDescr(&descr);
cusparseSetMatType(descr, CUSPARSE_MATRIX_TYPE_GENERAL);
cusparseSetMatIndexBase(descr, CUSPARSE_INDEX_BASE_ZERO);
// Transfer the input vectors and dense matrix A to the device
...
// Compute the number of non-zero elements in A
cusparseSnnz(handle, CUSPARSE_DIRECTION_ROW, M, N, descr, dA,M, dNnzPerRow, &totalNnz);
// Allocate device memory to store the sparse CSR representation of A
...
// Convert A from a dense formatting to a CSR formatting, using the GPU
cusparseSdense2csr(handle, M, N, descr, dA, M, dNnzPerRow,dCsrValA, dCsrRowPtrA, dCsrColIndA);
// Perform matrix-vector multiplication with the CSR-formatted matrix A
cusparseScsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,M, N, totalNnz, &alpha, descr, dCsrValA, dCsrRowPtrA,dCsrColIndA, dX, &beta, dY);
// Copy the result vector back to the host
cudaMemcpy(Y, dY, sizeof(float) * M, cudaMemcpyDeviceToHost);
上述代码的过程能够总结为：
- 使用cusparseCreate创建库的handle。
- 使用cudaMalloc分配device内存空间用来存储矩阵和向量，并分别使用dense和CSR两种格式存储。
- cusparseCreateMatDescr和cusparseSetMat*使用来配置矩阵属性的，cudaMemcpy用来拷贝数据到device，cusparseSdense2csr来生成CSR格式的数据，非零数据个数有cusparseSnnz计算得到。
- cusparseScsrmv是矩阵和向量乘操作函数。
- 再次使用cudaMemcpy将结果拷贝回host。
- 释放资源。
编译：
$ nvcc -lcusparse cusparse.cu –o cusparse
Important Topics in cuSPARSE Development
虽然cuSPARSE提供了一个相对来说最高速和简洁的方式以达到高性能的线性代数库，我们仍须要谨记cuSPARSE使用的一些关键点。
第一点就是，要保证正确的矩阵和向量的数据格式，cuSPARSE本身没有什么能力来检測出错误的或者不恰当的数据格式，而一次错误的格式操作就可能导致段错误，这也算是给自己debug提供一种方向吧，虽然段错误五花八门。
对于矩阵和向量规模比較小的情况下，手动验证其数据格式还是可行的。
我们能够将转换后的数据进行一次逆转换过程来和原始数据比对。
第二点是cuSPARSE的默认异步行为。当然这对于GPU编程来说。已经习以为常了，可是对于传统的host端堵塞式的数学库来说，GPU的计算结果会非常有趣。
对于cuSPARSE来说，假设使用了cudaMemcpy拷贝数据后。host会自己主动堵塞住。等待device的计算结果。可是假设cuSPARSE库被配置来使用CUDA steam和cudaMemcpyAsync，我们就须要多留一个心眼，使用确保正确的同步行为来获取device的计算结果。
最后一点比較新鲜的是标量的使用，这里要使用标量的引用形式。例如以下代码中的beta变量：
float beta = 4.0f;
...
// Perform matrix-vector multiplication with the CSR-formatted matrix A
cusparseScsrmv(handle, CUSPARSE_OPERATION_NON_TRANSPOSE,
M, N, totalNnz, &alpha, descr, dCsrValA, dCsrRowPtrA,
dCsrColIndA, dX, &beta, dY);
假设不小心直接传递了beta这个參数。APP会报错（SEGFAULT）。不注意的话这个bug非常不好查。除此外。当标量作为输出參数时，能够使用指针。cuSPARSE提供了cusparseSetPointMode这个API来调整是否使用指针来获取计算结果。
THE cuBLAS LIBRARY
cuBLAS也是一个线代库。不同于cuSPARSE，cuBLAS传统线代库的接口。BLAS即Basic Linear Algebra Subprograms的意思。cuBLAS level1是专门的vector之间操作。level2是矩阵和向量之间的操作。
level3是矩阵和矩阵之间的操作。相对于cuSPARSE。cuBLAS不支持稀疏矩阵数据格式，他仅仅支持并且善于稠密矩阵和向量的使用。
由于BLAS库最初是由FORTRAN语言编写。他就是用了column-major和one-based的方式存储数据，而cuSPARSE则是使用的row-major。
下图是这样的方式的存储格式。一看便明：
我们能够比較下row-major和column-major将二维转化为一维的过程公式：
为了考虑兼容性，cuBLAS也使用了column-major的方式存储。所以，对于习惯C/C++的人来说，这可能比較让人困惑吧。
还有一方面。就像C和其他语言那样，one-based索引意味着数组中第一个元素的引用使用1而不是0，也就是说。一个有N个元素的数组，其最后一个值的索引是N而不是N-1。
可是，cuBLAS没有办法决定C/C++（cuBLAS使用C/C++）编程的语境，所以他就必须使用zero-based索引，这就导致了一个奇怪的混乱情况，要满足FORTRAN的column-major，但one-based却不行。
cuBLAS提出了两个API，cuBLASASLegacy API是cuBLAS最開始的一个实现。已经废弃。当前使用cuBLAS API。二者差异非常小。
看过接下来的内容，你会发现。cuBLAS的使用流程跟cuSPARSE有非常多同样之处，所以对于这些库代码编写基本能够触类旁通。
Managing cuBLAS Data
相较于cuSPARSE，cuBLAS的数据格式要简单的多，全部操作都作用在稠密向量或矩阵。同样是使用cudaMalloc来分配device内存空间，可是使用cublasSetVector/cublasGetVector和cubalsSetMartix/cublasGetMartix在device和host之间传送数据（事实上相对cuSPARSE也没多大区别）。本质上，这些API底层都是调用cudaMemcpy，并且，他们对Strided和unstrided数据都有非常好的优化，比方以下的代码：
cublasStatus_t cublasSetMatrix(int rows, int cols, int elementSize,const void *A, int lda, void *B, int ldb);
这些參数大部分看名字就知道什么意思了，当中lda和ldb指明了源矩阵A和目的矩阵B的主维度（leading dimension）。所谓主维就是矩阵的行总数，这个參数仅仅在须要host矩阵一部分数据的时候非常实用。也就是说。当须要完整的矩阵时，lda和ldb都应该是M。
假设我们使用一个稠密的二维column-major的矩阵A。其元素是单精度浮点类型，矩阵大小为MxN，则使用以下的函数传输矩阵：
cublasSetMatrix(M, N, sizeof(float), A, M, dA, M);
也能够例如以下传输一个仅仅有一列的矩阵A到一个向量dV：
cublasStatus_t cublasSetVector(int n, int elemSize, const void *x, int incx,void *y, int incy)；
x是host上源起始地址，y是device上目的起始地址，n是要传送数据的总数，elemSize是每一个元素的大小，单位是byte。incx/incy是要传送的元素之间地址间隔。或者叫步调。传送一个单列长度M的column-major 矩阵A到向量dV例如以下：
cublasSetVector(M, sizeof(float), A, 1, dV, 1);
也能够例如以下传送一个单行的矩阵A到一个向量dV：
cublasSetVector(N, sizeof(float), A, M, dV, 1);
通过这些样例能够发现，使用cuBLAS要比cuSPARSEeasy的多。所以除非我们的APP对稀疏矩阵需求比較大，一般都是用cuBLAS，保证性能的同一时候，还能提高开发效率。
Demonstrating cuBLAS
这部分代码主要关注cuBLAS的一些统一使用并理解他为什么易于使用。得益于GPU的高性能计算，表现要比在CPU上的BLAS号15倍，并且cuBLAS的开发也就比传统的BLAS略微费事儿。
// Create the cuBLAS handle
cublasCreate(&handle);
// Allocate device memory
cudaMalloc((void **)&dA, sizeof(float) * M * N);
cudaMalloc((void **)&dX, sizeof(float) * N);
cudaMalloc((void **)&dY, sizeof(float) * M);
// Transfer inputs to the device
cublasSetVector(N, sizeof(float), X, 1, dX, 1);
cublasSetVector(M, sizeof(float), Y, 1, dY, 1);
cublasSetMatrix(M, N, sizeof(float), A, M, dA, M);
// Execute the matrix-vector multiplication
cublasSgemv(handle, CUBLAS_OP_N, M, N, &alpha, dA, M, dX, 1,&beta, dY, 1);
// Retrieve the output vector from the device
cublasGetVector(M, sizeof(float), dY, 1, Y, 1);
使用cuBLAS比較直观，易于理解，其使用步骤基本例如以下：
- 使用cublasCreateHandle创建handle。
- 使用cudaMalloc分配device资源。
- 使用cublasSetVector和cublasSetMartix向device填充数据。
- 使用cublasSgemv运行矩阵和向量的乘操作。
- 使用cublasGetVector获取计算结果。
- 释放资源。
编译命令：
$ nvcc -lcublas cublas.cu
Porting from BLAS
将一个传统的C实现的APP（使用BLAS库）转化为cuBLAS也是比較直观的，基本能够归纳为以下几步：
- 添加device Memory的分配操作（cudaMalloc）和其资源释放操作。
- 添加device和host之间数据传送的过程。
- 变更BLAS的API为等价的cuBLAS API。这一步比較麻烦，这里以之前的代码为例：
// Allocate device memory
cudaMalloc((void **)&dA, sizeof(float) * M * N);
cudaMalloc((void **)&dX, sizeof(float) * N);
cudaMalloc((void **)&dY, sizeof(float) * M);
// Transfer inputs to the device
cublasSetVector(N, sizeof(float), X, 1, dX, 1);
cublasSetVector(M, sizeof(float), Y, 1, dY, 1);
cublasSetMatrix(M, N, sizeof(float), A, M, dA, M);
// Execute the matrix-vector multiplication
cublasSgemv(handle, CUBLAS_OP_N, M, N, &alpha, dA, M, dX, 1,&beta, dY, 1);
// Retrieve the output vector from the device
cublasGetVector(M, sizeof(float), dY, 1, Y, 1);
其等价的BLAS代码是：
void cblas_sgemv(const CBLAS_ORDER order, const CBLAS_TRANSPOSE TransA,const MKL_INT M, const MKL_INT N, const float alpha, const float *A,const MKL_INT lda, const float *X, const MKL_INT incX, const float beta, float *Y,const MKL_INT incY);
二者还是有非常多相似之处的。不同的是。BLAS有个order參数来使用户能够指定输入数据是row-major还是column-major。
还有就是BLAS的beta和alpha没有使用引用形式。
4. 最后就是在实现功能后调节性能了，比方：
- 复用device资源而不是释放。
- device和host之间传输数据尽量降低冗余数据。
- 使用stream-based运行来实现异步传输。
Important Topics in cuBLAS Development
相较于cuSPARSE。假设大家对BLAS熟悉的话，cuBLAS更easy上手。只是须要注意的是，虽然cuBLAS的行为更easy理解，可是有时候恰恰是这份理所当然的理解会造成一些认识误区，毕竟cuBLAS并不等于BLAS。
对于大部分习惯于row-major的编程语言，使用cuBLAS就得分外小心了，我们可能非常熟悉将一个row-major的多维数组展开，可是过度到column-major会有点不适应，以下的宏定义能够帮我们实现row-major到column-major的转换：
#define R2C(r, c, nrows) ((c) * (nrows) + (r))
只是，当使用上述的宏时，仍然须要一些循环的顺序问题，对于C/C++程序员来说，会经经常使用以下的代码：
for (int r = 0; r < nrows; r++) {
for (int c = 0; c < ncols; c++) {
A[R2C(r, c, nrows)] = ...
}
}
代码当然没什么问题。可是却不是最优的，由于他在訪问A的时候，不是线性扫描内存空间的。
假设nrows非常大的话。cache命中率基本为零了。因此。我们须要以下这样的代码：
for (int c = 0; c < ncols; c++) {
for (int r = 0; r < nrows; r++) {
A[R2C(r, c, nrows)] = ...
}
}
所以，做优化要步步小心，由于一个没注意。就有可能导致非常差的cache命中。