在GPU上,on-board memory包含以下类型:
on-chip memory包含以下类型:
HOST内存函数
DEVICE内存函数
请注意,这里函数只返回状态。所以分配的内存地址作为函数参数。
HOST《-》DEVICE互相拷贝
cudaMemcpy( 目的内存地址,源内存地址,内存大小,cudaMemcpyHostToDevice/cudaMemcpyDeviceToHost/cudaMemcpyDeviceToDevice/cudaMemcpyHostToHost)
以矩阵乘为例:
CPU的做法是嵌套循环,如上图所示。
GPU的做法应该是使用 index( blockIdx和 threadIdx的组合公式)替换原来的下标i,j。
这也是一般CUDA程序的套路——把for loop展开成每个线程处理其中的一步。
那么,如何使用CUDA将坐标拆开呢?将二维坐标(矩阵)改为 在全局中的索引:需要找到每个线程需要处理元素的位置。
ty=线程在y方向的坐标
tx=线程在x方向的坐标
ty=blockIdx.y*blockDim.y + threadIdx.y
tx=blockIdx.x*blockDim.x + threadIdx.x
nx=x方向有多少数据。
index = ty * nx + tx
目的是将高维降为低维。
矩阵乘的每个核函数的算法如下:
典型的核函数算法代码如下:
需要注意:
矩阵乘 矩阵M是 mXn,矩阵N是 nXk,这里面需要 矩阵M和矩阵N都有n。否则无法相乘。
上代码:
matrix_mul.cu
#include <stdio.h> #include <math.h> #define BLOCK_SIZE 16 //使用GPU进行矩阵计算 __global__ void gpu_matrix_mult(int *a,int *b, int *c, int m, int n, int k) { int row = blockIdx.y * blockDim.y + threadIdx.y; int col = blockIdx.x * blockDim.x + threadIdx.x; int sum = 0; if( col < k && row < m) { for(int i = 0; i < n; i++) { sum += a[row * n + i] * b[i * k + col]; } c[row * k + col] = sum; } } //使用CPU进行矩阵计算 void cpu_matrix_mult(int *h_a, int *h_b, int *h_result, int m, int n, int k) { for (int i = 0; i < m; ++i) { for (int j = 0; j < k; ++j) { int tmp = 0.0; for (int h = 0; h < n; ++h) { tmp += h_a[i * n + h] * h_b[h * k + j]; } h_result[i * k + j] = tmp; } } } int main(int argc, char const *argv[]) { /* 矩阵A mXn,矩阵B nXk --》矩阵乘计算的结果是 mXk */ int m=3; int n=4; int k=5; int *h_a, *h_b, *h_c, *h_cc; //分配原矩阵的内存 h是host memory cudaMallocHost((void **) &h_a, sizeof(int)*m*n); cudaMallocHost((void **) &h_b, sizeof(int)*n*k); //分配 CPU结果内存 cudaMallocHost((void **) &h_c, sizeof(int)*m*k); //分配 GPU结果内存 cudaMallocHost((void **) &h_cc, sizeof(int)*m*k); //初始化矩阵A(mxn) srand(time(0)); printf("---------------h_a------------------\n"); for (int i = 0; i < m; ++i) { for (int j = 0; j < n; ++j) { h_a[i * n + j] = rand() % 1024; printf("%d", h_a[i * n + j] ); printf(" "); } printf("\n"); } //初始化矩阵B(nxk) printf("---------------h_b------------------\n"); for (int i = 0; i < n; ++i) { for (int j = 0; j < k; ++j) { h_b[i * k + j] = rand() % 1024; printf("%d", h_b[i * k + j] ); printf(" "); } printf("\n"); } int *d_a, *d_b, *d_c; //分配 原矩阵的GPU内存 d是device memory cudaMalloc((void **) &d_a, sizeof(int)*m*n); cudaMalloc((void **) &d_b, sizeof(int)*n*k); //分配 目的矩阵的GPU内存 cudaMalloc((void **) &d_c, sizeof(int)*m*k); // copy matrix A and B from host to device memory cudaMemcpy(d_a, h_a, sizeof(int)*m*n, cudaMemcpyHostToDevice); cudaMemcpy(d_b, h_b, sizeof(int)*n*k, cudaMemcpyHostToDevice); unsigned int grid_rows = (m + BLOCK_SIZE - 1) / BLOCK_SIZE; unsigned int grid_cols = (k + BLOCK_SIZE - 1) / BLOCK_SIZE; dim3 dimGrid(grid_cols, grid_rows); dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE); //GPU计算,结果放入h_c gpu_matrix_mult<<<dimGrid, dimBlock>>>(d_a, d_b, d_c, m, n, k); cudaMemcpy(h_c, d_c, sizeof(int)*m*k, cudaMemcpyDeviceToHost); //cudaThreadSynchronize(); //CPU计算,结果直接放入h_cc cpu_matrix_mult(h_a, h_b, h_cc, m, n, k); int ok = 1; for (int i = 0; i < m; ++i) { for (int j = 0; j < k; ++j) { // 比较大小的时候使用 a-b<0.0000000001 if(fabs(h_cc[i*k + j] - h_c[i*k + j])>(1.0e-10)) { ok = 0; } } } printf("---------------h_c cpu result------------------\n"); for(int i=0;i<m;i++) { for(int j=0;j<k;j++) { //矩阵小的时候还可以打印,大的时候就别打了 printf("%d",h_c[i*k + j] ); printf(" "); } printf("\n"); } printf("---------------h_cc gpu result----------------\n"); for(int i=0;i<m;i++) { for(int j=0;j<k;j++) { //矩阵小的时候还可以打印,大的时候就别打了 printf("%d",h_cc[i*k + j] ); printf(" "); } printf("\n"); } if(ok) { printf("Pass!!!\n"); } else { printf("Error!!!\n"); } // free memory cudaFree(d_a); cudaFree(d_b); cudaFree(d_c); cudaFreeHost(h_a); cudaFreeHost(h_b); cudaFreeHost(h_c); return 0; }
代码中张小白加上了注释,已经介绍得比较清楚了。
我们执行下看看:
代码以 3X4和4X5的矩阵相乘,得到了3X5的矩阵结果。
这个结果跟CPU计算的结果做了对比。显示Pass表示结果是一致的(其实张小白把两个结果都打印的出来,当然也是一致的)
这里面有个小TIPS,就是在调用rand()生成随机数的时候,可以使用srand(time(0)) 做随机数种子,这样下次调用的时候跟这次生成的内容就会不一样。如果去掉这句话,每次执行的结果都是一样的。
当然,如果在同一秒同时执行,srand(time(0)) 也会导致同时生成的随机数是一样的。这点需要注意。