利用Opencl加速Eigen矩阵(二)

2024-02-05 18:27:52

经过实验发现如果计算量不够大，利用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[i] * B[i] * A[i].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库，精度高一些。

码农公寓

相关文章