1. CUDA學習步驟

1. CPU實現 a*b = c 的矩陣乘法（矩陣尺寸是n*m的，n和m大於1000）
2. 將cpu程式碼移植到cuda。將CPU值傳入GPU，使用cuda計算，與cpu結果對比。
3. 優化思路1：將矩陣分塊進行計算
4. 優化思路2：使用share memory進行優化
5. 優化思路3：將資料繫結在texture上

2. CPU實現的矩陣乘法

廢話不多說，直接上原始碼

/* CPUMatMultiply：CPU下矩陣乘法
 * a:第一個矩陣指標，表示a[M][N]
 * b:第二個矩陣指標，表示b[N][S]
 * result:結果矩陣，表示為result[M][S]
 */
void CPUMatMultiply(const int * a,const int *
 
 b, int *result,const int M,const int N,const int S)
{
    for (int i = 0; i < M; i++)
    {
        for (int j = 0; j < S; j++)
        {
            int index = i * S + j;
            result[index] = 0;

            //迴圈計算每一個元素的結果
            for (int k = 0; k < N; k++)
            {
                result[index] += a[i * N + k] * b[k * S + j];
 

            }
        }
    }
}

3. CUDA實現的矩陣乘法

下面直接進入正題，矩陣乘法的移植。
從CPU上直接移植矩陣乘法到GPU上是非常簡單的，不需要for迴圈，直接通過CUDA執行緒的id號，即threadIdx.x和threadIdx.y即可操作相應的資料。

gpu矩陣乘法核函式-原始碼：

/* gpuMatMultKernel：GPU下矩陣乘法核函式
*  a:第一個矩陣指標，表示a[M][N]
*  b:第二個矩陣指標，表示b[N][S]
*  result:結果矩陣，表示result[M][S]
*/
__global__ void gpuMatMultKernel(const int *a, const int *
 
b, int *result, const int M, const int N, const int S)
{
    int threadId = (blockIdx.y * blockDim.y + threadIdx.y) * gridDim.x * blockDim.x 
                    + blockIdx.x * blockDim.x + threadIdx.x;
    if (threadId < M * S)
    {
        int row = threadId / S;
        int column = threadId % S;

        result[threadId] = 0;
        for (int i = 0; i < N; i++)
        {
            result[threadId] += a[row * N + i] * b[i * S + column];
        }
    }
}

其中，blockIdx、blockDim、threadIdx、gridDim都是CUDA的內建變數。
blockIdx、threadIdx：分別CUDA中執行緒塊的ID、執行緒的ID。
blockDim、gridDim：分別是CUDA中執行緒塊的維度，執行緒網格的維度。
if語句是為了判斷是否是對矩陣中的資料進行操作，防止在當執行緒threadId超過了M*S時，使用result[threadId]出現越界行為。
後續的操作則是矩陣乘法的簡單實現。

4. 程式優化1：共享記憶體分塊運算

CUDA C支援共享記憶體。關鍵字__shared__新增到宣告中，這將使這個變數駐留在共享記憶體中。當宣告變數為shared變數時，執行緒塊中的每一個執行緒都共享這塊記憶體，使得一個執行緒塊中的多個執行緒能夠在計算上進行通訊和協作。而且，共享記憶體緩衝區駐留在物理GPU上，而不是駐留在GPU之外的系統記憶體中。因此，在訪問共享記憶體時的延遲要遠遠低於訪問普通緩衝區的延遲，提高了效率。

使用共享記憶體的核函式：

/* gpuMatMultWithSharedKernel：GPU下使用shared記憶體的矩陣乘法
*  a:第一個矩陣指標，表示a[height_A][width_A]
*  b:第二個矩陣指標，表示b[width_A][width_B]
*  result:結果矩陣，表示result[height_A][width_B]
*/
template<int BLOCK_SIZE>
__global__ void gpuMatMultWithSharedKernel(const int *a, const int *b, int *result, const int height_A, const int width_A, const int width_B)
{
    int block_x = blockIdx.x;
    int block_y = blockIdx.y;
    int thread_x = threadIdx.x;
    int thread_y = threadIdx.y;

    if ((thread_y + block_y * blockDim.y) * width_B + block_x * blockDim.x + thread_x >= height_A * width_B)
    {
        return;
    }

    const int begin_a = block_y * blockDim.y * width_A;
    const int end_a = begin_a + width_A - 1;
    const int step_a = blockDim.x;

    const int begin_b = block_x * blockDim.x;
    const int step_b = blockDim.y * width_B;

    int result_temp = 0;

    for (int index_a = begin_a, int index_b = begin_b;
        index_a < end_a; index_a += step_a, index_b += step_b)
    {
        __shared__ int SubMat_A[BLOCK_SIZE][BLOCK_SIZE];
        __shared__ int SubMat_B[BLOCK_SIZE][BLOCK_SIZE];

        SubMat_A[thread_y][thread_x] = a[index_a + thread_y * width_A + thread_x];
        SubMat_B[thread_y][thread_x] = b[index_b + thread_y * width_B + thread_x];

        __syncthreads();

        for (int i = 0; i < BLOCK_SIZE; i++)
        {
            result_temp += SubMat_A[thread_y][i] * SubMat_B[i][thread_x];
        }

        __syncthreads();
    }

    int begin_result = block_y * blockDim.y * width_B + begin_b;
    result[begin_result + thread_y * width_B + thread_x] = result_temp;
}

