cuda(2)---方阵乘法

　　在GPU加速计算中，矩阵乘法是绕不开的话题，所以笔者就从简单的方阵乘法入手，进行cuda的入门学习

所以就写了一个简单cuda 方阵乘法的内核。

　　矩阵乘法的原理，我就不在这里赘述了，学过线性代数，都应该比较清楚。直接就上代码：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#include <stdio.h>
#include <iostream>
#include<time.h>

typedef float DataType;

#define WIDTH 32
#define Matrix_h 2
#define Matrix_w 2

#define Arr_Size(arr)  (sizeof(arr)/sizeof(arr[0][0]))
 
__global__ void MatrixMulKernel(DataType *M, DataType *N, DataType *P, int width) {
    //calcualte the row index of the P element and M
    int row = blockIdx.x*blockDim.x + threadIdx.x;
    //calcualte the col index of the P element and N
    int col = blockIdx.y*blockDim.y + threadIdx.y;

    if (row < width && col < width) {
        float pvalue = 0;
        //each thread computes one element of the block sub-matrix
        for (int k = 0; k < width; ++k) {
            pvalue += M[row*width + k] * N[k*width + col];
        }
        P[row*width + col] = pvalue;
    }
}

void MatrixMulCPU(DataType M[][Matrix_h], DataType N[][Matrix_h], DataType P[][Matrix_h], int width) {
    
    for (int i = 0; i < width; ++i) {
        for (int j = 0; j < width; ++j) {
            DataType pvalue = 0;
            for (int k = 0; k < width; ++k) {
                pvalue += M[i][k] * N[k][j];
            }
            P[i][j] = pvalue;
        }
        
    }
}


int main(void) {

    DataType A[Matrix_w][Matrix_h]={
        {1,2},
        {3,4}
    };

    DataType B[Matrix_w][Matrix_h] = {
        {5,6},
        {7,8}

    };
    DataType C[Matrix_w][Matrix_h] = { 0 };

    DataType*d_A, *d_B, *d_C;


    cudaMalloc((void**)&d_A, sizeof(DataType)*Arr_Size(A));
    cudaMalloc((void**)&d_B, sizeof(DataType)*Arr_Size(B));
    cudaMalloc((void**)&d_C, sizeof(DataType)*Arr_Size(C));
    cudaMemcpy(d_A, A, sizeof(DataType)*Arr_Size(A), cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, B, sizeof(DataType)*Arr_Size(B), cudaMemcpyHostToDevice);

    dim3 Block(WIDTH, WIDTH, 1);
    dim3 Grid((Matrix_w - 1) / WIDTH + 1, (Matrix_h - 1) / WIDTH + 1, 1);

    clock_t time_start, time_end;
    time_start = clock();
    MatrixMulKernel << <Grid, Block >> > (d_A, d_B, d_C, Matrix_w);
    cudaThreadSynchronize();
    time_end = clock();
    cudaMemcpy(C, d_C, sizeof(DataType)*Arr_Size(C), cudaMemcpyDeviceToHost);
    std::cout << "GPU time=" << (double)(time_end - time_start) / CLOCKS_PER_SEC * 1000 << "ms" << std::endl;
    std::cout << "GPU_result:" << std::endl;
    for (int i = 0; i < Matrix_w; ++i) {
        for (int j = 0; j < Matrix_h; ++j) {
            printf("C[%d][%d]=%f ", i, j, C[i][j]);
        }
        printf("
");
    }

    time_start = clock();
    MatrixMulCPU(A, B, C, Matrix_w);
    time_end = clock();

    std::cout << "CPU time=" << (double)(time_end - time_start) / CLOCKS_PER_SEC * 1000 << "ms" << std::endl;
    std::cout << "CPU_result:" << std::endl;
    for (int i = 0; i < Matrix_w; ++i) {
        for (int j = 0; j < Matrix_h; ++j) {
            printf("C[%d][%d]=%f ", i, j, C[i][j]);
        }
        printf("
");
    }

    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
    return 0;
}

　　运行结果如下图：