import pycuda.autoinit
import pycuda.driver as drv
import numpy

from pycuda.compiler import SourceModule

mod = SourceModule("""
__global__ void get_distance(float* dest, float *a, float *b)
{
  
  const int i = threadIdx.x;
  dest[i] = (a[i] - b[i])*(a[i] - b[i]);

}
""")

get_distance = mod.get_function("get_distance")

a = numpy.random.randn(4000).astype(numpy.float32)
b = numpy.random.randn(4000).astype(numpy.float32)

dest = numpy.zeros_like(a)

get_distance(
        drv.Out(dest), drv.In(a), drv.In(b),
        block=(400,1,1), grid=(1,1))

assert numpy.allclose(dest, (a-b)*(a-b))
print(dest)
