如何解决等效于C ++ Cuda中的一维numpy.interpCUDA中为Lerp
我有两个向量“ xp”和“ fp”,分别对应于数据的x和y值。第三个向量“ x”是我要评估插值的x坐标。我使用NumPy的interp函数在python中得到的结果与预期的一样。import numpy as np
xp = np.array([1.0,1.5,2.0,2.5,3.5,4.0,4.5,7.0,8.0,9.0,10.0,14.0,17.0,18.0,20.0])
yp = xp**2
x = np.array([3,5])
np.interp(x,xp,yp) #array([ 9.25,26. ])
我的问题是如何在cuda内核中复制此算法?
这是我的尝试:
大小-> len(xp)== len(fp), x_size-> len(x)。
__global__ void lerp_kernel(double* xp,double* yp,double* output_result,double* x,unsigned int size,unsigned int x_size)
{
for( unsigned int idx = blockIdx.x*blockDim.x + threadIdx.x ; idx < size ; idx +=
blockDim.x*gridDim.x )
{
if(idx > size - 2){
idx = size - 2;
}
double dx = xp[idx + 1] - xp[idx];
double dy = yp[idx + 1] - yp[idx];
double dydx = dy/dx ;
double lerp_result = yp[idx] + (dydx * (x[idx] - xp[idx]));
output_result[idx] = lerp_result;
}
}
我认为我犯的一个错误是我没有在xp中搜索包含x的索引范围(使用python中的numpy.searchsorted之类的东西)。我不确定如何在CUDA中实现这一部分。如果有人知道在cuda中做lerp的更好方法,请告诉我。
我在Cg(https://developer.download.nvidia.com/cg/lerp.html)文档中看到了lerp函数,但是对于x向量,它们需要权重向量(0-1之间的分数)。我不确定如何重新缩放x,以便可以使用权重向量来解决此问题。
解决方法
要模仿numpy.interp
的行为,将需要几个步骤。我们至少要做出一个简化的假设:numpy.interp
函数expects将xp
数组增加(我们也可以说“已排序”)。否则,它特别提到需要进行(内部)排序。如此处所示,我们将跳过这种情况,并假设您的xp
数组在增加。
从我所见,numpy函数还允许x
数组大致上是任意的。
为了进行适当的插值,我们必须找到每个xp
值所属的x
的“段”。我唯一想到的方法是执行binary search。 (还请注意,thrust具有便捷的二进制搜索功能)
该过程将是:
- 使用
x
中每个元素的二进制搜索,在xp
中找到相应的“段”(即端点) - 根据识别出的线段,使用直线方程(y = mx + b)计算端点之间的插值
这是一个例子:
$ cat t40.cu
#include <iostream>
typedef int lt;
template <typename T>
__device__ void bsearch_range(const T *a,const T key,const lt len_a,lt *idx){
lt lower = 0;
lt upper = len_a;
lt midpt;
while (lower < upper){
midpt = (lower + upper)>>1;
if (a[midpt] < key) lower = midpt +1;
else upper = midpt;
}
*idx = lower;
return;
}
template <typename T>
__global__ void my_interp(const T *xp,const T *yp,const lt xp_len,const lt x_len,const T *x,T *y){
for (lt i = threadIdx.x+blockDim.x*blockIdx.x; i < x_len; i+=gridDim.x*blockDim.x){
T val = x[i];
if ((val >= xp[0]) && (val < xp[xp_len - 1])){
lt idx;
bsearch_range(xp,val,xp_len,&idx);
T xlv = xp[idx-1];
T xrv = xp[idx];
T ylv = yp[idx-1];
T yrv = yp[idx];
// y = m * x + b
y[i] = ((yrv-ylv)/(xrv-xlv)) * (val-xlv) + ylv;
}
}
}
typedef float mt;
const int nTPB = 256;
int main(){
mt xp[] = {1.0,1.5,2.0,2.5,3.5,4.0,4.5,7.0,8.0,9.0,10.0,14.0,17.0,18.0,20.0};
lt xp_len = sizeof(xp)/sizeof(xp[0]);
mt *yp = new mt[xp_len];
for (lt i = 0; i < xp_len; i++) yp[i] = xp[i]*xp[i];
mt x[] = {3,5};
lt x_len = sizeof(x)/sizeof(x[0]);
mt *y = new mt[x_len];
mt *d_xp,*d_x,*d_yp,*d_y;
cudaMalloc(&d_xp,xp_len*sizeof(xp[0]));
cudaMalloc(&d_yp,xp_len*sizeof(yp[0]));
cudaMalloc(&d_x,x_len*sizeof( x[0]));
cudaMalloc(&d_y,x_len*sizeof( y[0]));
cudaMemcpy(d_xp,xp,xp_len*sizeof(xp[0]),cudaMemcpyHostToDevice);
cudaMemcpy(d_yp,yp,xp_len*sizeof(yp[0]),cudaMemcpyHostToDevice);
cudaMemcpy(d_x,x,x_len*sizeof(x[0]),cudaMemcpyHostToDevice);
my_interp<<<(x_len+nTPB-1)/nTPB,nTPB>>>(d_xp,d_yp,x_len,d_x,d_y);
cudaMemcpy(y,d_y,x_len*sizeof(y[0]),cudaMemcpyDeviceToHost);
for (lt i = 0; i < x_len; i++) std::cout << y[i] << ",";
std::cout << std::endl;
}
$ nvcc -o t40 t40.cu
$ cuda-memcheck ./t40
========= CUDA-MEMCHECK
9.25,26,========= ERROR SUMMARY: 0 errors
$
我并不是说上面的代码没有缺陷,也不适合任何特定目的。我的目的是演示一种方法,而不提供经过全面测试的代码。特别是,可能需要对边缘情况进行仔细测试(例如,值在整个范围的边缘或范围之外)。
版权声明:本文内容由互联网用户自发贡献,该文观点与技术仅代表作者本人。本站仅提供信息存储空间服务,不拥有所有权,不承担相关法律责任。如发现本站有涉嫌侵权/违法违规的内容, 请发送邮件至 dio@foxmail.com 举报,一经查实,本站将立刻删除。