ocl_prototype.cpp

Go to the documentation of this file.

#if defined(__cplusplus) && !defined(_MSC_VER)
extern "C" {
#endif
#include <clm4rm/ocl_prototype.h>
#include <m4ri/graycode.h>
#if defined(__cplusplus) && !defined(_MSC_VER)
}
#endif

#include <assert.h>
#include <gtest/gtest.h>
#include <omp.h>
#include <atomic>

//#define CEILCOLS(i)   ((i+clm4rm_radix-1)/clm4rm_radix)
//#define row(M,i)      (M + (i)*M ## _rowstride)
//#define at(M,row,col) *(M + (col)*M ## _nrows + row)

int nblocks(mzd_t const* A) {
    int blockrows = __M4RI_MAX_MZD_BLOCKSIZE / A->rowstride;
    return (A->nrows + blockrows - 1) / blockrows;
}

int adjust_k(int k, rci_t A_nrows) {
    size_t shared_mem_size = 32 * 1024;
    size_t shared_mem_words = shared_mem_size/sizeof(gpuword);
    //  TODO read actual size from OpenCL settings
    if (shared_mem_words < POW2(2) + 3 * A_nrows) {
        //  TODO when A_nrows becomes to large, we must split the matrix !!
        EXPECT_FALSE(true) << "too many rows " << A_nrows;
        return 2;
    }

    if (k==0) {
        //  find optimal k for given Cache size
       
        //  We need
        //  shared_mem_words >= 2^k + 3 * A_nrows;
        k = log2(shared_mem_words - 3 * A_nrows);
    }
    if (k < 2)
        k = 2;
//   else if (k > 8)
//        k = 8;
    
    EXPECT_GE(shared_mem_words, POW2(k) + 3 * A_nrows);
    return k;
}


bool assertEquals(const mzd_t* M, const gpuword* G, int padded_rows)
{
    int width = CEILCOLS(M->ncols);

    for (int row = 0; row < M->nrows; ++row)
    {
        word* Mrow = M->rows[row];
        for (int col = 0; col < width; col += 2)
        {
            if (col+1 < width)
                EXPECT_EQ(Mrow[col>>1], (word)G[col*padded_rows + row] | ((word)G[(col+1)*padded_rows + row]) << 32);
            else
                EXPECT_EQ(Mrow[col>>1], (word)G[col*padded_rows + row]);
        }
    }
    return true;
}


mzd_t* proto_bool_mul_m4rm(mzd_t* Cm, mzd_t const* Am, mzd_t const* Bm, int k)
{
    EXPECT_EQ(Am->ncols,Bm->nrows);
    EXPECT_TRUE(Cm==NULL || Am->nrows==Cm->nrows);
    EXPECT_TRUE(Cm==NULL || Bm->ncols==Cm->ncols);

    if (Cm == NULL)
        Cm = mzd_init(Am->nrows, Bm->ncols);
    else
        mzd_set_ui(Cm, 0);
    //  TODO add_mul

    //  We don't want blocking (resp: we must copy all blocks to GPU memory)
    EXPECT_LE(nblocks(Am), 1);
    EXPECT_LE(nblocks(Bm), 1);
    EXPECT_LE(nblocks(Cm), 1);
    
    k = adjust_k(k, Am->nrows);

    gpuword* A = copy_matrix_data(NULL,Am,Am->nrows);
    gpuword* B = copy_matrix_data(NULL,Bm,Bm->nrows);
    gpuword* C = copy_matrix_data(NULL,Cm,Cm->nrows);
    
    EXPECT_TRUE(assertEquals(Am, A,Am->nrows));
    EXPECT_TRUE(assertEquals(Bm, B,Bm->nrows));
    EXPECT_TRUE(assertEquals(Cm, C,Cm->nrows));

    proto_mul_m4rm( C,A,B, k,
                    Am->nrows, Am->ncols, Bm->ncols);

    copy_back_matrix_data(Cm,C,Cm->nrows);
    EXPECT_TRUE(assertEquals(Cm, C,Cm->nrows));

    delete A;
    delete B;
    delete C;

    return Cm;
}

