矩阵编程,cuda编程入门
一、先了解什么是矩阵
在数学中,矩阵(Matrix)是一个按照长方阵列排列的复数或实数集合[1],最早来自于方程组的系数及常数所构成的方阵。
二、矩阵的基本参数
typedef struct
{
int width;
int height;
float* elements;
} Matrix;
三、在CPU上使用c/c++进行矩阵编程
- 先来个简单点的方阵举例
for(i = 0; i < n; i++) { for(j = 0; j < n; j++) { C[i][j] = 0; //将和矩阵初始化 for(k = 0; k < n; k++) { C[i][j] += A[i][k] * B[k][j]; } } }
- 然后一般矩阵举例
void MatrixMulCPU(float* C,const float *A,const float *B,int wa,int ha,int wb)//wa表示a的width,ha表示a的height,其实这里wb表示的是B矩阵的高 { float sum = 0; for (int i = 0; i < ha; ++i) { for (int j = 0; j < wb; ++j) { sum = 0; for (int k = 0; k < wa; ++k) { sum += (float)_A[i*wa+k]*(float)B[k*wb+ j]; } C[i*wb+j] = (float)sum; } } }
- 综合一下,cpu进行矩阵计算需使用三次for循环,时间复杂度O(n^3);
四、GPU使用cuda编程
简单例子举例:
请仔细观察图中的参数,关注红线的指引!
这里简单介绍一下cuda的编程结构
- cuda编程结构分为两部分,第一部分Device、第二部分为Host
- Device中最重要的是两部分,1、线程的分配(进行线程编号),2、运算指令描述
- Host大致分为三部分,1、变量申明、内存分配,2、数据传输、调用kernel函数,3释放内存空间
Device部分:
先看代码逐步解释:
__global__ void MatMulKernel(Matrix A, Matrix B, Matrix C)
{
// Each thread computes one element of C
// by accumulating results into Cvalue
float Cvalue = 0;
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
for (int e = 0; e < A.width; ++e)
Cvalue += A.elements[row * A.width + e]
* B.elements[e * B.width + col];
C.elements[row * C.width + col] = Cvalue;
}
这里的row与col让人难受,与我们习以为常的x与y轴有点区别,让人头大不要担心我会带你们慢慢理顺
我们习以为常的row主列的矩阵,
int row = blockIdx.y * blockDim.y + threadIdx.y;
刚开始我特别不理解为什么不是.x而是.y,反向思考如果确定了y,只有x可变时不就是代表着一行了嘛!(一下子转不过来不要紧,多思考几次就习惯了)
经过上述cpu矩阵运算示例,我们大致知道需要三次for循环才能完成矩阵运算!你们有没有注意到这里只使用了一次for循环就搞定了!是不是感觉哇塞,好酷!
接下来我们来解读一下这一个for循环:
for (int e = 0; e < A.width; ++e)
Cvalue += A.elements[row * A.width + e] * B.elements[e * B.width + col];
相当简单的一个for循环语句,一共两行代码就描述了
运算指令(与之前的cuda编程结构相对应)
row*A.width+e这代表着什么,开始一直困惑着我和小伙伴讨论过后弄明白了,如果你仔细看了我刚刚提到的x与y的关系时,你会醍醐灌顶,恍然大悟!确定了x值,变化着y值,表示这A矩阵的一列,e*B.width+col与之类似,自行思考;
Host部分:
void MatMul(const Matrix A, const Matrix B, Matrix C)
{
// Load A and B to device memory
Matrix d_A;
d_A.width = A.width; d_A.height = A.height;
size_t size = A.width * A.height * sizeof(float);
cudaMalloc(&d_A.elements, size);
cudaMemcpy(d_A.elements, A.elements, size,
cudaMemcpyHostToDevice);
Matrix d_B;
d_B.width = B.width; d_B.height = B.height;
size = B.width * B.height * sizeof(float);
cudaMalloc(&d_B.elements, size);
cudaMemcpy(d_B.elements, B.elements, size,
cudaMemcpyHostToDevice);
// Allocate C in device memory
Matrix d_C;
d_C.width = C.width; d_C.height = C.height;
size = C.width * C.height * sizeof(float);
cudaMalloc(&d_C.elements, size);
// Invoke kernel
dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y);
MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C);
// Read C from device memory
cudaMemcpy(C.elements, Cd.elements, size,
cudaMemcpyDeviceToHost);
// Free device memory
cudaFree(d_A.elements);
cudaFree(d_B.elements);
cudaFree(d_C.elements);
}
看见这么长的一段代码,莫怕(我带你慢慢分解):
第一步:申明变量,分配空间,赋值
Matrix d_A;
d_A.width = A.width;
d_A.height = A.height;
size_t size = A.width * A.height * sizeof(float);
cudaMalloc(&d_A.elements, size);
cudaMemcpy(d_A.elements, A.elements, size,cudaMemcpyHostToDevice);
第二步:调用kernel函数:
dim3 dimBlock(BLOCK_SIZE, BLOCK_SIZE);
dim3 dimGrid(B.width / dimBlock.x, A.height / dimBlock.y);
MatMulKernel<<<dimGrid, dimBlock>>>(d_A, d_B, d_C);
第三步:释放空间:
cudaFree(d_A.elements);
cudaFree(d_B.elements);
cudaFree(d_C.elements);
这只是我的一点点小小心得,如有错误欢迎您指出来,谢谢!