def test_take(self):
        idx = gpuarray.arange(0, 10000, 2, dtype=np.uint32)
        for dtype in [np.float32, np.complex64]:
            a = gpuarray.arange(0, 600000, dtype=np.uint32).astype(dtype)
            a_host = a.get()
            result = gpuarray.take(a, idx)

            assert (a_host[idx.get()] == result.get()).all()

    def test_fmod(self):
        """tests if the fmod function works"""
        for s in sizes:
            a = gpuarray.arange(s, dtype=np.float32)/10
            a2 = gpuarray.arange(s, dtype=np.float32)/45.2 + 0.1
            b = cumath.fmod(a, a2)

            a = a.get()
            a2 = a2.get()
            b = b.get()

            for i in range(s):
                assert math.fmod(a[i], a2[i]) == b[i]

    def test_ldexp(self):
        """tests if the ldexp function works"""
        for s in sizes:
            a = gpuarray.arange(s, dtype=np.float32)
            a2 = gpuarray.arange(s, dtype=np.float32)*1e-3
            b = cumath.ldexp(a, a2)

            a = a.get()
            a2 = a2.get()
            b = b.get()

            for i in range(s):
                assert math.ldexp(a[i], int(a2[i])) == b[i]

def bptrs(a):
    """
    Pointer array when input represents a batch of matrices.
    """

    return gpuarray.arange(a.ptr,a.ptr+a.shape[0]*a.strides[0],a.strides[0],
                dtype=cublas.ctypes.c_void_p)

    def test_arange(self):
        """test the arrangement of the array"""
        a = gpuarray.arange(12)

        res = a.get()

        for i in range(12):
            self.assert_(res[i] ==i)

    def test_abs(self):
        """test if the abs function works"""
        a = gpuarray.arange(111)
        a = a * -1

        res = a.get()

        for i in range (111):
            self.assert_(res[i] <= 0)

        a = abs(a)

        res = a.get()

        for i in range (111):
            self.assert_(res[i] >= 0)
            self.assert_(res[i] == i)


        for i in range(100,200):
            a = gpuarray.arange(500 * i)
            self.assert_(a[len(a)-1] == len(a)-1)

    def test():
        gpu_func = getattr(cumath, name)
        cpu_func = getattr(np, numpy_func_names.get(name, name))

        for s in sizes:
            for dtype in dtypes:
                args = gpuarray.arange(a, b, (b-a)/s, dtype=np.float32)
                gpu_results = gpu_func(args).get()
                cpu_results = cpu_func(args.get())

                max_err = np.max(np.abs(cpu_results - gpu_results))
                assert (max_err <= threshold).all(), \
                        (max_err, name, dtype)

    def test_abs(self):
        a = -gpuarray.arange(111, dtype=np.float32)
        res = a.get()

        for i in range(111):
            assert res[i] <= 0

        a = abs(a)

        res = a.get()

        for i in range(111):
            assert abs(res[i]) >= 0
            assert res[i] == i

    def test_frexp(self):
        """tests if the frexp function works"""
        for s in sizes:
            a = gpuarray.arange(s, dtype=np.float32)/10
            significands, exponents = cumath.frexp(a)

            a = a.get()
            significands = significands.get()
            exponents = exponents.get()

            for i in range(s):
                sig_true, ex_true = math.frexp(a[i])

                assert sig_true == significands[i]
                assert ex_true == exponents[i]

    def test_modf(self):
        """tests if the modf function works"""
        for s in sizes:
            a = gpuarray.arange(s, dtype=np.float32)/10
            fracpart, intpart = cumath.modf(a)

            a = a.get()
            intpart = intpart.get()
            fracpart = fracpart.get()

            for i in range(s):
                fracpart_true, intpart_true = math.modf(a[i])

                assert intpart_true == intpart[i]
                assert abs(fracpart_true - fracpart[i]) < 1e-4

    def test_sum_allocator(self):
        import pycuda.tools
        pool = pycuda.tools.DeviceMemoryPool()

        rng = np.random.randint(low=512,high=1024)

        a = gpuarray.arange(rng,dtype=np.int32)
        b = gpuarray.sum(a)
        c = gpuarray.sum(a, allocator=pool.allocate)

        # Test that we get the correct results
        assert b.get() == rng*(rng-1)//2
        assert c.get() == rng*(rng-1)//2

        # Test that result arrays were allocated with the appropriate allocator
        assert b.allocator == a.allocator
        assert c.allocator == pool.allocate

 def test_arange(self):
     a = gpuarray.arange(12, dtype=np.float32)
     assert (np.arange(12, dtype=np.float32) == a.get()).all()

def compute_pi(n):
    h = 1.0 / n
    x = h * (ga.arange(1, n, dtype=np.float32) + 0.5)
    s = ga.sum(4.0 / (1.0 + x**2), dtype=np.float32)
    return s.get() * h

ndev = 2
devlist = range(ndev)

# Setup the pycuda side
drv.init()
ctxs = [drv.Device(i).retain_primary_context() for i in devlist]

# Setup the communicator object
nc = NCCLComm(devlist)

