标签：Eigen int double Opencl 矩阵 Number 64 sizeof
经过实验发现如果计算量不够大，利用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
本站声明： 1. iCode9 技术分享网（下文简称本站）提供的所有内容，仅供技术学习、探讨和分享； 2. 关于本站的所有留言、评论、转载及引用，纯属内容发起人的个人观点，与本站观点和立场无关； 3. 关于本站的所有言论和文字，纯属内容发起人的个人观点，与本站观点和立场无关； 4. 本站文章均是网友提供，不完全保证技术分享内容的完整性、准确性、时效性、风险性和版权归属；如您发现该文章侵犯了您的权益，可联系我们第一时间进行删除； 5. 本站为非盈利性的个人网站，所有内容不会用来进行牟利，也不会利用任何形式的广告来间接获益，纯粹是为了广大技术爱好者提供技术内容和技术思想的分享性交流网站。