/*
//  return the least significant bit
#define lsb(i) ((i) & -(i))

gpuword proto_read_bits(gpuword a0, gpuword a1, int spot, int n)
{
    int spill = spot + n - 32;
    gpuword temp;
    if (spill <= 0)
        temp = a0 << -spill;
    else
        temp = (a1 << (32 - spill)) | (a0 >> spill);
    return temp >> (32 - n);
}

void swap(gpuword** A, gpuword** B) {
    gpuword* aux = *A;
    *A = *B;
    *B = aux;
}

gpuword combinate(gpuword x, gpuword* T)
{
    gpuword result = 0;
    for (gpuword y = 1; x; y <<= 1, x >>= 1)
        result |= (x & 1) * T[y];
    //  Note: avoid if()
    return result;
}
/*
class Barrier {
private:
    int count;
    pthread_barrier_t barrier;
public:
    Barrier(int acount) : count(acount) {
        pthread_barrier_init(&barrier,NULL,count);
    }
    ~Barrier() {
        pthread_barrier_destroy(&barrier);
    }

    void reset() {
        pthread_barrier_destroy(&barrier);
        pthread_barrier_init(&barrier,NULL,count);
    }

    void operator() (void) {
        pthread_barrier_wait(&barrier);
    }
};
*/

#define A_width CEILCOLS(A_ncols)
#define C_ncols B_ncols
#define C_width CEILCOLS(C_ncols)
#define B_nrows A_ncols
#define C_nrows A_nrows

#define read(M,row,col)     M[col*M##_nrows+row]
#define write(M,row,col,x)  M[col*M##_nrows+row]=x

void proto_mul_m4rm_block(gpuword *C, const gpuword *A, const gpuword *B,
                          int k,
                          int A_nrows, int A_ncols, int B_ncols,
                          int r0, int r1);

void proto_mul_m4rm(gpuword *C, const gpuword *A, const gpuword *B, 
    int k,
    int A_nrows, int A_ncols, int B_ncols)
{
    k = 8;
    max_group_size = 256;

    int row_offset = A_nrows - (A_nrows % max_group_size);
    if (row_offset > 0)
        proto_mul_m4rm_block(C, A, B, k,
                         A_nrows, A_ncols, B_ncols,
                         0,row_offset);
    if (row_offset < A_nrows)
        proto_mul_m4rm_block(C, A, B, k,
                         A_nrows, A_ncols, B_ncols,
                         row_offset,A_nrows);
}

