clm4rm_multiplication.cpp

Go to the documentation of this file.

//
// Created by nightrider on 21.09.18.
//

#include <clm4rm.h>
#include <stdlib.h>
#include <qdebug.h>
extern cl_kernel clm4rm_mul_kernel;

extern cl_kernel clcubic_mul_kernel[MAX_TILE_M+1];
extern cl_kernel clutri_mul_kernel[MAX_TILE_M+1];

void print_event_info(cl_event event);

void clm4rm_mul_block(clmatrix_t* C, clmatrix_t* A, clmatrix_t* B,
                          int r0, int r1,
                          cl_command_queue queue, clm4rm_conditions* cond);

void clm4rm_mul(clmatrix_t* C, clmatrix_t* A, clmatrix_t* B,
                    cl_command_queue queue, clm4rm_conditions* cond)
{
    int row_offset = A->nrows - (A->nrows % max_group_size);
    cl_event events[2];
    int e=0;
    if (row_offset > 0) {
        clm4rm_mul_block(C, A, B, 0, row_offset, queue, cond);
    }
    if (row_offset < A->nrows) {
        clm4rm_mul_block(C, A, B, row_offset, A->nrows, queue, cond);
    }
}

void clm4rm_mul_block(
                    clmatrix_t* C, clmatrix_t* A, clmatrix_t* B,
                    int r0, int r1,
                    cl_command_queue queue, clm4rm_conditions* cond)
{
    Q_ASSERT(A->ncols==B->nrows);
    Q_ASSERT(A->nrows==C->nrows);
    Q_ASSERT(B->ncols==C->ncols);
    Q_ASSERT(A->width*32 == B->padded_rows);
    Q_ASSERT(A->padded_rows==C->padded_rows);

    //  1 Work-Item = 1 Word in C
    cl_uint work_dim = 2;
    size_t work_size[2]  = { (size_t)r1-r0, (size_t)C->width };
    size_t group_size[2] = { MIN(work_size[0],max_group_size), 1 };

    //  group size must divide work size
    Q_ASSERT((work_size[0] % group_size[0]) == 0);
    Q_ASSERT((work_size[1] % group_size[1]) == 0);

    //  Work-Group = 1 Word-Column.
    //  Each Work-Group manages their Look-Up-Tables.

    int k=8; //TODO auto-choose

    /*
        __kernel void clm4rm_mul(
            __global gpuword* C,
            __global gpuword* A,
            __global gpuword* B,
            __local  gpuword T [2^k],
            int k, int r0,
            int A_nrows, int A_ncols, int B_ncols)
    */

    clm4rm_error = clSetKernelArg(clm4rm_mul_kernel, 0, sizeof(cl_mem), &C->data);
    clm4rm_error = clSetKernelArg(clm4rm_mul_kernel, 1, sizeof(cl_mem), &A->data);
    clm4rm_error = clSetKernelArg(clm4rm_mul_kernel, 2, sizeof(cl_mem), &B->data);
    clm4rm_error = clSetKernelArg(clm4rm_mul_kernel, 3, sizeof(gpuword)*POW2(k), NULL);
    clm4rm_error = clSetKernelArg(clm4rm_mul_kernel, 4, sizeof(int), &k);
    clm4rm_error = clSetKernelArg(clm4rm_mul_kernel, 5, sizeof(int), &r0);
    clm4rm_error = clSetKernelArg(clm4rm_mul_kernel, 6, sizeof(int), &A->padded_rows);
    clm4rm_error = clSetKernelArg(clm4rm_mul_kernel, 7, sizeof(int), &A->ncols);
    clm4rm_error = clSetKernelArg(clm4rm_mul_kernel, 8, sizeof(int), &B->ncols);

    clm4rm_error = clEnqueueNDRangeKernel(queue, clm4rm_mul_kernel,
                                          work_dim, NULL, work_size, group_size,
                                          pre_count(cond),pre_events(cond),push_event(cond));
    Q_ASSERT(clm4rm_error==CL_SUCCESS);
    Q_ASSERT(pushed_event(cond) != NULL);
}

void clcubic_mul_enqeue(clmatrix_t* C, clmatrix_t* A, clmatrix_t* B,
                            int tile_n, int tile_m,
                            size_t work_offset[2],
                            size_t work_size[2],
                            int uptri,
                            cl_command_queue queue, clm4rm_conditions* cond);

