OpenCL 条件语句问题
OpenCL conditional statement issue
我试图并行化一组条件语句,但在包含内核的第一个循环执行后输出与现有实现不匹配(mapI 是一个 135 的 int 数组,并且在第 60 个索引上第二个循环它失败了,总共调用了 195 次 mapI)。我已经通过将它们与现有实现进行比较来检查所有数组是否正确地传入和传出内核,并且对为什么这个计算没有 return 正确的结果感到困惑,就像它对第一个循环所做的那样代码。所有 OpenCL 开销函数 return CL_SUCCESS。
cl_mem Q1cl, Q3cl;
Q1cl = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, sizeof(double)*um->Npts*um->Nel, Q1, &err);
Q3cl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double)*um->Npts*um->Nel, Q3, &err);
nxcl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double)*um->Nel*um->Nfaces*um->Nfq, nx, &err);
nycl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double)*um->Nel*um->Nfaces*um->Nfq, ny, &err);
mapIcl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(int)*(um->Nfq+um->Ninflow), mapI, &err);
mapFcl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(int)*(um->Nfq+um->Nfar), mapF, &err);
fluxQ1cl = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(double)*um->Nel*um->Nfaces*um->Nfq, *fluxQ1check, &err);
err = clSetKernelArg(kernel[7], 0, sizeof(cl_mem), (void*)&mapIcl);
err = clSetKernelArg(kernel[7], 1, sizeof(cl_mem), (void*)&nxcl);
err = clSetKernelArg(kernel[7], 2, sizeof(cl_mem), (void*)&nycl);
err = clSetKernelArg(kernel[7], 3, sizeof(cl_mem), (void*)&mapFcl);
err = clSetKernelArg(kernel[7], 4, sizeof(cl_mem), (void*)&Q1cl);
err = clSetKernelArg(kernel[7], 5, sizeof(cl_mem), (void*)&Q3cl);
err = clSetKernelArg(kernel[7], 6, sizeof(cl_mem), (void*)&fluxQ1cl);
globalWorkSize[0] = Ninflow; //Old implementation, now NEL
globalWorkSize[1] = Nfar; //Old implementation, now NFACES*NFQ
err = clEnqueueNDRangeKernel(queue[0], kernel[7], 2, NULL, globalWorkSize, NULL, 0, NULL, NULL);
err = clEnqueueReadBuffer(queue[0], fluxQ1cl, CL_TRUE, 0, sizeof(double)*um->Nel*um->Nfaces*um->Nfq, *fluxQ1check, 0, NULL, NULL);
内核代码:
__kernel void umBC(__global int* mapI,
__global double* nx,
__global double* ny,
__global int* mapF,
__global double* Q1,
__global double* Q3,
__global double* fluxQ1)
{
int id, idF;
double e[9][2] = { {0, 0}, {1, 0}, {0, 1}, {-1, 0}, {0, -1}, {1, 1}, {-1, 1}, {-1, -1}, {1, -1}};
double t_1 = 1. / 9.;
double uf = 0.;
double vf = 0.;
int globalx = get_global_id(0);
int globaly = get_global_id(1);
id = mapI[globalx];
fluxQ1[id] = ((e[1][0]*nx[id] + e[1][1]*ny[id]) < 0)*((Q1[id]-Q3[id] -2.*t_1*1.*(e[1][0]*uf+e[1][0]*vf)*3.) * (e[1][0]*nx[id] + e[1][1]*ny[id])) + 0.;
uf = 0.01;
vf = 0.;
idF = mapF[globaly];
fluxQ1[idF] = ((e[1][0]*nx[idF] + e[1][1]*ny[idF]) < 0)*((Q1[idF]-Q3[idF] -2.*t_1*1.*(e[1][0]*uf+e[1][0]*vf)*3.) * (e[1][0]*nx[idF] + e[1][1]*ny[idF])) + 0.;
}
编辑:下面是工作实现,再次感谢 doqtor 和 Lee 的帮助。为了实现这一点,我需要更改 mapI 和 mapF 的工作方式以匹配 fluxQ 的大小。
__kernel void umBC(__global int* mapI,
__global double* nx,
__global double* ny,
__global int* mapF,
__global double* Q1,
__global double* Q3,
__global double* fluxQ1)
{
double e[9][2] = { {0, 0}, {1, 0}, {0, 1}, {-1, 0}, {0, -1}, {1, 1}, {-1, 1}, {-1, -1}, {1, -1}};
double t_1 = 1. / 9.;
double uf = 0.;
double vf = 0.;
int globalx = get_global_id(0); //NEL
int globaly = get_global_id(1); //NFACES*NFQ
if(mapI[globalx*NFACES*NFQ+globaly] != NEL*NFACES*NFQ+1000){
fluxQ1[globalx*NFACES*NFQ+globaly] = 0.0;
if ((e[1][0]*nx[globalx*NFACES*NFQ+globaly] + e[1][1]*ny[globalx*NFACES*NFQ+globaly]) < 0){
fluxQ1[globalx*NFACES*NFQ+globaly] = (Q1[globalx*NFACES*NFQ+globaly]-Q3[globalx*NFACES*NFQ+globaly] -2.*t_1*1.*(e[1][0]*uf+e[1][0]*vf)*3.) * (e[1][0]*nx[globalx*NFACES*NFQ+globaly] + e[1][1]*ny[globalx*NFACES*NFQ+globaly]);
}
}
uf = 0.01;
vf = 0.;
if(mapF[globalx*NFACES*NFQ+globaly] != NEL*NFACES*NFQ+1000){
fluxQ1[globalx*NFACES*NFQ+globaly] = 0.0;
if ((e[1][0]*nx[globalx*NFACES*NFQ+globaly] + e[1][1]*ny[globalx*NFACES*NFQ+globaly]) < 0){
fluxQ1[globalx*NFACES*NFQ+globaly] = (Q1[globalx*NFACES*NFQ+globaly]-Q3[globalx*NFACES*NFQ+globaly] -2.*t_1*1.*(e[1][0]*uf+e[1][0]*vf)*3.) * (e[1][0]*nx[globalx*NFACES*NFQ+globaly] + e[1][1]*ny[globalx*NFACES*NFQ+globaly]);
}
}
}
您在每个维度中都有多个工作项,但您正在写入一个仅由一个索引的位置:
int globalx = get_global_id(0);
int globaly = get_global_id(1);
id = mapI[globalx];
fluxQ1[id] =
所以这里的多个工作项将共享同一个globalx。这意味着他们将从 mapI 读取相同的 id,并将写入 fluxQ1 中的相同位置。
你得到错误的结果是因为你的内核没有正确实现。 @Lee 已经说过了,我会尝试进一步解释原因。
为了简单起见,我们考虑使用 2x2 全局范围执行内核。
现在,每个工作项的 (globalx, globaly) 将是:
(0, 0)
(0, 1)
(1, 0)
(1, 1)
所以:
work item (0, 0) writes to fluxQ1[mapI[0]] and fluxQ1[mapF[0]]
work item (0, 1) writes to fluxQ1[mapI[0]] and fluxQ1[mapF[1]]
work item (1, 0) writes to fluxQ1[mapI[1]] and fluxQ1[mapF[0]]
work item (1, 1) writes to fluxQ1[mapI[1]] and fluxQ1[mapF[1]]
假设 mapI[0]、mapI[1]、mapF[0] 和 mapF[1] return 唯一 ID,那么 2 个工作项并行写入同一位置,这显然不能给出你纠正结果。如果 mapI 和 mapF 可以 return 相同的 ID,情况会变得更糟。
有些结果是靠运气得到的。
更改全局或本地工作项的数量无济于事,使用或不使用 select
也无济于事。您需要确保工作项不会同时写入同一位置!
我试图并行化一组条件语句,但在包含内核的第一个循环执行后输出与现有实现不匹配(mapI 是一个 135 的 int 数组,并且在第 60 个索引上第二个循环它失败了,总共调用了 195 次 mapI)。我已经通过将它们与现有实现进行比较来检查所有数组是否正确地传入和传出内核,并且对为什么这个计算没有 return 正确的结果感到困惑,就像它对第一个循环所做的那样代码。所有 OpenCL 开销函数 return CL_SUCCESS。
cl_mem Q1cl, Q3cl;
Q1cl = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_COPY_HOST_PTR, sizeof(double)*um->Npts*um->Nel, Q1, &err);
Q3cl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double)*um->Npts*um->Nel, Q3, &err);
nxcl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double)*um->Nel*um->Nfaces*um->Nfq, nx, &err);
nycl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_USE_HOST_PTR, sizeof(double)*um->Nel*um->Nfaces*um->Nfq, ny, &err);
mapIcl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(int)*(um->Nfq+um->Ninflow), mapI, &err);
mapFcl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(int)*(um->Nfq+um->Nfar), mapF, &err);
fluxQ1cl = clCreateBuffer(context, CL_MEM_READ_WRITE | CL_MEM_USE_HOST_PTR, sizeof(double)*um->Nel*um->Nfaces*um->Nfq, *fluxQ1check, &err);
err = clSetKernelArg(kernel[7], 0, sizeof(cl_mem), (void*)&mapIcl);
err = clSetKernelArg(kernel[7], 1, sizeof(cl_mem), (void*)&nxcl);
err = clSetKernelArg(kernel[7], 2, sizeof(cl_mem), (void*)&nycl);
err = clSetKernelArg(kernel[7], 3, sizeof(cl_mem), (void*)&mapFcl);
err = clSetKernelArg(kernel[7], 4, sizeof(cl_mem), (void*)&Q1cl);
err = clSetKernelArg(kernel[7], 5, sizeof(cl_mem), (void*)&Q3cl);
err = clSetKernelArg(kernel[7], 6, sizeof(cl_mem), (void*)&fluxQ1cl);
globalWorkSize[0] = Ninflow; //Old implementation, now NEL
globalWorkSize[1] = Nfar; //Old implementation, now NFACES*NFQ
err = clEnqueueNDRangeKernel(queue[0], kernel[7], 2, NULL, globalWorkSize, NULL, 0, NULL, NULL);
err = clEnqueueReadBuffer(queue[0], fluxQ1cl, CL_TRUE, 0, sizeof(double)*um->Nel*um->Nfaces*um->Nfq, *fluxQ1check, 0, NULL, NULL);
内核代码:
__kernel void umBC(__global int* mapI,
__global double* nx,
__global double* ny,
__global int* mapF,
__global double* Q1,
__global double* Q3,
__global double* fluxQ1)
{
int id, idF;
double e[9][2] = { {0, 0}, {1, 0}, {0, 1}, {-1, 0}, {0, -1}, {1, 1}, {-1, 1}, {-1, -1}, {1, -1}};
double t_1 = 1. / 9.;
double uf = 0.;
double vf = 0.;
int globalx = get_global_id(0);
int globaly = get_global_id(1);
id = mapI[globalx];
fluxQ1[id] = ((e[1][0]*nx[id] + e[1][1]*ny[id]) < 0)*((Q1[id]-Q3[id] -2.*t_1*1.*(e[1][0]*uf+e[1][0]*vf)*3.) * (e[1][0]*nx[id] + e[1][1]*ny[id])) + 0.;
uf = 0.01;
vf = 0.;
idF = mapF[globaly];
fluxQ1[idF] = ((e[1][0]*nx[idF] + e[1][1]*ny[idF]) < 0)*((Q1[idF]-Q3[idF] -2.*t_1*1.*(e[1][0]*uf+e[1][0]*vf)*3.) * (e[1][0]*nx[idF] + e[1][1]*ny[idF])) + 0.;
}
编辑:下面是工作实现,再次感谢 doqtor 和 Lee 的帮助。为了实现这一点,我需要更改 mapI 和 mapF 的工作方式以匹配 fluxQ 的大小。
__kernel void umBC(__global int* mapI,
__global double* nx,
__global double* ny,
__global int* mapF,
__global double* Q1,
__global double* Q3,
__global double* fluxQ1)
{
double e[9][2] = { {0, 0}, {1, 0}, {0, 1}, {-1, 0}, {0, -1}, {1, 1}, {-1, 1}, {-1, -1}, {1, -1}};
double t_1 = 1. / 9.;
double uf = 0.;
double vf = 0.;
int globalx = get_global_id(0); //NEL
int globaly = get_global_id(1); //NFACES*NFQ
if(mapI[globalx*NFACES*NFQ+globaly] != NEL*NFACES*NFQ+1000){
fluxQ1[globalx*NFACES*NFQ+globaly] = 0.0;
if ((e[1][0]*nx[globalx*NFACES*NFQ+globaly] + e[1][1]*ny[globalx*NFACES*NFQ+globaly]) < 0){
fluxQ1[globalx*NFACES*NFQ+globaly] = (Q1[globalx*NFACES*NFQ+globaly]-Q3[globalx*NFACES*NFQ+globaly] -2.*t_1*1.*(e[1][0]*uf+e[1][0]*vf)*3.) * (e[1][0]*nx[globalx*NFACES*NFQ+globaly] + e[1][1]*ny[globalx*NFACES*NFQ+globaly]);
}
}
uf = 0.01;
vf = 0.;
if(mapF[globalx*NFACES*NFQ+globaly] != NEL*NFACES*NFQ+1000){
fluxQ1[globalx*NFACES*NFQ+globaly] = 0.0;
if ((e[1][0]*nx[globalx*NFACES*NFQ+globaly] + e[1][1]*ny[globalx*NFACES*NFQ+globaly]) < 0){
fluxQ1[globalx*NFACES*NFQ+globaly] = (Q1[globalx*NFACES*NFQ+globaly]-Q3[globalx*NFACES*NFQ+globaly] -2.*t_1*1.*(e[1][0]*uf+e[1][0]*vf)*3.) * (e[1][0]*nx[globalx*NFACES*NFQ+globaly] + e[1][1]*ny[globalx*NFACES*NFQ+globaly]);
}
}
}
您在每个维度中都有多个工作项,但您正在写入一个仅由一个索引的位置:
int globalx = get_global_id(0);
int globaly = get_global_id(1);
id = mapI[globalx];
fluxQ1[id] =
所以这里的多个工作项将共享同一个globalx。这意味着他们将从 mapI 读取相同的 id,并将写入 fluxQ1 中的相同位置。
你得到错误的结果是因为你的内核没有正确实现。 @Lee 已经说过了,我会尝试进一步解释原因。
为了简单起见,我们考虑使用 2x2 全局范围执行内核。
现在,每个工作项的 (globalx, globaly) 将是:
(0, 0)
(0, 1)
(1, 0)
(1, 1)
所以:
work item (0, 0) writes to fluxQ1[mapI[0]] and fluxQ1[mapF[0]]
work item (0, 1) writes to fluxQ1[mapI[0]] and fluxQ1[mapF[1]]
work item (1, 0) writes to fluxQ1[mapI[1]] and fluxQ1[mapF[0]]
work item (1, 1) writes to fluxQ1[mapI[1]] and fluxQ1[mapF[1]]
假设 mapI[0]、mapI[1]、mapF[0] 和 mapF[1] return 唯一 ID,那么 2 个工作项并行写入同一位置,这显然不能给出你纠正结果。如果 mapI 和 mapF 可以 return 相同的 ID,情况会变得更糟。 有些结果是靠运气得到的。
更改全局或本地工作项的数量无济于事,使用或不使用 select
也无济于事。您需要确保工作项不会同时写入同一位置!