OpenCL Error Computing Matrix Multiplication during Runtime

Question

I have been debugging for the past few days and cannot get this OpenCL matrix multiplication kernel to run. Whenever I run the program, the output from the GPU results in large negative numbers similar to -198746573.0000. I was wondering if someone with HPC experience could point out an error in my code or if it is an error with the driver.

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <string.h>

#define widthA 2
#define heightA 2

#define widthB heightA
#define heightB 2

#define widthC widthA
#define heightC heightB

#ifdef __APPLE__
#include < OpenCL/opencl.h >
#else
#include <opencl.h>
#endif

#define MEM_SIZE (128)
#define MAX_SOURCE_SIZE (0x100000)

int main()
{
  float * A = (float *)malloc(sizeof(float)*widthA*heightA);
  float * B = (float *)malloc(sizeof(float)*widthB*heightB);
  float * C = (float *)malloc(sizeof(float)*widthC*heightC);
  float * Res = (float *)malloc(sizeof(float)*widthC*heightC);
  float * D= (float *)malloc(sizeof(float)*widthC*heightC);

  float ref[widthC][heightC];

  int i, j, k;

   FILE * fp1 = fopen("matAdata.txt", "w");
  if (!fp1) {
    fprintf(stderr, "Failed to open matAdata.\n");
    exit(1);
   }

  for(i = 0;i < widthA; i++)
  {
        for(j=0;j < heightA; j++)       {
            float p=(rand()%100)/7.0;
            //*(A+i*heightA+j)=rand()%100 + p;
            *(A+i*heightA+j)=4.0;
            fprintf(fp1, "%f ",*(A+i*heightA+j));
        }
        fprintf(fp1, "\n");
   }
   fclose(fp1);

   fp1 = fopen("matBdata.txt", "w");
   if (!fp1) {
    fprintf(stderr, "Failed to open matAdata.\n");
    exit(1);
   }


    for(i = 0;i < widthB; i++)
    {
        for(j=0; j < heightB; j++)      {
            float p=(rand()%100)/7.0;
            //*((B+i*heightB+j))=rand()%100 + p;
            *((B+i*heightB+j))=4.0;
            fprintf(fp1, "%f ",*(B+i*heightA+j));
        }
        fprintf(fp1, "\n");
    }
    fclose(fp1);

  cl_device_id device_id = NULL;
  cl_context context = NULL;
  cl_command_queue command_queue = NULL;
  cl_mem memobjA = NULL;
  cl_mem memobjB = NULL;
  cl_mem memobjC = NULL;
  cl_mem rowA = NULL;
  cl_mem colC = NULL;
  cl_program program = NULL;
  cl_kernel kernel = NULL;
  cl_platform_id platform_id[10];
  cl_platform_id platform = NULL;
  cl_uint ret_num_devices;
  cl_uint ret_num_platforms;
  cl_int ret;
  cl_event GPUDone[0];
  //char string[MEM_SIZE];

  FILE *fp;
  char fileName[] = "matrixMultiplication.cl";
  char *source_str;
  size_t source_size;
  int row = widthA;
  int col = heightC;
  /* Load the source code containing the kernel*/
  fp = fopen(fileName, "r");
  if (!fp) {
    fprintf(stderr, "Failed to load kernel.\n");
    exit(1);
  }
  source_str = (char*)malloc(MAX_SOURCE_SIZE);
  source_size = fread( source_str, 1, MAX_SOURCE_SIZE, fp);
  fclose( fp );

  /* Get Platform and Device Info */
  ret = clGetPlatformIDs(10, platform_id, &ret_num_platforms);

  char cBuffer[1024];
  cl_uint c;

  for(c = 0; c < ret_num_platforms; c++)
  {
    clGetPlatformInfo(platform_id[c], CL_PLATFORM_NAME, 1024, &cBuffer, NULL);
    if (strstr(cBuffer, "NVIDIA") != NULL)
    {
        platform = platform_id[c];
        break;
    }

  }

  printf("Found Platform %s\n", cBuffer);

  ret = clGetDeviceIDs( platform, CL_DEVICE_TYPE_GPU, 1, &device_id, &ret_num_devices);

  printf("Found %d devices.\n", ret_num_devices);

  /* Create OpenCL context */
  context = clCreateContext( NULL, 1, &device_id, NULL, NULL, &ret);

  /* Create Command Queue */
  command_queue = clCreateCommandQueue(context, device_id, 0, &ret);

  /* Create Memory Buffer */
  memobjA = clCreateBuffer(context, CL_MEM_READ_ONLY, widthA * heightA * sizeof(float), NULL, &ret);
  memobjB = clCreateBuffer(context, CL_MEM_READ_ONLY, widthB * heightB * sizeof(float), NULL, &ret);
  memobjC = clCreateBuffer(context, CL_MEM_READ_WRITE, widthC * heightC * sizeof(float), NULL, &ret);
  rowA = clCreateBuffer(context, CL_MEM_READ_ONLY,  sizeof(int), NULL, &ret);
  colC = clCreateBuffer(context, CL_MEM_READ_ONLY,  sizeof(int), NULL, &ret);

  // Copy the lists A and B to their respective memory buffers
    ret = clEnqueueWriteBuffer(command_queue,memobjA, CL_TRUE, 0,
           widthA * heightA * sizeof(float), A, 0, NULL, NULL);
    ret = clEnqueueWriteBuffer(command_queue, memobjB, CL_TRUE, 0,
            widthB * heightB * sizeof(float), B, 0, NULL, NULL);
    ret = clEnqueueWriteBuffer(command_queue, rowA, CL_TRUE, 0, sizeof(int), &row, 0, NULL, NULL);
    ret = clEnqueueWriteBuffer(command_queue, colC, CL_TRUE, 0, sizeof(int), &col, 0, NULL, NULL);

  /* Create Kernel Program from the source */
  program = clCreateProgramWithSource(context, 1, (const char **)&source_str,
                                      (const size_t *)&source_size, &ret);

  /* Build Kernel Program */
  ret = clBuildProgram(program, 1, &device_id, NULL, NULL, NULL);

  /* Create OpenCL Kernel */
  kernel = clCreateKernel(program, "matrixMultiplication", &ret);

  /* Set OpenCL Kernel Arguments */
  ret = clSetKernelArg(kernel, 0, sizeof(cl_mem), (void *)&memobjA);
  ret = clSetKernelArg(kernel, 1, sizeof(cl_mem), (void *)&memobjB);
  ret = clSetKernelArg(kernel, 2, sizeof(cl_mem), (void *)&memobjC);
  ret = clSetKernelArg(kernel, 3, sizeof(int), (void *)&row);
  ret = clSetKernelArg(kernel, 4, sizeof(int), (void *)&col);
  /* Execute OpenCL Kernel */

  //ret = clEnqueueTask(command_queue, kernel, 0, NULL,NULL);
  size_t globalThreads[2] = {widthA, heightB};
  size_t localThreads[2] = {16,16};

  clEnqueueNDRangeKernel(command_queue, kernel, 2, NULL, globalThreads, localThreads, 0, NULL, NULL);
  //clFlush(command_queue);
  //clFinish(command_queue);

  /* Copy results from the memory buffer */
  ret = clEnqueueReadBuffer(command_queue, memobjC, CL_TRUE, 0,
                            widthA * heightC * sizeof(float), Res, 0, NULL, &GPUDone[0]);

  printf("Buffer Read ended with %d.\n", ret);
  clWaitForEvents(1, GPUDone);

  fp1 = fopen("matGPURes.txt", "w");
  if (!fp1) {
    fprintf(stderr, "Failed to open matAdata.\n");
    exit(1);
  }

  printf("\nResult\n");
    for(i = 0;i < widthA; i++)
    {
        for(j=0;j < heightC; j++)
        {

            fprintf(fp1, "%f ",*(Res+i*heightC+j));
            ref[i][j] = *(Res+i*heightC+j);
            printf("GPU Output: %f\n", *(Res+i*heightC+j));
        }
        fprintf(fp1, "\n");
    }
    fclose(fp1);

  ret = clFlush(command_queue);
  ret = clFinish(command_queue);
  ret = clReleaseKernel(kernel);
  ret = clReleaseProgram(program);
  ret = clReleaseMemObject(memobjA);
  ret = clReleaseMemObject(memobjB);
  ret = clReleaseMemObject(memobjC);
  ret = clReleaseCommandQueue(command_queue);
  ret = clReleaseContext(context);
  ret = clReleaseEvent(GPUDone[0]);

  free(source_str);

  float sum=0.0;

  for(i = 0;i < widthA; i++)
    {
        for(j = 0; j < heightC; j++)
        {
            sum = 0;
            for(k = 0; k < widthB; k++)
            {
                sum += A[i*col+k] * B[k*row+j];
                printf("Multiplying A: %f, B: %f\n", A[i*col+k], B[k*row+j]);
            }
        D[i*heightC+j] = sum;
        }

    }

    fp1 = fopen("matNormalMultiplicationRes.txt", "w");

  if (!fp1) {
    fprintf(stderr, "Failed to open matNormalMultiplicationRes.txt\n");
    exit(1);
  }

    for(i = 0; i<widthA; i++)
    {
        for(j = 0; j<heightA; j++)
        {
            if (ref[i][j] != D[i*heightA+j])
            {
                printf("Calculation error[ CPU: %f, GPU: %f ]\n", D[i*heightA+j], ref[i][j]);
            }
        }
    }

  printf("\nResult\n");
    for(i = 0;i < widthA; i++)
    {
        for(j=0;j < heightC; j++)
        {
            fprintf(fp1, "%f ",*(D+i*heightC+j));

        }
        fprintf(fp1, "\n");
    }
   free(A);
   free(B);
   free(C);
   free(D);
   free(Res);
  return 0;
}

