OpenCL C Matrix Multiplication over Multiple Runs

Question

I am attempting to translate a dgemm / MPI matrix multiplier onto the GPU through OpenCL C. My issue is that the code below gives the correct output for a 9x6 matrix being multiplied by a 6x450 matrix until the 6x450 matrix changes slightly and the kernel output is shifted by 315 indeces. Included is the dgemm that I am trying to emulate for reference.

DGEMM:

    umAxBtrans(A, 450, 6, B, 9, 6, &fQ1);
    umAxBtrans(double **A, int Arows, int Acols, double *B, int Brows, int Bcols, double **C){
    op(A) = 'T';
    op(B) = 'F';
    dgemm(&opA, &opB, &rowsB, &rowsA, &colsA, &one, *B, &colsB, *A, &colsA, &zero, *C, &rowsB);

dgemm man page: http://www.math.utah.edu/software/lapack/lapack-blas/dgemm.html

Now if that dgemm implementation isn't a slight bit confusing then I don't know what is. I have taken care of the transpose of the B matrix outside of this matrix multiplication code because B is a static matrix. Matrix A, however, updates on every loop of the function call, meaning a new A is passed by value every time to openCL using the below function and kernel code.

OpenCL:

    clAxBtrans(A, 6, 450, Bcl, 9, 6, &C, context, program, queue, device_counter);
    clAxBtrans(double* A, int Arows, int Acols, double* Bcl, int Brows, int Bcols, double** C, cl_context context, cl_program program, cl_command_queue* queue, cl_uint device_counter){

In this function I pass the value of A with the dimensions used in the previous dgemm call (I've tried both passing A by reference and by value, no difference), and the static Bcl which is already transposed. I then setup the kernel and define:

    globalWorkSize[0]=450;
    globalWorkSize[1]=9;

and then write the kernel like so:

    __kernel void clAxBtrans(__global double* A,
                     __global double* B,
                     int rowsA,
                     int colsA,
                     int rowsB,
                     int colsB,
                     __global double* C)
    {
    int globalx = get_global_id(0);
    int globaly = get_global_id(1);
    double tmp = 0;
    double tmp0 = 0;
    double tmp1 = 0;
    double tmp2 = 0;
    double tmp3 = 0;
    double tmp4 = 0;
    double tmp5 = 0;
    tmp0 = B[globaly * colsB + 0] * A[0 * colsA + globalx];
    tmp1 = B[globaly * colsB + 1] * A[1 * colsA + globalx];
    tmp2 = B[globaly * colsB + 2] * A[2 * colsA + globalx];
    tmp3 = B[globaly * colsB + 3] * A[3 * colsA + globalx];
    tmp4 = B[globaly * colsB + 4] * A[4 * colsA + globalx];
    tmp5 = B[globaly * colsB + 5] * A[5 * colsA + globalx];
    tmp = tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;
    barrier(CLK_GLOBAL_MEM_FENCE);
    C[globaly * colsA + globalx] = tmp;
    }

The host code has blocking calls everywhere because this is in an MPI parallelized function and the localWorkSize is set as null to allow the gpu to choose the most "optimal" size.

The output is a 6x450 C matrix that matches the previous implementation until A updates. By that I mean for the first 5 runs A is a static matrix but then changes by a very small amount each run after depending on the resultant C matrix. This leads me to believe that I am interpreting matrix A incorrectly in memory. I am writing it to the kernel using the clBuffer call below:

    Acl = clCreateBuffer(context, CL_MEM_READ_ONLY | CL_MEM_COPY_HOST_PTR, sizeof(double)*Arows*Acols, A, &err);

For reference matrix B is the 9x6 below, A is an entire 6x450 filled with .111111 doubles, and the resultant matrix C is a 9x450 filled with .111111 doubles for the first 5 runs.

Matrix B:

    [1, 0 ,0 ,0 ,0 ,0]
    [0, 1 ,0 ,0 ,0 ,0]
    [0, 0 ,0 ,1 ,0 ,0]
    [0, 1 ,0 ,0 ,0 ,0]
    [0, 0 ,1 ,0 ,0 ,0]
    [0, 0 ,0 ,0 ,1 ,0]
    [1, 0 ,0 ,0 ,0 ,0]
    [0, 0 ,1 ,0 ,0 ,0]
    [0, 0 ,0 ,0 ,0 ,1]

Because I am copying the host ptr, I believe my problem either lies in how I think dgemm interprets the dimension flipping on the matrices, with my actual kernel code(most likely), or with how I am writing it as a buffer. I have checked and the very first new A is the same on both the opencl and previously working code, but the resultant C matrices are different(pointing to bad kernel code or memory management); I can include the resulting C comparison if required, but the first C results are shifted incorrectly by 315 indeces which causes the divergent solutions to really grow out of control.

I'm very grateful for any help you can afford me, even if it's pointing out that I've done something incredibly stupid! Thank you!

Answer 1

The answer to my question is that dgemm operates in a column major format, where C and the openCL kernel that I was implementing are row major. I changed my kernel code to the one below and that fixed everything. I found this out because I saw a Cblas_dgemm call where the first operator was CBLAS_ROW_MAJOR, telling the cblas solver to work in a row major format. This prompted me to do further research into the actual way dgemm works and found out that there is no option for making it row major in this implementation.

Kernel:

__kernel void clAxBtrans(__global double* A,
                     __global double* B,
                     int rowsA,
                     int colsA,
                     int rowsB,
                     int colsB,
                     __global double* C)
{
int globalx = get_global_id(0);
int globaly = get_global_id(1);
double tmp = 0;
double tmp0 = 0;
double tmp1 = 0;
double tmp2 = 0;
double tmp3 = 0;
double tmp4 = 0;
double tmp5 = 0;
tmp0 = B[0 * rowsB + globaly] * A[globalx * rowsA + 0];
tmp1 = B[1 * rowsB + globaly] * A[globalx * rowsA + 1];
tmp2 = B[2 * rowsB + globaly] * A[globalx * rowsA + 2];
tmp3 = B[3 * rowsB + globaly] * A[globalx * rowsA + 3];
tmp4 = B[4 * rowsB + globaly] * A[globalx * rowsA + 4];
tmp5 = B[5 * rowsB + globaly] * A[globalx * rowsA + 5];
tmp = tmp0 + tmp1 + tmp2 + tmp3 + tmp4 + tmp5;
C[globalx * rowsB + globaly] = tmp;
}

where the global variables are the same as described in the question.