2
votes

In OpenCL, if I want to add two N-dimension vectors, the global work group size (globalSize) should satisfy globalSize = ceil(N/localSize) * localSize, where localSize is the local work group size. Is this correct? If N = 1000, and localSize = 128, globalSize should be 1024? Can we always set globalSize some multiple of localSize and larger than needed?

I tried many times and it worked well for 1-dimension problems.

However, when it comes to 2d problems, for example, multiply two matrices of dimension m*n and n*p, the result matrix is of order m*p, things get more complicated.

The max work group size on my device is 128, so I set localSize [2] = {16,8} and globalSize [2] = {ceil(m/16)*16,ceil(p/8)*8}.

It is similar to the 1-dimension case but the result is wrong!

If I set localSize [2] = {1,128} and change the globalSize accordingly, I can get the correct result. So where is the problem? Can anyone tell me why?

In addition, I find out the indices where the matrix element is wrong.

It seems that the result is wrong at (i,j) where i*p + j = n * some constant (n = 1,2,3...)

Why?

Here is my kernel function:

kernel void mmult(const int Mdim, const int Ndim, const int Pdim,
                  global float *A, global float *B, global float *C)
{

    int i = get_global_id(1);
    int j = get_global_id(0);
    if(i < 0 || j < 0 || i > Mdim || j > Pdim) return;
    else 
    {
        float tmp = 0;
        for(int k = 0; k < Ndim; k++)
            tmp += A[i*Ndim+k] * B[k*Pdim+j];

        C[i*Pdim + j] = tmp;
    }
}

And then it is the host program:

#define __NO_STD_VECTOR // Use cl::vector instead of STL version
#define __CL_ENABLE_EXCEPTIONS
#include <CL/cl.hpp>
#include <utility>
#include <iostream>
#include <fstream>
#include <string>
#include <cmath>
using namespace cl;

int main()
{
    // Create the two input matrices
    int m = 1000;
    int n = 1000;
    int p = 1000;
    float *A = new float[m*n];
    float *B = new float[n*p];
    for(int i = 0; i < m*n; i++)
    {
        A[i] = i;

    }
    for(int i = 0; i < n*p; i++)
    {
        B[i] = i;
    }
    try
    {
        // Get available platforms
        vector<Platform> platforms;
        Platform::get(&platforms);

        // Select the default platform and create a context using this platform and the GPU
        cl_context_properties cps[3] =
        {
            CL_CONTEXT_PLATFORM,
            (cl_context_properties)(platforms[0])(),
            0
        };
        Context context( CL_DEVICE_TYPE_GPU, cps);

        // Get a list of devices on this platform
        vector<Device> devices = context.getInfo<CL_CONTEXT_DEVICES>();

        // Create a command queue and use the first device
        CommandQueue queue = CommandQueue(context, devices[0]);

        // Read source file
        std::ifstream sourceFile("mmul.cl");
        std::string sourceCode(
            std::istreambuf_iterator<char>(sourceFile),
            (std::istreambuf_iterator<char>()));
        Program::Sources source(1, std::make_pair(sourceCode.c_str(), sourceCode.length()+1));

        // Make program of the source code in the context
        Program program = Program(context, source);

        // Build program for these specific devices

        program.build(devices);

        // Make kernel
        Kernel kernel(program, "mmult");

        // Create memory buffers
        Buffer bufferA = Buffer(context, CL_MEM_READ_ONLY, m*n * sizeof(float));
        Buffer bufferB = Buffer(context, CL_MEM_READ_ONLY, p*n * sizeof(float));
        Buffer bufferC = Buffer(context, CL_MEM_WRITE_ONLY, m*p * sizeof(float));

        // Copy lists A and B to the memory buffers
        queue.enqueueWriteBuffer(bufferA, CL_TRUE, 0, m * n * sizeof(float), A);
        queue.enqueueWriteBuffer(bufferB, CL_TRUE, 0, p * n * sizeof(float), B);

        // Set arguments to kernel
        kernel.setArg(0, m);
        kernel.setArg(1, n);
        kernel.setArg(2, p);
        kernel.setArg(3, bufferA);
        kernel.setArg(4, bufferB);
        kernel.setArg(5, bufferC);

        // Run the kernel on specific ND range

        NDRange global((ceil((float)(p)/16))*16,(ceil((float)(m)/8))*8);

        NDRange local(16,8);
        queue.enqueueNDRangeKernel(kernel, NullRange, global, local);

        // Read buffer C into a local list
        float *C = new float[m*p];
        queue.enqueueReadBuffer(bufferC, CL_TRUE, 0, m*p * sizeof(float), C);


        // check the correctness of the result
        float *c = new float[m*p];
        for(int i = 0; i < m; i++)
            for(int j = 0; j < p; j++)
            {
                float z = 0.0;
                for(int k = 0; k < n; k++)
                {
                    z += A[i*n+k] * B[k*p+j];
                }
                c[i*p+j] = z;
            }

        for(int i = 0; i < m*p; i++)
        {
            if(fabs(c[i]-C[i])>0.001)
                std::cout<<i<<" "<<c[i]<<" "<<C[i]<<std::endl;
        }

        delete []A;
        delete []B;
        delete []C;
    }
    catch(Error error)
    {
        std::cout << error.what() << "(" << error.err() << ")" << std::endl;
    }

    return 0;
}
2
Does your kernel check the global indices to ensure they are within the correct range? - jprice
I have if conditonal check, if the indices is out of proper range, then I just let the function return without doing any calculation. Like x = get_global_id(0), y = get_global_id(1); And only x between 0 and m and y between 0 and n, then the function does its calculation, otherwise it just returns. - iceiceice

2 Answers

2
votes

Your bounds checking code inside your OpenCL kernel is incorrect. Instead of this:

if(i < 0 || j < 0 || i > Mdim || j > Pdim) return;

You should have this:

if(i < 0 || j < 0 || i >= Mdim || j >= Pdim) return;
0
votes

Let's assume, that you have float matrix of size 1000x1000:

const int size = 1000;
// Whatever
float* myMatrix = (float*)calloc(size * size, sizeof(*myMatrix));

Determine size of Local Group first:

size_t localSize[] = {16, 8};

Then determine, how many Local Groups do you need:

size_t numLocalGroups[] = {ceil(size/localSize[0]), ceil(size/localSize[1])};

Finally, determine NDRange size:

size_t globalSize[] = {localSize[0] * numLocalGroups[0], localSize[1] * numLocalGroups[1]};

Don't forget to handle out-of-bounds access in right-most Local Groups.