利用Opencl加速Eigen矩阵(二)SHSHSHSH1的博客-

经过实验发现如果计算量不够大，利用opencl反而浪费时间。所以本实验进行加速的原来代码如下：

Matrix<double,8,8>A[64],B[64],D[64]; //…… //A,B,D初始化之后进行计算 for(int i=0;i<64;i++) {  D[i].noalias() +=A * B * A.transpose(); } 

接下来进行opencl进行加速，对64个8*8的矩阵计算并行计算。
第一步：用Eigen声明矩阵，并赋初值。这里就是简单地对矩阵乘法重复64次，所以初始化一组矩阵就行了。

Matrix<double, 8, 8> eigMatA,eigMatB;//两个计算矩阵 Matrix<double, 8, 8> eigMatC1;//利用opencl计算出来的矩阵 Matrix<double, 8, 8> eigMatD1[64];//利用Eigen计算出来的矩阵  for (int i = 0; i < 8; i++)  {   for (int j = 0; j < 8; j++)   {    eigMatA(i, j) = i * j+1;    eigMatB(i, j) = i * j+2;     eigMatC(i, j) = 1;     eigMatD(i, j) = 1;   }  }  

第二步：声明指针

//声明指针 //Number=64  double* pA = (double*)malloc(sizeof(double) * eigMatA.size() * Number);//对应A  double* pB = (double*)malloc(sizeof(double) * eigMatB.size() * Number);//对应B  double* pC1 = (double*)malloc(sizeof(double) * eigMatC1.size() * Number);//对应C 

第三步：利用Map进行映射（这里相关原理在（一）中已经说明）

//Map进行映射，这里默认64个相同矩阵，实际中只需通过指针移动指向相应矩阵即可 for (int i = 0; i < Number; i++)  {   Map<Matrix<double, 8, 8>, RowMajor>(pA + i * eigMatA.rows()* eigMatA.cols(), eigMatA.rows(), eigMatA.cols()) = eigMatA;         Map<Matrix<double, 8, 8>, RowMajor>(pB + i * eigMatB.rows()* eigMatB.cols(), eigMatB.rows(), eigMatB.cols()) = eigMatB;         Map<Matrix<double, 8, 8>, RowMajor>(pC1 + i * eigMatC1.rows()* eigMatC1.cols(), eigMatC1.rows(), eigMatC1.cols()) = eigMatC1;  } 

第四步：创建内存对象，这一步就是OpenCL的知识了，这里假设大家已经了解OpenCL的总流程了，本文重点讲解OpenCL加速矩阵计算。

// 创建输入内存对象 memInutBuffer1 = clCreateBuffer(   Context,   CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,  // 输入内存为只读，并可以从宿主机内存复制到设备内存   Number * sizeof(double) * eigMatA.size(),    // 输入内存空间大小   (void*)pA,   NULL);   memInutBuffer2 = clCreateBuffer(   Context,   CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR,  // 输入内存为只读，并可以从宿主机内存复制到设备内存   Number * sizeof(double) * eigMatB.size(),    // 输入内存空间大小   (void*)pB,   NULL);   // 创建输出内存对象  memOutputBuffer1 = clCreateBuffer(   Context,   CL_MEM_READ_WRITE| CL_MEM_COPY_HOST_PTR,     // 输出内存能读写   Number * sizeof(double) * eigMatC1.size(), // 输出内存空间大小   (void*)pC1,   NULL); 

第五步：创建内核函数

//创建核函数  cl_kernel  Kernel = NULL;    // 内核对象  Kernel = clCreateKernel(Program, "add", NULL);// cl文件中的入口函数  if (NULL == Kernel)  {   cout << "Error: Can not create kernel" << endl;   return 0;  } 

介绍内核函数add