Here is the kernel

#define BLOCK_SIZE 16

__kernel
void matrixMultiplication(__global float* A, __global float* B, __global float* C,  int wA, int wB )
{
    //int i = get_global_id(0);
    //int j = get_global_id(1);

    float Csub = 0.0f;        

    int bx = get_group_id(0);
    int by = get_group_id(1);

    int tx = get_local_id(0);
    int ty = get_local_id(1);

    int aBegin = wA * BLOCK_SIZE * by;
    int aEnd = aBegin + wA - 1;
    int aStep = BLOCK_SIZE;

    int bBegin = BLOCK_SIZE * bx;
    int bStep = BLOCK_SIZE * wB;

    for (int a = aBegin, b=bBegin;
        a <= aEnd;
        a += aStep, b+=bStep)
    {
        __local float As[BLOCK_SIZE][BLOCK_SIZE];
        __local float Bs[BLOCK_SIZE][BLOCK_SIZE];

        As[ty][tx] = A[a + wA * ty + tx];
        Bs[ty][tx] = B[b + wB * ty + tx];
        barrier(CLK_LOCAL_MEM_FENCE);

        for( int k = 0; k < BLOCK_SIZE; ++k)
            Csub += As[ty][k] * Bs[k][tx];
        barrier(CLK_LOCAL_MEM_FENCE);

    }

    int c = wB * BLOCK_SIZE * by + BLOCK_SIZE * bx;
    C[c + wB * ty + tx] = Csub;
    /*
    float value=0;
    for ( int k = 0; k < widthA; k++)
    {
        value = value + A[k + j * widthA] * B[k*widthB + i];
    }
    C[i + widthA * j] = value;
    */
}

