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要比cuSPARSE容易的多，所以除非我們的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更容易上手。不過需要注意的是，雖然cuBLAS的行為更容易理解，但是有時候恰恰是這份理所當然的理解會造成一些認識誤區，畢竟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命中。
cuFFT
未完待續~~~
參考書：《professional cuda c programming》