__kernel void add(const int Ndim,const int Mdim,const int Pdim, const int Nframes,  __global const double* A, __global const double* B, __global double* C1) { //这里约定：矩阵A大小是Ndim*Mdim,矩阵B大小是Mdim*Pdim，矩阵C1大小是Ndim*Pdim，一个矩阵的大小为Nframes //这里理论上输入时64*8*8的三维输入，但是这里选择输入为64*8的两维输入，每次计算一行  int m = get_global_id(0);//获取第一维度位置  int i = get_global_id(1);//获取第二维度位置  int mj, mk;  double tmp;  double ctmp;  double tmpAline[8];//将A矩阵的某一行先读到设备的缓存中，减少数据搬移时间  double tmpCline[8];//暂存A*B的的一行数据  //计算C1 if (i < Ndim)//判读是否超界  {  //这里要特别注意一个问题  //切记opencl的设备读入数组均是按照列排序的,大家可以实验如果按照按行读入理解计算，算出来的是利用Eigen计算出来的转置  /*  假设一个矩阵Matrix为N*M，矩阵大小为numberNM,共有Number个这样的矩阵  由于指针是按列排序读入 所以假设当前是第m个矩阵  矩阵的第i行：  for(int j=0;j<M;j++)    Matrix[m*numberNM+j*N+i]    矩阵的第j列：   for(int i=0;i<N;i++)    Matrix[m*numberNM+i+N*j];  */   //首先先读取A的第i行数据   for (mk = 0; mk < Pdim; mk++)    {     tmpAline[mk] = A[m * Nframes + i + mk* Ndim];   } //然后得到A*B的第i行数据，用A的第i行分别乘以B的每一列  for (mj = 0; mj < Mdim; mj++)     {      tmp = 0.0;      for (mk = 0; mk < Pdim; mk++)         {        tmp += tmpAline[mk] * B[m * Nframes + mk + mj* Pdim];        }      tmpCline[mj] = tmp;     }    //再计算C1第i行数据，因为这里乘以A的转置，实际上就是乘以A的每一行    for (mj = 0; mj < Ndim; mj++)     {      tmp = 0.0;      for (mk = 0; mk < Pdim; mk++)      {        tmp += tmpCline[mk] * A[m * Nframes + mj + mk* Pdim];      }     tmp = tmp + C1[m * Nframes + i + mj*Ndim];      C1[m * Nframes + i  + mj*Ndim] = tmp;     }   } } 

第六步：设置内核参数

//设置内核参数 //Ndim=8,Mdim=8,Pdim=8,NMP=8*8=64  iStatus = clSetKernelArg(Kernel,  0, sizeof(int), &Ndim);  iStatus |= clSetKernelArg(Kernel, 1, sizeof(int), &Mdim);  iStatus |= clSetKernelArg(Kernel, 2, sizeof(int), &Pdim);  iStatus |= clSetKernelArg(Kernel, 3, sizeof(int), &NMP);  iStatus |= clSetKernelArg(Kernel, 4, sizeof(cl_mem), (void*)&memInutBuffer1);  iStatus |= clSetKernelArg(Kernel, 5, sizeof(cl_mem), (void*)&memInutBuffer2);  iStatus |= clSetKernelArg(Kernel, 6, sizeof(cl_mem), (void*)&memOutputBuffer1); 

第七步：运行内核

//运行内核  size_t   uiGlobal_Work_Size[2] = { 0 };  // 用于设定内核分布  uiGlobal_Work_Size[0] = Number;//64，代表64个矩阵  uiGlobal_Work_Size[1] = Ndim;  //8，代表每个矩阵的行数，计算的时候是一行一行计算的 // 利用命令队列使将再设备上执行的内核排队  iStatus = clEnqueueNDRangeKernel(   CommandQueue,   Kernel,   2,   NULL,   uiGlobal_Work_Size,  // 确定内核在设备上的多个处理单元间的分布   NULL,     // 确定内核在设备上的多个处理单元间的分布   0,   NULL,   &event);  

第八步：读取内核计算结果

//读取结果    iStatus = clEnqueueReadBuffer(   CommandQueue,  // 命令队列   memOutputBuffer1, // 输出内存对象   CL_TRUE,   // 内核读取结束之前该函数不会返回   0,   sizeof(double) * eigMatC1.size() * Number,   pC1,   0,   NULL,   NULL); 

第九步：将结果映射回矩阵，注意由于这里矩阵是按列读取的，所以在内核计算的时候就考虑了，这里这样读出后，就是正常的矩阵了。

for (int i = 0; i < Number; i++)  {   eigMatC1 = Map<Matrix<double, 8, 8>, RowMajor>(pC1 + i * 64, eigMatC1.rows(), eigMatC1.cols());   } 

第十步：利用Eigen直接计算矩阵

for (int i = 0; i < Number; i++)  {  eigMatD1[i].noalias() += eigMatA * eigMatB * eigMatA.transpose(); } 

最后大家可以比较，矩阵eigMatC1 与矩阵eigMatD1一样，同时可以比较两者时间。
比较时间时虽然Opencl提供了计算内核时间的函数，但是个人感觉没有意义，因为在实际应用中，需要考虑到数据搬移的时间。所以建议opencl的运行时间约定为从创建设备内存对象开始一直到将读取数据copy到矩阵为止。这时计算时间就可以通过获取时间戳相减就行了，建议使用chrono库，精度高一些。