void proto_mul_m4rm_block(gpuword *C, const gpuword *A, const gpuword *B,
                    int k,
                    int A_nrows, int A_ncols, int B_ncols,
                    int r0, int r1)
{
    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
    ASSERT_TRUE((work_size[0] % group_size[0]) == 0);
    ASSERT_TRUE((work_size[1] % group_size[1]) == 0);

    size_t group_id_0,group_id_1;
    size_t group_offset_0,group_offset_1;
//    size_t global_id_0,global_id_1;
    size_t local_size_0,local_size_1;
//    size_t local_id_0,local_id_1;
    size_t local_count;
#define get_group_id(i) group_id_##i
#define get_global_id(i) global_id_##i
#define get_local_size(i) local_size_##i
#define get_local_id(i) local_id_##i
    group_offset_0=0;
    get_group_id(0)=0;
    while(group_offset_0 < work_size[0])
    {
        get_local_size(0) = MIN(group_size[0],work_size[0]-group_offset_0);
        group_offset_1=0;
        get_group_id(1)=0;

        while(group_offset_1 < work_size[1])
        {
            get_local_size(1) = MIN(group_size[1],work_size[1]-group_offset_1);
            local_count = get_local_size(0) * get_local_size(1);

            //
            //  1 GROUP
            //
            //  emulate shared memory
            gpuword* T = new gpuword[POW2(k)];

//            for(int get_local_id(0)=0; get_local_id(0) < get_local_size(0); ++get_local_id(0))
//                for(int get_local_id(1)=0; get_local_id(1) < get_local_size(1); ++get_local_id(1))
#pragma omp parallel num_threads(local_count)
            {
                int thread_num = omp_get_thread_num();

                int get_local_id(0) = thread_num % get_local_size(0);
                int get_local_id(1) = thread_num / get_local_size(0);

                int get_global_id(0) = group_offset_0+get_local_id(0);
                int get_global_id(1) = group_offset_1+get_local_id(1);

                EXPECT_EQ(get_local_size(1),1);
                EXPECT_EQ(get_local_id(1),0);
                EXPECT_EQ(get_group_id(1), get_global_id(1));

                //
                //  1 ITEM
                //
                int group_size = get_local_size(0);
                int ci = get_group_id(1);
                int cj = r0 + get_global_id(0);

                gpuword Csum = 0;
                gpuword A0 = read(A, cj, 0);
                gpuword A1 = read(A, cj, 1);

                //  iterate one row of A and one block column of B
                int ablock = 0;
                int aspot = 0;

                for (int ai = 0; ai < A_ncols; ai += k)
                {
                    int k1 = MIN(k, A_ncols - ai);
                    //  Make one column of T
                    //  distribute the 1 bit loop over the first k items
                    T[0] = 0;
                    for (int sj = 0; sj < k1; sj += group_size) {
                        int tj = sj + cj;//get_local_id(0);
                        if (tj < k1)
                            T[POW2(tj)] = read(B, ai + tj, ci);
                    }

#pragma omp barrier

                    for (int sj = 0; sj < POW2(k1); sj += group_size) {
                        int tj = sj + cj;//get_local_id(0);
                        if (tj < POW2(k1))
                            T[tj] = combinate(tj, T);
                    }

#pragma omp barrier

                    gpuword a = proto_read_bits(A0, A1, aspot, k1);
                    Csum |= T[a];

                    aspot += k;
                    if (aspot >= 32) {
                        //  cross over to next A column
                        aspot -= 32;
                        ablock++;
                        A0 = A1;

                        if ((ablock + 1) < A_width)
                            A1 = read(A, cj, ablock + 1);
                    }

#pragma omp barrier
                }
                //  write results back to global memory
                write(C, cj, ci, Csum);
                //  Note: thanks to column-major storage for A,B,C
                //  global memory access is coalesced, i.e. consecutive
                //  for all items in the group=column.
            }

            delete T;

            group_offset_1 += get_local_size(1);
            get_group_id(1)++;
        }

        group_offset_0 += get_local_size(0);
        get_group_id(0)++;
    }

    /*  Memory access pattern:
        C written once (for ci, for cj)
        B read once (for ci, for ai)
        A read once in each group

        =  n^2 * (2 + n/32)

        (where is the four-russian speed-up? it's gone!)

        TODO forget about the four russian and calculate all combinations from B
        TODO improve A reads by making the groups wider (fit a chunk of B into shared memory)
                but we need enough groups to keep the GPU busy
     */

}

mzd_t* proto_bool_mul_cubic(mzd_t* Cm, mzd_t const* Am, mzd_t const* Bm, int)
{
    EXPECT_EQ(Am->ncols, Bm->nrows);
    EXPECT_TRUE(Cm == NULL || Am->nrows == Cm->nrows);
    EXPECT_TRUE(Cm == NULL || Bm->ncols == Cm->ncols);

    EXPECT_GE(Am->nrows, 1);
    EXPECT_GE(Am->ncols, 1);
    EXPECT_GE(Bm->ncols, 1);

    if (Cm == NULL)
        Cm = mzd_init(Am->nrows, Bm->ncols);
    else
        mzd_set_ui(Cm, 0);
    //  TODO add_mul

    //  We don't want blocking (resp: we must copy all blocks to GPU memory)
    EXPECT_LE(nblocks(Am), 1);
    EXPECT_LE(nblocks(Bm), 1);
    EXPECT_LE(nblocks(Cm), 1);

    //  A and B and C are padded to multiple of 32 rows
    int rowpadding = 32;
    gpuword* A = copy_matrix_data(NULL, Am, padded_rows(Am->nrows, 32));
    gpuword* B = copy_matrix_data(NULL, Bm, padded_rows(Bm->nrows, 32));
    gpuword* C = copy_matrix_data(NULL, Cm, padded_rows(Cm->nrows, 32));

    EXPECT_TRUE(assertEquals(Am, A, padded_rows(Am->nrows, 32)));
    EXPECT_TRUE(assertEquals(Bm, B, padded_rows(Bm->nrows, 32)));
    EXPECT_TRUE(assertEquals(Cm, C, padded_rows(Cm->nrows, 32)));

    proto_mul_cubic(C, A, B,
        padded_rows(Am->nrows, rowpadding),
        Am->ncols, Bm->ncols);

    copy_back_matrix_data(Cm, C, padded_rows(Cm->nrows, 32));
    EXPECT_TRUE(assertEquals(Cm, C, padded_rows(Cm->nrows, 32)));

    delete A;
    delete B;
    delete C;

    return Cm;
}

