利用Opencl加速Eigen矩阵(二)
作者:互联网
经过实验发现如果计算量不够大,利用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库,精度高一些。
标签:Eigen,int,double,Opencl,矩阵,Number,64,sizeof 来源: https://blog.csdn.net/SHSHSHSH1/article/details/105552640