矩陣乘法的並行運算，每次計算矩陣的一塊資料。利用共享記憶體的共享功能，每次將一塊資料儲存到共享記憶體中使得一個執行緒塊同時呼叫資料進行計算當前塊相對應得矩陣乘法結果值。

• 程式碼 __shared__ int SubMat_A中的__shared__宣告變數SubMat_A為共享記憶體中儲存的變數。然後將陣列中的資料提取到變數SubMat_A中儲存在共享記憶體。
• __syncthreads()對執行緒塊中的執行緒進行同步，確保對__shared__進行下面的操作時上面的操作已經完成。
• 兩個for迴圈完成了當前執行緒塊對應矩陣子塊的乘法結果計算。

注意：CUDA架構將確保，除非執行緒塊中的每個執行緒都執行了__syncthreads()，否則沒有任何執行緒能執行__syncthreads()之後的指令。例如：

if (cacheIndex < i){
    cache[cacheIndex] += cache[i];
    __syncthreads();
}

此時，如果有執行緒中的cacheIndex>i則會出現一直等待的情況。

5. 程式優化2：紋理記憶體的運用

紋理記憶體，CUDA C程式中的另一種只讀記憶體。紋理記憶體是專門為那些在記憶體訪問模式中存在大量空間區域性性的圖形應用程式而設計的。在某個計算應用程式中，這意味著一個執行緒讀取的位置可能與鄰近執行緒讀取的位置“非常接近”，如下圖所示。

從數學的角度看，這四個地址並非連續的，在一般的CPU快取模式中，這些地址將不會快取。但由於GPU紋理快取是專門為了加速這種訪問模式而設計的，因此如果在這種情況中使用紋理記憶體而不是全域性記憶體，那麼將會獲得性能提升。

使用紋理記憶體的核函式-原始碼：

/* gpuMatMultWithTextureKernel：GPU下使用texture記憶體的矩陣乘法
*  result：結果矩陣，表示為result[M][S];
*  M：表示為矩陣A與矩陣result的行數
*  N：表示矩陣A的列數，矩陣B的行數
*  S：表示矩陣B和矩陣result的列數
*/
__global__ void gpuMatMultWithTextureKernel(int * result, const int M, const int N, const int S)
{
    int x = threadIdx.x + blockIdx.x * blockDim.x;
    int y = threadIdx.y + blockIdx.y * blockDim.y;
    int offset = x + y * blockDim.x * gridDim.x;

    if (offset < M * S)
    {
        int a = 0, b = 0;
        int temp_result = 0;
        for (int i = 0; i < N; i++)
        {
            a = tex1Dfetch(texA, y * N + i);
            b = tex1Dfetch(texB, i * S + x);
            temp_result += a * b;
        }
        result[offset] = temp_result;
    }
}

紋理記憶體的運用與普通記憶體運用時的演算法大致相同，只不過資料是在核函式中呼叫tex1Dfetch從紋理中提取。

在使用紋理記憶體時，主要注意的是紋理記憶體的使用。
首先，需要將輸入的資料宣告為texture型別的引用。
注意，輸入的資料是什麼型別，相應的紋理也應該與之一致。並且紋理引用必須宣告為檔案作用域內的全域性變數。

//這些變數將位於GPU上
texture<int> texA;
//二維紋理引用，增加了代表維數的引數2
texture<float, 2> texB;

在為這兩個緩衝區分配了GPU記憶體後，需要通過cudaBindTexture將這些變數繫結到記憶體緩衝區。這相當於告訴CUDA執行時兩件事：

• 我們希望將指定的緩衝區作為紋理來使用。
• 我們希望將紋理引用作為紋理的“名字”。

cudaBindTexture(NULL, texA, dev_a, desc, M * N * sizeof(int));
cudaBindTexture(NULL, texB, dev_b, desc, N * S * sizeof(int));

在繫結紋理時，CUDA執行時要求提供一個cudaChannelFormatDesc。此時，需要呼叫cudaCreateChannelDesc<int>()。

最後，通過cudaUnbindTexture()函式來取消紋理的繫結。

6.總結