GPU程式設計(五):利用好sharedmemory

目錄

• 前言
• CPU矩陣轉置
• GPU實現

• 簡單移植
• 單block
• tile

• 利用率計算
• shared memory
• 最後

---

前言

之前在第三章對比過CPU和GPU, 差距非常大. 這一次來看看GPU自身的優化, 主要是shared memory的用法.

---

CPU矩陣轉置

矩陣轉置不是什麼複雜的事情. 用CPU實現是很簡單的:

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>

#define LOG_
#define N 1024

/* 轉置 */
void transposeCPU( float in[], float out[] )
{
    for ( int j = 0; j < N; j++ )
    {
        for ( int i = 0; i < N; i++ )
        {
            out[j * N + i] = in[i * N + j];
        }
    }
}


/* 列印矩陣 */
void logM( float m[] )
{
    for ( int i = 0; i < N; i++ )
    {
        for ( int j = 0; j < N; j++ )
        {
            printf( "%.1f ", m[i * N + j] );
        }
        printf( "
" );
    }
}


int main()
{
    int    size    = N * N * sizeof(float);
    float    *in    = (float *) malloc( size );
    float    *out    = (float *) malloc( size );

    /* 矩陣賦值 */
    for ( int i = 0; i < N; ++i )
    {
        for ( int j = 0; j < N; ++j )
        {
            in[i * N + j] = i * N + j;
        }
    }

    struct timeval    start, end;
    double        timeuse;
    int        sum = 0;
    gettimeofday( &start, NULL );

    transposeCPU( in, out );

    gettimeofday( &end, NULL );
    timeuse = end.tv_sec - start.tv_sec + (end.tv_usec - start.tv_usec) / 1000000.0;
    printf( "Use Time: %fs
", timeuse );

#ifdef LOG
    logM( in );
    printf( "
" );
    logM( out );
#endif

    free( in );
    free( out );
    return(0);
}

---

GPU實現

簡單移植

如果什麼都不考慮, 只是把程式碼移植到GPU:

#include <stdio.h>
#include <stdlib.h>
#include <sys/time.h>

#define N 1024
#define LOG_

/* 轉置 */
__global__ void transposeSerial( float in[], float out[] )
{
    for ( int j = 0; j < N; j++ )
        for ( int i = 0; i < N; i++ )
            out[j * N + i] = in[i * N + j];
}

/* 列印矩陣 */
void logM( float m[] ){...}

int main()
{
    int size = N * N * sizeof(float);

    float *in, *out;

    cudaMallocManaged( &in, size );
    cudaMallocManaged( &out, size );

    for ( int i = 0; i < N; ++i )
        for ( int j = 0; j < N; ++j )
            in[i * N + j] = i * N + j;

    struct timeval    start, end;
    double        timeuse;
    gettimeofday( &start, NULL );

    transposeSerial << < 1, 1 >> > (in, out);

    cudaDeviceSynchronize();

    gettimeofday( &end, NULL );
    timeuse = end.tv_sec - start.tv_sec + (end.tv_usec - start.tv_usec) / 1000000.0;
    printf( "Use Time: %fs
", timeuse );


#ifdef LOG
    logM( in );
    printf( "
" );
    logM( out );
#endif

    cudaFree( in );
    cudaFree( out );
}

不用想, 這裡肯定是還不如單執行緒的CPU的, 真的是完完全全的資源浪費. 實測下來, 耗時是CPU的20多倍, 大寫的丟人.

單block

單block最多可以開1024執行緒, 這裡就開1024執行緒跑下.

/* 轉置 */
__global__ void transposeParallelPerRow( float in[], float out[] )
{
    int i = threadIdx.x;
    for ( int j = 0; j < N; j++ )
        out[j * N + i] = in[i * N + j];
}

int main()
{
    ...
    transposeParallelPerRow << < 1, N >> > (in, out);
    ...
}

效率一下就提升了, 耗時大幅下降.

tile

但是的話, 如果可以利用多個block, 把矩陣切成更多的tile, 效率還會進一步提升.

/* 轉置 */
__global__ void transposeParallelPerElement( float in[], float out[] )
{
    int i = blockIdx.x * K + threadIdx.x;
    /* column */
    int j = blockIdx.y * K + threadIdx.y;
    /* row */
    out[j * N + i] = in[i * N + j];
}

int main()
{
    ...
    dim3 blocks( N / K, N / K );
    dim3 threads( K, K );

    ...
    
    transposeParallelPerElement << < blocks, threads >> > (in, out);
    ...
}

這些都是GPU的常規操作, 但其實利用率依舊是有限的.

---

利用率計算

利用率是可以粗略計算的, 比方說, 這裡的Memory Clock rate和Memory Bus Width是900Mhz和128-bit, 所以峰值就是14.4GB/s.

之前的最短耗時是0.001681s. 資料量是1024*1024*4(Byte)*2(讀寫). 所以是4.65GB/s. 利用率就是32%. 如果40%算及格, 這個利用率還是不及格的.

---

shared memory

那該如何提升呢? 問題在於讀資料的時候是連著讀的, 一個warp讀32個資料, 可以同步操作, 但是寫的時候就是散開來寫的, 有一個很大的步長. 這就導致了效率下降. 所以需要藉助shared memory, 由他轉置資料, 這樣, 寫入的時候也是連續高效的了.

/* 轉置 */
__global__ void transposeParallelPerElementTiled( float in[], float out[] )
{
    int    in_corner_i    = blockIdx.x * K, in_corner_j = blockIdx.y * K;
    int    out_corner_i    = blockIdx.y * K, out_corner_j = blockIdx.x * K;

    int x = threadIdx.x, y = threadIdx.y;

    __shared__ float tile[K][K];

    tile[y][x] = in[(in_corner_i + x) + (in_corner_j + y) * N];
    __syncthreads();
    out[(out_corner_i + x) + (out_corner_j + y) * N] = tile[x][y];
}

int main()
{

    ...
    dim3 blocks( N / K, N / K );
    dim3 threads( K, K );

    struct timeval    start, end;
    double        timeuse;
    gettimeofday( &start, NULL );

    transposeParallelPerElementTiled << < blocks, threads >> > (in, out);
    ...

}

這樣利用率就來到了44%, 及格了.

所以這就是依據架構來設計演算法, 回顧一下架構圖:

---

最後

但是44%也就是達到了及格線, 也就是說, 還有更深層次的優化工作需要做. 這些內容也就放在後續文章中了, 有意見或者建議評論區見~

---