cuda學習3-共享內存和同步

　　為什么要使用共享內存呢，因為共享內存的訪問速度快。這是首先要明確的，下面詳細研究。

　　cuda程序中的內存使用分為主機內存（host memory） 和 設備內存（device memory），我們在這里關注的是設備內存。設備內存都位於gpu之上，前面我們看到在計算開始之前，每次我們都要在device上申請內存空間，然后把host上的數據傳入device內存。cudaMalloc（）申請的內存，還有在核函數中用正常方法申請的變量的內存。這些內存叫做全局內存，那么還有沒有別的內存種類呢？常用的還有共享內存，常量內存，紋理內存，他們都用一些不正常的方法申請。

　　他們的申請方法如下：

　　共享內存：__shared__  變量類型 變量名；

　　常量內存：__constant__ 變量類型 變量名；

　　紋理內存：texture<變量類型> 變量名；

存儲類型	寄存器	共享內存	紋理內存	常量內存	全局內存
帶寬	～8TB/s	～1.5TB/s	～200MB/s	～200MB/s	～200MB/s
延遲	1個周期	1～32周期	400～600周期	400～600周期	400～600周期

　　他們在不同的情況下有各自的作用，他們最大的區別就是帶寬不同，通俗說就是訪問速度不同。后面三個看起來沒什么不同，但是他們在物理結構方面有差別，適用於不同的情況。

　　共享內存實際上是可受用戶控制的一級緩存。申請共享內存后，其內容在每一個用到的block被復制一遍，使得在每個block內，每一個thread都可以訪問和操作這塊內存，而無法訪問其他block內的共享內存。這種機制就使得一個block之內的所有線程可以互相交流和合作。下面的例子中就顯示了線程之間的交流和合作。

　　這個例子計算的是兩個向量的點積。

/*
 * Copyright 1993-2010 NVIDIA Corporation.  All rights reserved.
 *
 * NVIDIA Corporation and its licensors retain all intellectual property and 
 * proprietary rights in and to this software and related documentation. 
 * Any use, reproduction, disclosure, or distribution of this software 
 * and related documentation without an express license agreement from
 * NVIDIA Corporation is strictly prohibited.
 *
 * Please refer to the applicable NVIDIA end user license agreement (EULA) 
 * associated with this source code for terms and conditions that govern 
 * your use of this NVIDIA software.
 * 
 */


#include "../common/book.h"

#define imin(a,b) (a<b?a:b)

const int N = 33 * 1024;
const int threadsPerBlock = 256;
const int blocksPerGrid =
            imin( 32, (N+threadsPerBlock-1) / threadsPerBlock );


__global__ void dot( float *a, float *b, float *c ) {
    __shared__ float cache[threadsPerBlock];
    int tid = threadIdx.x + blockIdx.x * blockDim.x;
    int cacheIndex = threadIdx.x;

    float   temp = 0;
    while (tid < N) {
        temp += a[tid] * b[tid];
        tid += blockDim.x * gridDim.x;
    }
    
    // set the cache values
    cache[cacheIndex] = temp;
    
    // synchronize threads in this block
    __syncthreads();

    // for reductions, threadsPerBlock must be a power of 2
    // because of the following code
    int i = blockDim.x/2;
    while (i != 0) {
        if (cacheIndex < i)
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();
        i /= 2;
    }

    if (cacheIndex == 0)
        c[blockIdx.x] = cache[0];
}


