Pycuda:: 有一种方法可以将 float3 * 参数从具有 shape=(N,3) 的 numpy 数组传递给内核？答案

【问题标题】：Pycuda:: there is a way to pass a float3 * parameter to a kernel from a numpy array with shape=(N,3)?Pycuda:: 有一种方法可以将 float3 * 参数从具有 shape=(N,3) 的 numpy 数组传递给内核？
【发布时间】：2019-07-30 03:38:23
【问题描述】：

是否可以在 pycuda 内核中使用 shape (10, 3) 的 numpy 数组，就像 10 个 float3 的数组一样？

我正在尝试解决最近点的问题，在shape (10,3) 的点array_point 上有一个数组，其中10 个是点位置，例如array_point[0] 是[x,y,z]。

为了解决这个问题，我真的想向内核发送一个 float3* 参数，但我不知道该怎么做。

# for simplicity I will use a 4 point case with all handwritten directly
# only mockup script, not really working... actually is the question

import pycuda.driver as cuda
import pycuda.autoinit
from pycuda.compiler import SourceModule
from pycuda import gpuarray, tools

import numpy as np

data = np.array([[1,2,3],[3,4,5],[7,8,9], [10,11,12]], dtype=np.float32)
print(data.shape)
data_gpu = gpuarray.to_gpu(data)

// out must be like out_gpu[0] -> 1 means point 0 nearest point is point 1 ... I Hope be clear with the main idea

out_gpu = gpuarray.empty(4, np.int32)


mod = SourceModule("""
  __device__ float distance_not_sqrt(float3 p1, float3 p2)
  {
     return (p1.x - p2.x) * (p1.x - p2.x) + (p1.y - p2.y) * (p1.y - p2.y) + (p1.z - p2.z) * (p1.z - p2.z) ;
  }

  __global__ void find_closest(float3 *a, int*out)
  {
    int idx = threadIdx.x;
    int it;
    int it_min = -1;
    float dist_min = 1000.0; // more large than any real distance point
    for(it=0; it < 4; it++){
        if(it==idx)continue
        float dist = distance_not_sqrt(a[id], a[it])
        if(dist < dist_min){
            dist_min = dist;
            it_min = it;
        }
    }
    out[idx] = it_min;
  }
  """)

func = mod.get_function("find_closest")
func(data_gpu, out_gpu, block=(4,1,1))
print(out_gpu.get())

【问题讨论】：

标签： python pycuda

【解决方案1】：

传递 10 x 3 numpy 数组将自动适用于 float3 数组，但您需要确保您的连续维度（元素彼此相邻的维度）是其中包含 3 的维度。例如，像这样的 numpy 数组：

x = np.array([[1,2,3],[1,2,3],[1,2,3],[1,2,3]],dtype=np.float32,order='C')

与(1,2,3)的4个float3值相同，您可以使用np.ravel(order='K')进行检查

x.ravel(order='K')
array([1., 2., 3., 1., 2., 3., 1., 2., 3., 1., 2., 3.], dtype=float32)

但是如果改为使用 fortran 顺序（用“F”表示，其中“C”是 C 顺序），如果我们想要 4 (1,2,3)float3s，结果将不是我们所期望的。

x = np.array([[1,2,3],[1,2,3],[1,2,3],[1,2,3]],dtype=np.float32,order='F')

结果：

x.ravel(order='K')
array([1., 1., 1., 1., 2., 2., 2., 2., 3., 3., 3., 3.], dtype=float32)

这是因为我们的连续维度是 fortran 顺序中的第一个（即在x.shape == (4,3) 中，第一个维度 4 是我们在 fortran 顺序中的连续维度，最后一个维度是我们的C顺序的连续维度）

您的示例应该可以减去一些错误（例如 id 被使用但从未声明，假设您的意思是 idx？）

【讨论】：