【发布时间】：2012-10-13 10:39:07
【问题描述】：
我想在大型二维复杂数组（最终为 2*12x2*12 数据点）上运行一个函数。但是，pycuda 不能按预期工作。 ElementWise 函数不适用于二维数组，因此我使用了具有块大小的 SourceModule 函数。
现在的问题是，GPU 上的 C 代码给出的结果与 CPU 上的 numpy 计算结果不同。产生了非常大而奇怪的数字。
我正在使用以下代码。怎么了？
#!/usr/bin/env python
#https://github.com/lebedov/scikits.cuda/blob/master/demos/indexing_2d_demo.py
"""
Demonstrates how to access 2D arrays within a PyCUDA kernel in a
numpy-consistent manner.
"""
from string import Template
import pycuda
import pycuda.autoinit
import pycuda.gpuarray as gpuarray
from pycuda.compiler import SourceModule
import numpy as np
from matplotlib import pyplot as plt
# Set size
A = 2**3
B = 2**3
N = A*B
x_gpu = gpuarray.to_gpu(np.fromfunction(lambda x,y: (1.+x)*np.exp(1.j*y*np.pi/10), (A,B)) )
y_gpu = gpuarray.to_gpu(np.fromfunction(lambda x,y: 1.*x, (A,B)).astype(
x_gpu.dtype))
d_gpu = gpuarray.to_gpu(np.zeros_like(x_gpu.get()))#.astype(np.float32))
func_mod_template = Template("""
// Macro for converting subscripts to linear index:
#define INDEX(a, b) a*${B}+b
#include <pycuda-complex.hpp>
//__global__ void func(double *d,double *x,double *y, unsigned int N) {
__global__ void func(pycuda::complex<float> *d,pycuda::complex<float> *x,
pycuda::complex<float> *y)
{
// Obtain the linear index corresponding to the current thread:
// unsigned int idx = blockIdx.y*blockDim.y*gridDim.x +
blockIdx.x*blockDim.x*gridDim.y +threadIdx.x+threadIdx.y;
unsigned int block_num = blockIdx.x + blockIdx.y * gridDim.x;
unsigned int thread_num = threadIdx.y * blockDim.x + threadIdx.x;
unsigned int threads_in_block = blockDim.x * blockDim.y;
unsigned int idx = (threads_in_block * block_num + thread_num);
// Convert the linear index to subscripts:
unsigned int a = idx/${B};
unsigned int b = idx%${B};
// Use the subscripts to access the array:
// d[INDEX(a,b)] = x[INDEX(a,b)]+y[INDEX(a,b)];
pycuda::complex<float> j(0,arg(x[idx]));
pycuda::complex<float> i(abs(x[idx]),0);
d[idx] = i * exp(j);
}
""")
max_threads_per_block = pycuda.autoinit.device.get_attribute(pycuda._driver.device_attribute.MAX_THREADS_PER_BLOCK)
max_block_dim = (pycuda.autoinit.device.get_attribute(pycuda._driver.device_attribute.MAX_BLOCK_DIM_X),
pycuda.autoinit.device.get_attribute(pycuda._driver.device_attribute.MAX_BLOCK_DIM_Y),
pycuda.autoinit.device.get_attribute(pycuda._driver.device_attribute.MAX_BLOCK_DIM_Z))
max_grid_dim = (pycuda.autoinit.device.get_attribute(pycuda._driver.device_attribute.MAX_GRID_DIM_X),
pycuda.autoinit.device.get_attribute(pycuda._driver.device_attribute.MAX_GRID_DIM_Y),
pycuda.autoinit.device.get_attribute(pycuda._driver.device_attribute.MAX_GRID_DIM_Z))
max_blocks_per_grid = max(max_grid_dim)
block_dim = max_block_dim
block_dim = (max_block_dim[0],1,1)
grid_dim = (int(np.ceil(1.*x_gpu.shape[0]/block_dim[0])),
int(np.ceil(1.*x_gpu.shape[1]/block_dim[1])))
print block_dim,grid_dim, N
func_mod = \
SourceModule(func_mod_template.substitute(max_threads_per_block=max_threads_per_block,
max_blocks_per_grid=max_blocks_per_grid,
A=A, B=B))
func = func_mod.get_function('func')
func(d_gpu,x_gpu,y_gpu,
block=block_dim,
grid=grid_dim)
print d_gpu.get()/x_gpu.get()
#print 'Success status: ', np.allclose(x_np, x_gpu.get())
plt.imshow((d_gpu.get()/x_gpu.get()).real)
plt.colorbar()
plt.show()
【问题讨论】：
-
我不是专家，但我建议更明确地使用类型。确保 numpy 数据数组是 float32（您可以显式地为 numpy 数组指定类型），然后在内核中使用相同的类型。
-
我刚刚检查了代码，x_gpu 确实是使用 dtype complex128（默认）创建的。将其转换为 complex64 似乎可以解决问题。
-
根据@Bogdan，np.fromfunction(lambda x,y: (1.+x)*np.exp(1.jynp.pi/10), (A,B)).astype(np.complex64) 应该可以解决问题。始终，始终知道您传递到设备的数据的类型和顺序（C 与 Fortran）将使您在 Hotel PyCUDA 的逗留更加愉快。
-
请作为答案发布，而不是评论。