Pycuda：多次调用内核的最佳方式

Question

我正在使用 pycuda 制作相对论光线追踪器。基本上，对于大型 2D 数组中的每个“像素”，我们必须使用 Runge Kutta 求解一个包含 6 个 ODE 的系统。由于每个集成都独立于其他集成，因此应该非常容易。其他人已经使用 C/C++ CUDA 实现了它，并取得了很好的效果（参见this project）。

问题在于我不知道如何做到这一点的最佳方式。我正在编写一个执行一些 Runge Kutta 步骤的内核，然后将结果返回给 CPU。这个内核被调用了很多次，以便整合整个射线。问题是由于某种原因非常慢。当然，我知道内存传输确实是 CUDA 中的瓶颈，但由于这真的很慢，我开始认为我做错了什么。

如果你能推荐我这个案例的最佳编程实践，那就太好了。 （使用pycuda）。我正在徘徊的一些事情：

1. 是否需要在到达内核调用时创建新上下文？
2. 有一种方法不必将内存从 GPU 传输到 CPU，即启动一个内核，暂停它以获取 一些信息，重新启动它并重复.
3. 每个 RK4 迭代大约需要半秒，这太疯狂了（也与链接中执行一些类似操作的 CUDA 代码相比）。而且我认为这是由于我使用 pycuda 的方式有问题，所以如果您能以最佳方式解释执行此类操作的最佳方法，那就太好了！。

澄清一下：我必须暂停/重新启动内核的原因是看门狗。超过 10 秒的内核被看门狗杀死。

提前谢谢你！

Answer 1

您的主要问题似乎过于笼统，没有看到代码很难给出一些具体的建议。我会尝试回答您的子问题（不是实际答案，但评论有点长）

我需要在到达内核调用时创建一个新的上下文吗？

没有。

有一种方法不必将内存从 GPU 传输到 CPU，即启动内核，暂停它以获取一些信息，然后重新声明并重复。

取决于“获取一些信息”的含义。如果这意味着在 CPU 上用它做一些事情，那么，当然，你必须转移它。如果你想在另一个内核调用中使用，那么你不需要转移它。

每个 RK4 迭代大约需要半秒，这太疯狂了（也与链接中执行一些类似操作的 CUDA 代码相比）。

这实际上取决于方程式、线程数和您使用的视频卡。我可以想象一个 RK 步骤需要这么长时间的情况。

而且我认为这是由于我使用 pycuda 的方式有问题，所以如果您能以最佳方式解释执行此类操作的最佳方式，那就太好了！。

没有代码就无法确定。尝试创建一些最小的演示示例，或者至少发布一个指向可运行（即使它相当长）的代码片段来说明您的问题。至于 PyCUDA，它是对 CUDA 的一个非常薄的包装，适用于后者的所有编程实践也适用于前者。

Answer 2

我可能会帮助您处理内存，即在您的迭代过程中不必从 CPU 复制到 GPU。我正在使用 euler timestepping 随着时间的推移发展一个系统，我将所有数据保存在我的 GPU 上的方式如下所示。 但是，这样做的问题是，一旦启动了第一个内核，cpu 就会继续执行它之后的行。 IE。边界内核在时间演化步骤之前启动。

我需要的是一种同步事物的方法。我尝试过使用 strm.synchronize() （请参阅我的代码），但它并不总是有效。如果您对此有任何想法，我将非常感谢您的意见！谢谢！

def curveShorten(dist,timestep,maxit):
"""
iterates the function image on a 2d grid through an euler anisotropic
diffusion operator with timestep=timestep maxit number of times
"""
image = 1*dist
forme = image.shape
if(np.size(forme)>2):
    sys.exit('Only works on gray images')

aSize = forme[0]*forme[1]
xdim  = np.int32(forme[0])
ydim  = np.int32(forme[1])  


image[0,:]      = image[1,:]
image[xdim-1,:] = image[xdim-2,:]
image[:,ydim-1] = image[:,ydim-2]
image[:,0]      = image[:,1]

#np arrays  i need to store things on the CPU, image is the initial 
#condition and final is the final state
image = image.reshape(aSize,order= 'C').astype(np.float32)
final = np.zeros(aSize).astype(np.float32)

#allocating memory to GPUs
image_gpu = drv.mem_alloc(image.nbytes)
final_gpu = drv.mem_alloc(final.nbytes)

#sending data to each memory location
drv.memcpy_htod(image_gpu,image) #host to device copying
drv.memcpy_htod(final_gpu,final)

#block size: B := dim1*dim2*dim3=1024
#gird size : dim1*dimr2*dim3 = ceiling(aSize/B)
blockX     = int(1024)
multiplier = aSize/float(1024)   
if(aSize/float(1024) > int(aSize/float(1024)) ):
    gridX = int(multiplier + 1)
else:
    gridX = int(multiplier)
strm1 = drv.Stream(1)
ev1   = drv.Event()
strm2 = drv.Stream()
for k in range(0,maxit):

    Kern_diffIteration(image_gpu,final_gpu,ydim, xdim, np.float32(timestep), block=(blockX,1,1), grid=(gridX,1,1),stream=strm1)
    strm1.synchronize()

    if(strm1.is_done()==1):
     Kern_boundaryX0(final_gpu,ydim,xdim,block=(blockX,1,1), grid=(gridX,1,1))
     Kern_boundaryX1(final_gpu,ydim,xdim,block=(blockX,1,1), grid=(gridX,1,1))#,stream=strm1)
     Kern_boundaryY0(final_gpu,ydim,xdim,block=(blockX,1,1), grid=(gridX,1,1))#,stream=strm2)
     Kern_boundaryY1(final_gpu,ydim,xdim,block=(blockX,1,1), grid=(gridX,1,1))#,stream=strm1)

    if(strm1.is_done()==1): 
      drv.memcpy_dtod(image_gpu, final_gpu, final.nbytes)
    #Kern_copy(image_gpu,final_gpu,ydim,xdim,block=(blockX,1,1), grid=(gridX,1,1),stream=strm1)



drv.memcpy_dtoh(final,final_gpu) #device to host copying
#final_gpu.free()
#image_gpu.free()


return final.reshape(forme,order='C')