CUDA實作矩陣相乘

文章目錄

• 前言
• 1、簡單思路
• 分析
• 2、優化
• 總結

前言

本文主要借助CUDA實作矩陣相乘。

1、簡單思路

#include <stdio.h>

#define BLOCK_NUM  8
#define THREAD_NUM  32
#define R_SIZE BLOCK_NUM * THREAD_NUM
#define M_SIZE R_SIZE*R_SIZE

void __global__ matmul1(int *da, int *db, int *dres);

void __global__ matmul1(int *da, int *db, int *dres)
{
    // 擷取每一個線程的絕對編号，總共256條
    int tid = blockDim.x * blockIdx.x + threadIdx.x; 
    // 每一條線程計算結果矩陣一行的資料
    // 以tid = 0 為例，需要累加
    for(int c=0; c<R_SIZE; ++c)
    {
        for(int r=0; r<R_SIZE; ++r)
	    dres[tid*R_SIZE + c] += da[tid*R_SIZE+r] * db[r*R_SIZE+c];
    }
}


int main(int argc, char *argv[])
{
    //配置設定主機記憶體
    int *ha, *hb, *hres;
    ha = (int *) malloc (sizeof(int) * M_SIZE);
    hb = (int *) malloc (sizeof(int) * M_SIZE);
    hres = (int *) malloc(sizeof(int) * M_SIZE);

    //指派
    for(int i=0; i<R_SIZE; ++i)
    {
        for(int j=0; j<R_SIZE; ++j)
	{
	    ha[i*R_SIZE+j] = 1;
	    hb[i*R_SIZE+j] = 1;
	    hres[i*R_SIZE+j] = 0; 
	}
    }
    // 配置設定裝置内潤
    int *da, *db, *dres;
    cudaMalloc((void**)&da, sizeof(int)*M_SIZE);
    cudaMalloc((void**)&db, sizeof(int)*M_SIZE);
    cudaMalloc((void**)&dres, sizeof(int)*M_SIZE);

    // 拷貝資料
    cudaMemcpy(da,ha, sizeof(int)*M_SIZE, cudaMemcpyHostToDevice);
    cudaMemcpy(db,hb, sizeof(int)*M_SIZE, cudaMemcpyHostToDevice);
    cudaMemcpy(dres, hres, sizeof(int)*M_SIZE, cudaMemcpyHostToDevice);

    // 調用核函數
    matmul1<<<BLOCK_NUM,THREAD_NUM>>>(da,db,dres);

    // 拷貝資料
    cudaMemcpy(hres, dres, sizeof(int)*M_SIZE, cudaMemcpyDeviceToHost);
    
    // 列印看看
    printf("%d\n",hres[0]);

    //釋放記憶體
    free(ha);
    free(hb);
    free(hres);
    cudaFree(da);
    cudaFree(db);
    cudaFree(dres);

    return 0;
}


           

分析

 首先定義了256個線程，線程數量和矩陣的行數相等。在核函數中，變量tid擷取到了每一個線程的ID。即[0~255]。對應最終矩陣的256行。即一個線程需要計算一行的結果矩陣。假設tid =0，然後在分析核函數中的兩重循環，分别擷取da矩陣的行元素和db矩陣的列元素相乘并累加求和得到最終對應位置的解。

 後續會介紹矩陣乘法優化，根據合理的線程安排去掉一層for循環。

2、優化

#include <stdio.h>

#define BLOCK_NUM  8
#define THREAD_NUM  32
#define R_SIZE BLOCK_NUM * THREAD_NUM
#define M_SIZE R_SIZE*R_SIZE

void __global__ matmul2(int *da, int *db, int *dres);

void __global__ matmul2(int *da, int *db, int *dres)
{
    // 擷取每一個線程的ID, 編号ID:(row,col)。對應結果矩陣的 行 和 列
    int row = blockDim.y * blockIdx.y + threadIdx.y;
    int col = blockDim.x * blockIdx.x + threadIdx.x; 
    // 對應每一個的線程的結果,一個線程對應一個結果矩陣的一個元素
    for(int i=0; i<R_SIZE; ++i)
    {
        dres[row*R_SIZE + col] += da[row*R_SIZE+i] * db[i*row+col];
    }
}


int main(int argc, char *argv[])
{
    //配置設定主機記憶體
    int *ha, *hb, *hres;
    ha = (int *) malloc (sizeof(int) * M_SIZE);
    hb = (int *) malloc (sizeof(int) * M_SIZE);
    hres = (int *) malloc(sizeof(int) * M_SIZE);

    //指派
    for(int i=0; i<R_SIZE; ++i)
    {
        for(int j=0; j<R_SIZE; ++j)
	{
	    ha[i*R_SIZE+j] = 1;
	    hb[i*R_SIZE+j] = 1;
	    hres[i*R_SIZE+j] = 0; 
	}
    }
    // 配置設定裝置内潤
    int *da, *db, *dres;
    cudaMalloc((void**)&da, sizeof(int)*M_SIZE);
    cudaMalloc((void**)&db, sizeof(int)*M_SIZE);
    cudaMalloc((void**)&dres, sizeof(int)*M_SIZE);

    // 拷貝資料
    cudaMemcpy(da,ha, sizeof(int)*M_SIZE, cudaMemcpyHostToDevice);
    cudaMemcpy(db,hb, sizeof(int)*M_SIZE, cudaMemcpyHostToDevice);
    cudaMemcpy(dres, hres, sizeof(int)*M_SIZE, cudaMemcpyHostToDevice);
    
    // 調用核函數
    // 配置設定線程
    const dim3 grid_size(BLOCK_NUM, BLOCK_NUM);
    const dim3 block_size(THREAD_NUM, THREAD_NUM);

    matmul2<<<grid_size, block_size>>>(da,db,dres);

    // 拷貝資料
    cudaMemcpy(hres, dres, sizeof(int)*M_SIZE, cudaMemcpyDeviceToHost);
    
    // 列印看看
    printf("%d\n",hres[0]);

    //釋放記憶體
    free(ha);
    free(hb);
    free(hres);
    cudaFree(da);
    cudaFree(db);
    cudaFree(dres);

    return 0;
}
           

總結

多了解，線程是const dim3 block_size(8,8); 形式定義。