用Cuda计算两个矩阵相加
好了,现在我们来学习一个稍微复杂一点的程序,计算两个矩阵的和。Cuda被设计出来就是用来做并行计算的,如果只是像第一个例子那样计算两个数字的和,岂不是太浪费了?
//matrix_add.cu
//这个文件不推荐直接拷贝,最好一行一行地抄下来
#include <stdio.h>
#define N 32 //这里定义了矩阵的阶级,这里用一个32x32的方形矩阵做例子
//想要自己Debug的时候,遇到问题,可以先把N的值设置小一些,比如4
/*
* 本次的核函数,三个参数分别是两个NxN的输入矩阵和一个NxN的输出矩阵
*/
__global__ void matrix_add(const int a[][N], const int b[][N], int c[][N]) {
int idx = threadIdx.x;
int idy = threadIdx.y;
c[idx][idy] = a[idx][idy] + b[idx][idy];
}
int main(void) {
int *h_a, *h_b, *h_c;
int *dev_a, *dev_b, *dev_c;
//这里把一个block的里面的线程排列定义成二维的,这样会有两个维度的索引值,x和y
dim3 threads_in_block (N, N);
//err这个值是用来检查cuda的函数是否正常运行的
cudaError_t err = cudaSuccess;
h_a = (int *)malloc(sizeof(int) * N * N);
h_b = (int *)malloc(sizeof(int) * N * N);
h_c = (int *)malloc(sizeof(int) * N * N);
if (h_a == NULL || h_b == NULL || h_c == NULL) {
fprintf(stderr, "Malloc() failed.\n");
return -1;
}
err = cudaMalloc((void **)&dev_a, sizeof(int) * N * N);
if (err != cudaSuccess) {
fprintf(stderr, "cudaMalloc() failed.\n");
return -1;
}
err = cudaMalloc((void **)&dev_b, sizeof(int) * N * N);
if (err != cudaSuccess) {
fprintf(stderr, "cudaMalloc() failed.\n");
return -1;
}
err = cudaMalloc((void **)&dev_c, sizeof(int) * N * N);
if (err != cudaSuccess) {
fprintf(stderr, "cudaMalloc() failed.\n");
return -1;
}
for (int i = 0; i < N * N; i++) {
h_a[i] = 2 * i + 1;
h_b[i] = -1 * i + 5;
}
err = cudaMemcpy(dev_a, h_a, sizeof(int) * N * N, cudaMemcpyHostToDevice);
if (err != cudaSuccess) {
fprintf(stderr, "cudaMemcpy() failed.\n");
return -1;
}
err = cudaMemcpy(dev_b, h_b, sizeof(int) * N * N, cudaMemcpyHostToDevice);
if (err != cudaSuccess) {
fprintf(stderr, "cudaMemcpy() failed.\n");
return -1;
}
matrix_add<<<1, threads_in_block>>>((int (*)[N])dev_a, (int (*)[N])dev_b, (int (*)[N])dev_c);
err = cudaMemcpy(h_c, dev_c, sizeof(int) * N * N, cudaMemcpyDeviceToHost);
if (err != cudaSuccess) {
fprintf(stderr, "cudaMemcpy() failed.\n");
return -1;
}
for (int i = 0; i < N * N; i++) {
if (h_a[i] + h_b[i] != h_c[i]) {
fprintf(stderr, "a[%d]%d + b[%d]%d != c[%d]%d.\n", i, h_a[i], i, h_b[i], i, h_c[i]);
return -1;
}
}
printf("done.\n");
return 0;
}
上面的代码对比第一份最简单的代码,其实只有两个改变。一是多了一个错误判断,方便更好的定位问题。二是,这里计算的是两个NxN的矩阵相加。这里把一个Block定义成了NxN的大小,也就是说一个Block里面有NxN个线程(例子是32x32即是1024个线程),这些线程同时计算一个相加的步骤,但是它们相互之间没有任何联系。
这个例子把在CPU中要循环1024次的加法计算,分散到1024个cuda线程来执行。你可以试一下把N的值改得大一些,就会发现,程序出错了。这是因为每个block上面可以包含的线程数目不可能是无限的。但是如果我需要计算两个很大的矩阵的和那怎么办?那就先把程序改为多个block呗。
__global__ void matrix_add(const int a[][N], const int b[][N], int c[][N]) {
}