void clcubic_mul(clmatrix_t* C, clmatrix_t* A, clmatrix_t* B, size2_t max_tile,
    cl_command_queue queue, clm4rm_conditions* cond)
{
    //  A and B and C are expected to be padded to multiple of 32 rows
    Q_ASSERT(A->padded_rows % 32 == 0);
    Q_ASSERT(B->padded_rows % 32 == 0);
    Q_ASSERT(C->padded_rows % 32 == 0);

    Q_ASSERT(A->ncols == B->nrows);
    Q_ASSERT(A->nrows == C->nrows);
    Q_ASSERT(B->ncols == C->ncols);
    Q_ASSERT(A->padded_rows==C->padded_rows);
    Q_ASSERT(A->nrows==C->nrows);

    int tile_n = sqrt(max_group_size/32);
    tile_n = MAX(1, MIN(max_tile[0], tile_n));
    int tile_m = (sqrt(1 + 17*shared_mem_words)-1)/(68*tile_n);
    tile_m = MAX(1, MIN(max_tile[1], tile_m));
    Q_ASSERT(tile_m <= MAX_TILE_M);
    // limtied by shared memory size
    //  On Intel HD Graphis, n=2,m=2 seems to be OK, but m=3 is a disaster.
    //  Why is that so? Is the register file overbooked? Is shared memory overbooked?

    size2_t work_offset = { 0,0 };
    size2_t work_size_w, work_size_1, work_size_m, work_size_n, work_size, main_size_1;

    //  work-size in words
    work_size_w[0] = C->padded_rows;
    work_size_w[1] = C->width;
    //  work-size in 1x1 tiles (exact)
    work_size_1[0] = work_size_w[0]/32;
    work_size_1[1] = work_size_w[1];

    for(;;) {
        //  work-size in mxm tiles (small tiles)
        work_size_m[0] = work_size_1[0] / tile_m; // rounded down
        work_size_m[1] = work_size_1[1] / tile_m; // rounded down
        //  work-size in nxn tiles (big tiles)
        work_size_n[0] = work_size_m[0] / tile_n; // rounded down
        work_size_n[1] = work_size_m[1] / tile_n; // rounded down

        if ((work_size_n[0] > 0) && (work_size_n[1] > 0))
            break;

        //  else: reduce tile size
        if (tile_m > 1)
            tile_m--;
        else {
            Q_ASSERT(tile_n > 1);
            tile_n--;
        }
    }

    //  work-size in items
    work_size[0] = work_size_n[0]*32*tile_n;
    work_size[1] = work_size_n[1]*tile_n;
    //  group size is always { 32*tile_n,tile_n }

//    printf("   -- tile n=%i m=%i\n",tile_n, tile_m);
    clcubic_mul_enqeue(C, A, B, tile_n,tile_m, work_offset, work_size, 0, queue, cond);

    main_size_1[0] = work_size_n[0]*tile_n*tile_m;
    main_size_1[1] = work_size_n[1]*tile_n*tile_m;

    if (work_size_1[0] > main_size_1[0]) {
        // rest bottom with 1x1 tiles
        //  work offset in rows
        work_offset[0] = main_size_1[0]*32;
        work_offset[1] = 0;

        size2_t rest_work_size_1;
        rest_work_size_1[0] = work_size_1[0]-main_size_1[0];
        rest_work_size_1[1] = work_size_1[1];
        //  tile_n = tile_m = 1
        work_size[0] = rest_work_size_1[0]*32;
        work_size[1] = rest_work_size_1[1];

        clcubic_mul_enqeue(C, A, B, 1,1, work_offset, work_size, 0, queue, cond);
    }
    if (work_size_1[1] > main_size_1[1]) {
        // rest right with 1x1 tiles
        //  work offset in rows
        work_offset[0] = 0;
        work_offset[1] = main_size_1[1];

        size2_t rest_work_size_1;
        rest_work_size_1[0] = main_size_1[0];
        rest_work_size_1[1] = work_size_1[1]-main_size_1[1];

        work_size[0] = rest_work_size_1[0]*32;
        work_size[1] = rest_work_size_1[1];

        clcubic_mul_enqeue(C, A, B, 1,1, work_offset, work_size, 0, queue, cond);
    }
}