void proto_mul_cubic(gpuword *C, const gpuword *A, const gpuword *B,
                    int A_nrows, int A_ncols, int B_ncols)
{
    max_group_size = 256;
    size_t work_size[2] = { (size_t)A_nrows, (size_t)C_width };
    size_t group_size[2] = { 32, 1 };

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

    size_t group_id_0,group_id_1;
    size_t group_offset_0,group_offset_1;
//    size_t global_id[2];
    size_t local_size_0,local_size_1;
//    size_t local_id[2];
    size_t local_count;
//    global_id[0]=global_id[1]=0;
    group_offset_0=0;
    get_group_id(0)=0;
    while(group_offset_0 < work_size[0])
    {
        get_local_size(0) = MIN(group_size[0],work_size[0]-group_offset_0);
        group_offset_1=0;
        get_group_id(1)=0;

        while(group_offset_1 < work_size[1])
        {
            get_local_size(1) = MIN(group_size[1],work_size[1]-group_offset_1);
            local_count = get_local_size(0) * get_local_size(1);

            EXPECT_EQ(get_local_size(0),32);
            EXPECT_EQ(get_local_size(1),1);
            EXPECT_EQ(local_count,32);
            //
            //  1 GROUP
            //
            //  emulate shared memory
            gpuword* T = new gpuword[32];

#pragma omp parallel num_threads(local_count)
            {
                int thread_num = omp_get_thread_num();

                int get_local_id(0) = thread_num % get_local_size(0);
                int get_local_id(1) = thread_num / get_local_size(0);

                EXPECT_EQ(get_local_size(0),32);
                EXPECT_EQ(get_local_size(1),1);
                EXPECT_EQ(get_local_id(1),0);

                int get_global_id(0) = group_offset_0+get_local_id(0);
                int get_global_id(1) = group_offset_1+get_local_id(1);

                //
                //  1 ITEM
                //
                int ci = get_global_id(1);
                int cj = get_global_id(0);
                int lcj = get_local_id(0);

                gpuword Csum = 0;                
                for (int ai = 0; ai < A_width; ++ai)
                {
                    //  Buffer 32 words of B in shared memory
                    T[lcj] = read(B, 32*ai + lcj, ci);
#pragma omp barrier
                    gpuword a = read(A, cj, ai);
                    for (int y = 0; y < 32; ++y, a >>= 1)
                        Csum |= (a & 1) * T[y];
#pragma omp barrier
                }
                //  write results back to global memory
                write(C, cj, ci, Csum);
                //  Note: thanks to column-major storage for A,B,C
                //  global memory access is coalesced, i.e. consecutive
                //  for all items in the group=column.
            }

            delete T;

            group_offset_1 += get_local_size(1);
            get_group_id(1)++;
        }

        group_offset_0 += get_local_size(0);
        get_group_id(0)++;
    }

    /*  Memory access pattern:
        C written once (for ci, for cj)
        B read once (for ci, for ai)
        A read once in each group

        =  n^2 * (2 + n/32)

        (where is the four-russian speed-up? it's gone!)

        TODO forget about the four russian and calculate all combinations from B
        TODO improve A reads by making the groups wider (fit a chunk of B into shared memory)
                but we need enough groups to keep the GPU busy
     */

}