矩阵加法

▶ 直接的矩阵加法，没有优化

#include <stdio.h>
#include <stdlib.h>
#include <malloc.h>
#include <time.h>
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

#define WIDTH       (1024*1+19)
#define HEIGHT      (1024*1+17)
#define THREAD_SIZE 16
#define BLOCK_SIZE  256
#define DIMENSION   1   //使用指定维数的跳转
#define SEED        1   //(unsigned int)time(MULL)

void checkNULL(void *input)
{
    if (input == NULL)
    {
        printf("
	find a NULL!");
        exit(1);
    }
    return;
}

void checkCudaError(cudaError input)
{
    if (input != cudaSuccess)
    {
        printf("
	find a cudaError!");
        exit(1);
    }
    return;
}

int checkResult(float * in1, float * in2, const int length)
// 注意返回值为0（两向量相等）或者“值不等的元素下标加一”（防止0号元素就不相等），返回后若想使用该下标则需要减1
{
    int i;
    for (i = 0; i < length; i++)
    {
        if (in1[i] != in2[i])
            return i + 1;
    }
    return 0;
}

// 采用一维跳转
__global__ void add1(const float *a, const float *b, float *c)
{
    int id = blockIdx.x * blockDim.x + threadIdx.x;

    while (id < WIDTH * HEIGHT)
    {
        c[id] = a[id] + b[id];
        id += gridDim.x*blockDim.x;
    }
    return;
}

// 采用二维跳转
__global__ void add2(const float *a, const float *b, float *c)
{
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    int idy = blockIdx.y * blockDim.y + threadIdx.y;
    int id = idx + idy *blockDim.y *gridDim.y;

    while (id < WIDTH * HEIGHT)
    {
        c[id] = a[id] + b[id];
        id += gridDim.x * blockDim.x * gridDim.y * blockDim.y; 
    }
    return;
}

int main()
{
    int i;
    float *h_a, *h_b, *cpu_out, *gpu_out;
    float *d_a, *d_b, *d_c;
    clock_t timeCPU;
    cudaEvent_t start, stop;
    float timeGPU;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    printf("
	Start! Matrix dimension:	%d × %d", WIDTH, HEIGHT);

    checkNULL(h_a = (float *)malloc(sizeof(float)*WIDTH*HEIGHT));
    checkNULL(h_b = (float *)malloc(sizeof(float)*WIDTH*HEIGHT));
    checkNULL(cpu_out = (float *)malloc(sizeof(float)*WIDTH*HEIGHT));
    checkNULL(gpu_out = (float *)malloc(sizeof(float)*WIDTH*HEIGHT));
    checkCudaError(cudaMalloc((float **)&d_a, sizeof(float)*WIDTH*HEIGHT));
    checkCudaError(cudaMalloc((float **)&d_b, sizeof(float)*WIDTH*HEIGHT));
    checkCudaError(cudaMalloc((float **)&d_c, sizeof(float)*WIDTH*HEIGHT));

    // 初始化
    timeCPU = clock();
    srand(SEED);
    for (i = 0; i < WIDTH*HEIGHT; i++)
    {
        h_a[i] = (float)rand() / RAND_MAX;
        h_b[i] = (float)rand() / RAND_MAX;
    }
    printf("
	Initialized! Time:		%8.3f ms
", (float)(clock() - timeCPU));

    timeCPU = clock();
    for (i = 0; i < WIDTH*HEIGHT; i++)
        cpu_out[i] = h_a[i] + h_b[i];
    printf("
	CPU time:	%8.3f ms
", (float)(clock() - timeCPU)); 

    cudaMemset(d_c, 0, sizeof(float) * WIDTH * HEIGHT);
    cudaEventRecord(start, 0);
    cudaMemcpy(d_a, h_a, sizeof(float) * WIDTH * HEIGHT, cudaMemcpyHostToDevice);
    cudaMemcpy(d_b, h_b, sizeof(float) * WIDTH * HEIGHT, cudaMemcpyHostToDevice);
    if(DIMENSION == 1)
        add1 << < BLOCK_SIZE * BLOCK_SIZE, THREAD_SIZE * THREAD_SIZE >> > (d_a, d_b, d_c);
    else
        add2 << < dim3(BLOCK_SIZE, BLOCK_SIZE), dim3(THREAD_SIZE, THREAD_SIZE) >> > (d_a, d_b, d_c);
    cudaMemcpy(gpu_out, d_c, sizeof(float) * WIDTH * HEIGHT, cudaMemcpyDeviceToHost);
    cudaDeviceSynchronize();
    cudaEventRecord(stop,0);
    cudaEventSynchronize(stop);
    cudaEventElapsedTime(&timeGPU, start, stop);
    if (i = checkResult(cpu_out, gpu_out, WIDTH * HEIGHT))
        printf("
	Compute error at i = %d
	cpu_out[i] = %10d, gpu_out[i] = %10d
", i - 1, cpu_out[i - 1], gpu_out[i - 1]);
    else
        printf("
	GPU Compute correctly!
");
    printf("
	GPU Time:	%8.3f ms
", timeGPU);

    free(h_a);
    free(h_b);
    free(cpu_out);
    free(gpu_out);
    cudaFree(d_a);
    cudaFree(d_b);
    cudaFree(d_c);
    getchar();
    return 0;
}

▶输出结果，对于超大型矩阵的加法计算，左下图为一维跳转，右下图为二维跳转。矩阵较大时初始化矩阵费时较多。超大矩阵的计算上GPU相对CPU产生了一定的优势，一维跳转比二维跳转稍快，因为省去了复杂的下标映射。

▶ 对于较小的矩阵，GPU没有发挥出优势， 因为内存拷贝等方面耗时明显。