PyOpenCL Kronecker Product Kernel

0
votes

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 a and b respectively yields the same output.

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))

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

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)

Here is the contents of GPUTest.cl:

__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];
}

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?

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 a and b. Something along the lines of taking a single element of a, multiplying the entirity of b by it, and then storing the result in a tile of c.

2
pythonopenclgpgpumatrix-multiplicationpyopencl
You're using singles in the openCL bit, and explicitly float32 in pyopencl. However, your numpy example doesn't set a dtype, so it should be float64 by default. Are you taking care of this when manipulating the data with pyopencl, and if not, does this change anything? - Andras Deak
Changing the kernel parameters to be __global double*'s successfully fixes the first row to be the same as the normal kronecker product but the remaining 8 rows are all 0. - Mike Vandenberg

2 Answers

0
votes

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

0
votes

I developed a kernel for the the kronecker product as well. I will put it here just fore reference. For

A # B = C,

where # is the the kronecker product, A is a m-x-n-matrix, B is a p-x-q-matrix and C is a x-x-y-matrix, with x=mp and y=nq, the follwoing kernel will calculate C:

__kernel void kroneckerProdFast(__global float* res,
                                __global float* a,
                                __global float* b,
                                int p,
                                int q){
    int xi = get_global_id(0);
    int x = get_global_size(0); 
    int yi = get_global_id(1);
    int y = get_global_size(1); 

    int n = y / q;
    int mi = xi / p;
    int ni = yi / q;
    int pi = xi % p;
    int qi = yi % q;

    res[xi * y + yi] = a[mi * n + ni] * b[pi * q + qi];
}

The call from PyOpenCL would be:

program.kroneckerProdFast(queue,res.shape, None, resf_buf, a_buf, b_buf,np.int32(b.shape[0]),np.int32(b.shape[1]))