I have double checked over and over again but simply cannot find any errors. I want to make sure its not a code error before I conclude its a driver issue.

Thanks!

Answer 1

Case is probably closed already, but for the sake of google-comers: Shouldnt shared memory be explicitly declared on host and passed as kernel argument to the source? __local keyword is not the one you are looking for in this case.

See post on How to declare local memory in OpenCL? for the detailed explanation.

Answer 2

Check the functionality of your host. Here a few things to get you started ...

1) You don't need to create a buffer and enqueue it for a scalar constant Int like row and col. Just set it as a kernel arg.

2) Wait for the clEnqueueNDRangeKernel with an event. You want to be sure the calc has completed.

3) Add a printf statement in the kernel to print selected values to see that the input and output values are what you expect.

try

if ( get_local_id(0) % 8 == 0) {

printf    some  useful  value of a,b,c

}

3) Try the host code with a dumb kernel that copies an input array to an output array. That will confirm it you have the handling of buffer creation and the enqeue read/write code correct!

Answer 3

Do you really need a complex kernel like that ? if you really want to do simple matrix multiplication you can write a simple kernel like this, which is easy to debug.

  __kernel void matrixMultiplication (__global float* A, 
                                      __global float* B,
                                      __global float* C,

                                       int widthA, int widthB )
{
    //y direction
    int row = get_global_id(1);

    int col = get_global_id(0);

    float cSum = 0.0f;

    //calculate the result
    for (int i=0; i<widthA; i++)
    {
        cSum += A[row*widthA+ i] * B[i*widthB+col];
    }

    C[row*widthB+col] = cSum;
}