# Now create gpuarrays for sending/recv buffers
srcs, dsts, size = [], [], 10

# Create some test arrays
for ctx in ctxs:
    ctx.push()
    srcs.append(gpuarray.arange(100, 200, size, dtype='<f4'))
    dsts.append(gpuarray.zeros((size,), dtype='<f4'))
    ctx.pop()

# Perform the reduction
nc.all_reduce(size, srcs, dsts)
nc.sync()

# Look at the results
for c, i, o in zip(ctxs, srcs, dsts):
    c.push()
    print i.get()
    print o.get()
    c.pop()

time = int(round(22.0 * xres / yres))
#print time
xlen = time / float(size[0])
#print xlen

#initialize out
out = np.zeros(0)

#rgb aliases
r=0
g=1
b=2

for x in range(xres):
    #float32 degrades quality, but float64 not support by gpus
    t_gpu = gpuarray.arange(x*xlen, x*xlen + xlen, 1./44100, dtype=np.float32)
    tone_gpu = gpuarray.zeros(t_gpu.size, dtype=np.float32)
    print "{0}%".format(round(100.0 * x / xres, 2))
    for y in range(yres):
        p = d[x+xres*y]
        #keep playing with these values
        amplitude = 10**(1-5.25+4.25*(p[r]+p[g]+p[b])/(255*3))
#        print amplitude, math.log(amplitude+1)
#        amplitude = math.log(amplitude+1)# / math.log(255)
#        print x, y, amplitude
        if p[r] > 10 or p[g] > 10 and p[b] > 10:
            tone_gpu += oscillator(t_gpu,
                                   amp = amplitude,
                                   #amp=(p[r]+p[g]+p[b]),
                                   freq=yscale * (yres - y))
    tone_gpu = tone_gpu + 1

import pycuda.gpuarray as gpuarray
import pycuda.autoinit
import numpy
from pycuda.reduction import ReductionKernel

vector_length = 400

input_vector_a = gpuarray.arange(vector_length, dtype=numpy.int)
input_vector_b = gpuarray.arange(vector_length, dtype=numpy.int)
dot_product = ReductionKernel(numpy.int,
                       arguments="int *x, int *y",
                       map_expr="x[i]*y[i]",
                       reduce_expr="a+b", neutral="0")

dot_product = dot_product(input_vector_a, input_vector_b).get()

print("INPUT MATRIX A")
print input_vector_a

print("INPUT MATRIX B")
print input_vector_b

print("RESULT DOT PRODUCT OF A * B")
print dot_product

        if stats_callback is not None:
            stats_callback(size,  self,
                    kernel_rec.kernel.prepared_timed_call(vectors[0]._grid, results[0]._block, *args))
        else:
            kernel_rec.kernel.prepared_async_call(vectors[0]._grid, results[0]._block, self.stream, *args)

        return results




if __name__ == "__main__":
    test_dtype = numpy.float32

    import pycuda.autoinit
    from pymbolic import parse
    expr = parse("2*x+3*y+4*z")
    print expr
    cexpr = CompiledVectorExpression(expr,
            lambda expr: (True, test_dtype),
            test_dtype)

    from pymbolic import var
    ctx = {
        var("x"): gpuarray.arange(5, dtype=test_dtype),
        var("y"): gpuarray.arange(5, dtype=test_dtype),
        var("z"): gpuarray.arange(5, dtype=test_dtype),
        }

    print cexpr(lambda expr: ctx[expr])

c = 5*a+6*b
e = time()
print 'cpu elapsed time: %f \n' % (e-s)

###################
# 4) map/reduce kernel

print '\n map/reduce kernel\n'
print '--------------------\n'

from pycuda.reduction import ReductionKernel

sz = 7000

# on device
a = gpuarray.arange(sz, dtype=numpy.float32)
b = gpuarray.arange(sz, dtype=numpy.float32)

krnl = ReductionKernel(numpy.float32, neutral="0",
        reduce_expr="a+2*b+b*b", map_expr="x[i] + y[i]",
        arguments="float *x, float *y")

# device perf
s = time()
my_dot_prod = krnl(a, b).get()
e = time()
print 'kernel time: %f' % (e-s)

# on host
a2 = arange(sz, dtype=numpy.float32)
b2 = arange(sz, dtype=numpy.float32)

 def test_take(self):
     idx = gpuarray.arange(0, 200000, 2, dtype=np.uint32)
     a = gpuarray.arange(0, 600000, 3, dtype=np.float32)
     result = gpuarray.take(a, idx)
     assert ((3*idx).get() == result.get()).all()

from pycuda.reduction import ReductionKernel
import pycuda.gpuarray as gpuarray
import pycuda.autoinit
import numpy

a = gpuarray.arange(400, dtype=numpy.float32)
b = gpuarray.arange(400, dtype=numpy.float32)

dot = ReductionKernel(dtype_out=numpy.float32, neutral="0",reduce_expr="a+b", map_expr="x[i]*y[i]",arguments="const float *x, const float *y")

a_dot_b = dot(a, b).get()
a_dot_b_cpu = numpy.dot(a.get(), b.get())