int main( void ) {
    float   *a, *b, c, *partial_c;
    float   *dev_a, *dev_b, *dev_partial_c;

    // allocate memory on the cpu side
    a = (float*)malloc( N*sizeof(float) );
    b = (float*)malloc( N*sizeof(float) );
    partial_c = (float*)malloc( blocksPerGrid*sizeof(float) );

    // allocate the memory on the GPU
    HANDLE_ERROR( cudaMalloc( (void**)&dev_a,
                              N*sizeof(float) ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_b,
                              N*sizeof(float) ) );
    HANDLE_ERROR( cudaMalloc( (void**)&dev_partial_c,
                              blocksPerGrid*sizeof(float) ) );

    // fill in the host memory with data
    for (int i=0; i<N; i++) {
        a[i] = i;
        b[i] = i*2;
    }

    // copy the arrays 'a' and 'b' to the GPU
    HANDLE_ERROR( cudaMemcpy( dev_a, a, N*sizeof(float),
                              cudaMemcpyHostToDevice ) );
    HANDLE_ERROR( cudaMemcpy( dev_b, b, N*sizeof(float),
                              cudaMemcpyHostToDevice ) ); 

    dot<<<blocksPerGrid,threadsPerBlock>>>( dev_a, dev_b,
                                            dev_partial_c );

    // copy the array 'c' back from the GPU to the CPU
    HANDLE_ERROR( cudaMemcpy( partial_c, dev_partial_c,
                              blocksPerGrid*sizeof(float),
                              cudaMemcpyDeviceToHost ) );

    // finish up on the CPU side
    c = 0;
    for (int i=0; i<blocksPerGrid; i++) {
        c += partial_c[i];
    }

    #define sum_squares(x)  (x*(x+1)*(2*x+1)/6)
    printf( "Does GPU value %.6g = %.6g?\n", c,
             2 * sum_squares( (float)(N - 1) ) );

    // free memory on the gpu side
    HANDLE_ERROR( cudaFree( dev_a ) );
    HANDLE_ERROR( cudaFree( dev_b ) );
    HANDLE_ERROR( cudaFree( dev_partial_c ) );

    // free memory on the cpu side
    free( a );
    free( b );
    free( partial_c );
}

　　我們首先關注核函數dot。__shared__ float cache[threadsPerBlock];就是這節重點，申請cache數組時，由於使用了共享內存，則每一個block里面都有一份cache，使得block內的thread都可以訪問和操作其各自的cache數組。

while (tid < N) {
        temp += a[tid] * b[tid];
        tid += blockDim.x * gridDim.x;
    }

這一段我們相當熟悉，每個線程計算若干對a，b的乘積，然后相加。然后這樣cache[cacheIndex] = temp;將結果存入cache中。這時，每一個線程的結果都被存在了cache數組中，我們知道接下來要對數組求和，然而這里有潛在的危險，那就是我們不知道所有線程是否已經將數據寫入了cache，也就是說，是否每一個線程都已經執行完了第39行。這里就需要等待，等待所有線程執行到同一位置，這就是 __syncthreads();的作用。這個函數稱為同步函數，即在所有線程全部執行到__syncthreads()為止，誰也不許動，其后任何代碼都無法執行。

　　因此，我們可以很清楚的明白所有線程全部執行完了第39行，然后同步解除，大家再一起往前走。做加法。

int i = blockDim.x/2;
    while (i != 0) {
        if (cacheIndex < i)
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();
        i /= 2;
    }

    if (cacheIndex == 0)
        c[blockIdx.x] = cache[0];

　　這段就不難理解了，逐對相加，最后cache【0】位置的數就是結果。將其值存入c數組，准備導出。

剩下的main函數部分是如下幾步操作（和前面學習的差不多）：

1.為輸入輸出數組分配內存

2.將a，b數組付初值，然后復制給device中，cudaMemcpy（）

3.調用核函數執行並行計算。

4.device值返回后數組c求和。

　　很明顯，由於我們使用了共享內存存儲cache數組，使得在操作cache數組時的速度有了大幅提高（相比於全局內存）。共享內存的意義也就在此。

現在，請觀察下面的兩組代碼：

 while (i != 0) {
        if (cacheIndex < i)
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();
        i /= 2;
    }

 while (i != 0) {
        if (cacheIndex < i)
        {
            cache[cacheIndex] += cache[cacheIndex + i];
        __syncthreads();
         }
        i /= 2;
    }    

下面的代碼中由於if的存在，只有部分線程包含同步操作。代碼似乎得到了優化。但是真的如此嗎

當然不是的，上面的紅字“所有線程全部執行到__syncthreads()為止”，所有很重要，<<<>>>中launch了多少個threadperblock，那么就必須要等待所有的線程，一個都不能少。由於if的存在，上例中部分線程永遠都不可能執行到cache[cacheIndex] += cache[cacheIndex + i];這一步，因此就要永遠等待下去，因而程序無法執行。

總結：在能用共享內存的時候盡量用，進而提高block內的執行效率，但是在同步問題上一定要慎重。。。