cuda矩阵之心得

矩阵编程，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 Cvaluefloat 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 memoryMatrix 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 memoryMatrix 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 kerneldim3 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 memorycudaMemcpy(C.elements, Cd.elements, size,cudaMemcpyDeviceToHost);// Free device memorycudaFree(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);

这只是我的一点点小小心得，如有错误欢迎您指出来，谢谢！