PyOpenCL Kronecker Product Kernel
<p>I have the following code for confirming a hand written method for computing the kronecker product of two square matrices. The first portion indeed validates that my method of repeating and tiling <code>a</code> and <code>b</code> respectively yields the same output.</p>

<pre><code>import pyopencl as cl
import numpy
from time import time

N = 3
num_iter = 1

a = numpy.random.rand(N,N)
b = numpy.random.rand(N,N)
c = numpy.kron(a, b)

abig = numpy.repeat(numpy.repeat(a,N,axis=1),N,axis=0)
bbig = numpy.tile(b,(N,N))
cbig = abig*bbig

print(numpy.allclose(c,cbig))
</code></pre>

<p>I then attempt to port this multiplication over to the GPU using PyOpenCL. I first allocate <code>biga</code> and <code>bigb</code> as <code>d_a</code> and <code>d_b</code> respectively on the GPU memory. I also allocate <code>h_C</code> as an empty array on the host and <code>d_C</code> as the same size on the device.</p>

<pre><code>context = cl.create_some_context()
queue = cl.CommandQueue(context)

h_C = numpy.empty(cbig.shape)

d_a = cl.Buffer(context, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=abig)
d_b = cl.Buffer(context, cl.mem_flags.READ_ONLY | cl.mem_flags.COPY_HOST_PTR, hostbuf=bbig)
d_c = cl.Buffer(context, cl.mem_flags.WRITE_ONLY, h_C.nbytes)

kernelsource = open("../GPUTest.cl").read()
program = cl.Program(context, kernelsource).build()
kronecker = program.kronecker
kronecker.set_scalar_arg_dtypes([numpy.int32, None, None, None])

for i in range(num_iter):
    kronecker(queue, (N**2, N**2), None, N**2, d_a, d_b, d_c)
    queue.finish()

cl.enqueue_copy(queue, h_C, d_c)

print(h_C)
</code></pre>

<p>Here is the contents of GPUTest.cl:</p>

<pre><code>__kernel void kronecker(const int N,__global float* A,__global float*B,__global float* C)
{
    int i = get_global_id(0);
    int j = get_global_id(1);

    C[i,j] = A[i,j]*B[i,j];
}
</code></pre>

<p>However, my output is no where close. I believe my mistakes lie in how I'm handling the thread id's. From reading another example on matrix dot products, I was under the impression that the id's were essentially the location of the element within the block and since this is elementwise, I would only need to pull the element at the same location from A and B to multiply them together. Do these id's need to be combined into a single index to better address the way that the memory is actually allocated? </p>

<p>And only slightly related, but is there a way to utilize a tiling or memory sharing method? This was only a naiive attempt at the simplest way to do the calculation, I'm hoping to get to an algorithm that does not need the repeated/tiled versions of <code>a</code> and <code>b</code>. Something along the lines of taking a single element of <code>a</code>, multiplying the entirity of <code>b</code> by it, and then storing the result in a tile of <code>c</code>.</p>

<p>The issue is that the kernel does not take indices for the input and output memory addresses. The arguments should be <code>C[i+j*N]</code> in order to move throughout the whole block of memory appropriately.</p>

python opencl gpgpu matrix-multiplication pyopencl