void print_event_info(cl_event event)
{
    cl_ulong start, end;

    clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_END,
        sizeof(cl_ulong), &end, NULL);
    clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_START,
        sizeof(cl_ulong), &start, NULL);

    float executionTimeInMilliseconds = (end - start) * 1.0e-6f;
    printf("  -- GPU clock = %lf ms\n", executionTimeInMilliseconds);
}

void clcubic_mul_enqeue(clmatrix_t* C, clmatrix_t* A, clmatrix_t* B,
                        int tile_n, int tile_m,
                        size_t work_offset[2],
                        size_t work_size[2],
                        int utri,
                        cl_command_queue queue, clm4rm_conditions* cond)
{
    int tile_width = tile_n*tile_m;
    int tile_ncols = tile_width*32;
    int tile_nrows = tile_n*tile_m*32;
    int col_stride = 34*tile_n*tile_m+1;

    cl_uint work_dim = 2;
    size2_t group_size = { (size_t)tile_n*32, (size_t)tile_n };

//  B buffer is column-major
//  A buffer is row-major
    size_t buf_bytes =  sizeof(gpuword)*tile_width*col_stride;
    Q_ASSERT(2*buf_bytes <= shared_mem_bytes);
    //  Note: buffered tiles have odd (prime) colstride to minimize bank conflicts (right?)

    /*
    __kernel void clcubic_mul(
            write_only_global C,
            read_only_global A,
            read_only_global B,
            __local  gpuword* A_buf,
            __local  gpuword* B_buf,
            int A_nrows, int A_ncols)
    */
    cl_kernel kernel = utri ? clutri_mul_kernel[tile_m] : clcubic_mul_kernel[tile_m];

    int p = 0;
    clm4rm_error = clSetKernelArg(kernel, p++, sizeof(cl_mem), &C->data);
    clm4rm_error = clSetKernelArg(kernel, p++, sizeof(cl_mem), &A->data);
    clm4rm_error = clSetKernelArg(kernel, p++, sizeof(cl_mem), &B->data);
#if BUFFERED
    clm4rm_error = clSetKernelArg(kernel, p++, buf_bytes, NULL);
    clm4rm_error = clSetKernelArg(kernel, p++, buf_bytes, NULL);
#endif
    clm4rm_error = clSetKernelArg(kernel, p++, sizeof(int), &A->padded_rows);

    if (!utri) {
        clm4rm_error = clSetKernelArg(kernel, p++, sizeof(int), &A->ncols);
    }

    //  TODO how about clEnqueueNDRangeKernel(...global_work_offset ? )
    //  and get_global_offset(0,1)
    //  Note that the cubic_mul kernel does not query get_global_id() at all.
    //  It applies an offset to get_group_id()*X instead.
    clm4rm_error = clEnqueueNDRangeKernel(queue, kernel,
                                          work_dim, work_offset, work_size, group_size,
                                          pre_count(cond),pre_events(cond),push_event(cond));
    switch (clm4rm_error) {
    case CL_MEM_OBJECT_ALLOCATION_FAILURE:
        printf( "OpenCL: CL_MEM_OBJECT_ALLOCATION_FAILURE. "
                " Allocated %s memory = %li \n", (IMAGE2D ? "texture":"global"), allocated_size);
        break;
    default:
        printf("OpenCL ERROR: %i \n", clm4rm_error);
        break;
    case CL_SUCCESS:
        break;
    }
    
    Q_ASSERT(clm4rm_error == CL_SUCCESS);
    Q_ASSERT(pushed_event(cond) != NULL);
}


