Элементарная функция в массиве pycuda::complex

Я хочу запустить функцию для большого двухмерного сложного массива (в конечном итоге 2*12x2*12 точек данных). Однако pycuda не работает должным образом. Функция ElementWise не работает с двумерными массивами, поэтому я использовал функцию SourceModule с размерами блоков.

Теперь проблема в том, что код C на GPU не дает того же результата, что и вычисление numpy на CPU. Получаются очень большие и странные числа.

Я использую следующий код. Что происходит не так?

#!/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()

cuda pycuda
person Vasco    schedule 23.10.2012    source источник
comment
я не эксперт, но я бы посоветовал быть более точным с типами. убедитесь, что массив данных numpy имеет тип float32 (вы можете явно указать типы для массивов numpy), скажем, а затем используйте то же самое в ядре.   -  person andrew cooke    schedule 23.10.2012
comment
Я только что проверил код, и x_gpu действительно создан с dtype complex128 (по умолчанию). Приведение его к комплексу64, похоже, решает проблему.   -  person fjarri    schedule 23.10.2012
comment
Согласно @Bogdan, np.fromfunction(лямбда x,y: (1.+x)*np.exp(1.jynp.pi/10), (A,B)).astype (np.complex64) должен помочь. Всегда, всегда зная, какого типа и в каком порядке (C или Fortran) данные, которые вы передаете в устройство, сделают ваше пребывание в Hotel PyCUDA более приятным.   -  person Ahmed Fasih    schedule 23.10.2012
comment
Пожалуйста, напишите как ответ, а не комментарий.   -  person harrism    schedule 24.10.2012

Ответы (1)


arrow_upward
2
arrow_downward

Как фактический ответ: изменение строки x_gpu на

x_gpu = gpuarray.to_gpu(np.fromfunction(
    lambda x,y: (1.+x)*np.exp(1.j*y*np.pi/10), (A,B)).astype(np.complex64) )

кажется, решить проблему. Кроме того, хотя ElementwiseKernel не работает с 2d-массивами, вы все равно используете преобразование 2d->1d, поэтому ничто не мешает вам писать

func = ElementwiseKernel(
    "pycuda::complex<float> *d, pycuda::complex<float> *x, pycuda::complex<float> *y",

    Template("""
    // Convert the linear index to subscripts:
    unsigned int a = i/${B};
    unsigned int b = i%${B};

    // Use the subscripts to access the array:
    //d[INDEX(a,b)] = x[INDEX(a,b)]+y[INDEX(a,b)];
    pycuda::complex<float> angle(0,arg(x[i]));
    pycuda::complex<float> module(abs(x[i]),0);
    d[i] = module * exp(angle);
    """).substitute(A=A, B=B),

    preamble=Template("""
    #define INDEX(a, b) a*${B}+b
    """).substitute(A=A, B=B))

func(d_gpu, x_gpu, y_gpu)

Таким образом, вам не нужно жонглировать размерами блоков/сеток, потому что PyCUDA сделает это за вас.

person fjarri    schedule 24.10.2012