pyopencl倍增问题

Question

这似乎是一个奇怪但非常基本的问题。我试着用pyopencl做了一个简单的操作。下面给出的代码中，如果我将我的位置与exp(-f_/sigma2/sigma2)相乘，我会得到0(即使我的两个位置以及sigma都有非零值)，但当我将值相加时，我会得到正确的结果。

kernelsource = """  __kernel void forceFinder(
const int N,
const int dim,
const float sigma,
const float resistant,
__global float* datacl,
__constant float* poscl,
__global float* res

)
{
int i = get_global_id(0);
float f_sum ;
int k;
float sigma2 = sigma * sigma;
float tempo;
if (i < N ) {

    f_sum = 0;

    for (k = 0; k < dim; k++)
    {

        f_sum += pown((poscl[k] - datacl[i * dim + k]), 2);

    }




     for (k = 0; k < dim; k++)
    { 

        res[i * dim + k] = (datacl[i * dim + k] - poscl[k]) * exp(-f_sum/sigma2)/sigma2;

    }
  }

}
"""     

如果我将最后一个循环中的"*“替换为"+”，则会得到输出

Answer 1

我试着通过写一个完整的例子来理解你的内核在做什么(见下文)。虽然我不能完全理解发生了什么，但如果我在我的机器上运行下面的代码，我会收到一个10x10矩阵，所有条目的值都是-0.0024。如果您需要进一步的帮助，请提供完整的示例或更多信息作为您期望的结果。

另外:是否可以通过使用exp(x)*exp(x)=exp(x+x)来消除第一个for循环？

import pyopencl as cl
import numpy as np

kernelsource = """  
__kernel void forceFinder(  const int dim,
                            const float sigma,
                            __global float* datacl,
                            __constant float* poscl,
                            __global float* res ){
    int i = get_global_id(0);
    float f_sum = 0;
    float sigma2 = sigma * sigma;
    for (int k = 0; k < dim; k++){
        f_sum += pown((poscl[k] - datacl[i * dim + k]), 2);
    }

    for (int k = 0; k < dim; k++){ 
        res[i * dim + k] = (datacl[i * dim + k] - poscl[k]) * exp(-f_sum/sigma2)/sigma2;
    }
}
"""

device = cl.get_platforms()[0].get_devices()[0]
context = cl.Context([device])
program = cl.Program(context, kernelsource).build()
queue = cl.CommandQueue(context)

sigma = 20
dim = 10
N = 5
poscl_local = np.ones(dim).astype(np.float32) * 2.
datacl_local = np.ones((N,dim)).astype(np.float32)
res_local = np.zeros(datacl_local.shape).astype(np.float32)

poscl_buf = cl.Buffer(context, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=poscl_local)
datacl_buf = cl.Buffer(context, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=datacl_local)
res_buf = cl.Buffer(context, cl.mem_flags.WRITE_ONLY, res_local.nbytes)
program.forceFinder(queue,(N,), None, np.int32(dim), np.float32(sigma),datacl_buf,poscl_buf,res_buf)
cl.enqueue_copy(queue, res_local, res_buf)

print("result: {}".format(res_local))