void clutri_mul(clmatrix_t* C, clmatrix_t* A, clmatrix_t* B,
                size2_t max_tile,
                cl_command_queue queue, clm4rm_conditions* cond)
{
    //  A and B and C are expected to be padded to multiple of 32 rows
    Q_ASSERT(A->padded_rows % 32 == 0);
    Q_ASSERT(B->padded_rows % 32 == 0);
    Q_ASSERT(C->padded_rows % 32 == 0);

    Q_ASSERT(A->nrows==A->ncols);
    Q_ASSERT(B->nrows==B->ncols);
    Q_ASSERT(C->nrows==C->ncols);

    Q_ASSERT(A->ncols == B->nrows);
    Q_ASSERT(A->nrows == C->nrows);
    Q_ASSERT(B->ncols == C->ncols);
    Q_ASSERT(A->padded_rows==C->padded_rows);
    Q_ASSERT(A->nrows==C->nrows);

    //  fill C with Zero (not strictly necesary, but...
    gpuword zero=0;
    clEnqueueFillBuffer(
            queue, C->data, &zero, sizeof(zero),
            0, DATA_BYTES(C), 0, NULL, NULL);

    int tile_n = sqrt(max_group_size/32);
    tile_n = MAX(1, MIN(max_tile[0], tile_n));

    int tile_m = (sqrt(1 + 17*shared_mem_words)-1)/(68*tile_n);
    tile_m = MAX(1, MIN(max_tile[1], tile_m));
    Q_ASSERT(tile_m <= MAX_TILE_M);
    // limtied by shared memory size

    size2_t work_offset = { 0,0 };
    size2_t work_size_w, work_size_1, work_size_m, work_size_n, work_size, main_size_1;

    //  work-size in words
    work_size_w[0] = C->padded_rows;
    work_size_w[1] = C->width;
    //  work-size in 1x1 tiles (exact)
    work_size_1[0] = work_size_w[0]/32;
    work_size_1[1] = work_size_w[1];

    for(;;) {
        //  work-size in mxm tiles (small tiles)
        work_size_m[0] = work_size_1[0] / tile_m; // rounded down
        work_size_m[1] = work_size_1[1] / tile_m; // rounded down
        //  work-size in nxn tiles (big tiles)
        work_size_n[0] = work_size_m[0] / tile_n; // rounded down
        work_size_n[1] = work_size_m[1] / tile_n; // rounded down

        if ((work_size_n[0] > 0) && (work_size_n[1] > 0))
            break;

        //  else: reduce tile size
        if (tile_m > 1)
            tile_m--;
        else {
            Q_ASSERT(tile_n > 1);
            tile_n--;
        }
    }

    //  work-size in items
    work_size[0] = (work_size_n[0] - (work_size_n[0]-1)/2)*32*tile_n;
//  work_size[0] = work_size_n[0]*32*tile_n;
    //  fold lower part
    work_size[1] = work_size_n[1]*tile_n;
    //  group size is always { 32*tile_n,tile_n }

//  printf("   -- tile n=%i m=%i\n",tile_n, tile_m);
    clcubic_mul_enqeue(C, A, B, tile_n,tile_m, work_offset, work_size, 1, queue, cond);

    main_size_1[0] = work_size_n[0]*tile_n*tile_m;
    main_size_1[1] = work_size_n[1]*tile_n*tile_m;

    if (work_size_1[1] > main_size_1[1]) {
        // rest right with 1x1 tiles
        //  work offset in rows
        work_offset[0] = 0;
        work_offset[1] = main_size_1[1];

        size2_t rest_work_size_1;
        rest_work_size_1[0] = work_size_1[0];
        rest_work_size_1[1] = work_size_1[1]-main_size_1[1];

        work_size[0] = rest_work_size_1[0]*32;
        work_size[1] = rest_work_size_1[1];

        clcubic_mul_enqeue(C, A, B, 1,1, work_offset, work_size, 0, queue, cond);
    }
}


void printb(uint32_t x, int k) {
    uint32_t mask = (1<<k);
    for( ; mask; mask = mask>>1)
        printf("%i", (x & mask)?1:0);
}

void print3(uint32_t x, int k)
{
    printb(x,k);
    printf(" = ");
    //  LSB
    uint32_t lsb = x & -x;
    printb(lsb,k);
    printf(" + ");
    //  rest
    printb(x-lsb,k);
    printf("\n");
}


 void create_index_tables(uint32_t k)
{
    k++;
    for(int i=1; i < k; ++i) {
         //  generate all entries with Hamming weight i
         uint32_t val = (1 << i) - 1;
         print3(val,k+1);
         uint32_t stop = val << (k-i-1);
         //print3(stop,k+1);

         while(val != stop) {
             uint32_t c = val & -val; // lsb
             uint32_t r = val + c;
             val = (((r ^ val) >> 2) / c) | r;
             print3(val,k+1);
         }
     }
}