clm4rm_mul.cl

Go to the documentation of this file.

typedef unsigned int gpuword;

#define CEILCOLS(i)     ((i+31)/32)
#if IMAGE2D
//
// Matrix stored in texture memory
//
# define read_only_global     __read_only image2d_t
# define write_only_global    __write_only image2d_t 
// Note: column-major format
// a matrix colum is actually a row (y-coordinate) in Image2D
// a matrix row is actually a column (x-coordinate) in Image2D
// Pixel contains only one (red) component
# define read(M,row,col)      read_imageui(M,(int2)(row,col)).x
# define write(M,row,col,x)   write_imageui(M,(int2)(row,col),(uint4)(x,0,0,0))
#else
//
// Matrix stored in __global memory
//
# define read_only_global     __global gpuword*
# define write_only_global    __global gpuword*                  
# define read(M,row,col)      M[(col)*M ## _nrows + row]
# define write(M,row,col,x)   M[(col)*M ## _nrows + row]=x
#endif

#define MIN(x,y)        (((x) < (y)) ? (x) : (y))
#define POW2(x)         (((gpuword)1) << x)

gpuword 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);
}

gpuword combinate(gpuword x, int k, __local gpuword* T)
{
    gpuword result = 0;
#pragma unroll
    for (int y = 0; y < k; ++y, x >>= 1)
        result |= (x & 1) * T[POW2(y)];
    return result;
}

 __kernel void clm4rm_mul(
            write_only_global C,
            read_only_global A,
            read_only_global B,
            __local  gpuword* T,
            int k, int r0,
            int A_nrows, int A_ncols, int B_ncols)
 {
#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

    //  work-group = column of C
    int group_size = get_local_size(0);
    // assert(group_size==A_nrows)
    int ci = get_group_id(1);

    //  work-item = one word of C
    int cj = r0 + get_global_id(0);
    int lcj = get_local_id(0);

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

    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+lcj;
            if (tj < k1)
                T[POW2(tj)] = read(B, ai+tj,ci);
        }

        //  read/write access to T must be guarded by barriers
        barrier(CLK_LOCAL_MEM_FENCE);

        //  Then calcluate the remaining (2^k-k) combinations; distribute over all items
        for (int sj=0; sj < POW2(k1); sj += group_size) {
            int tj = sj+lcj;
            if (tj < POW2(k1))
                T[tj] = combinate(tj, k1, T);
        }

        barrier(CLK_LOCAL_MEM_FENCE);

        //  apply table
        gpuword a = 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);
        }

        barrier(CLK_LOCAL_MEM_FENCE);
    }

    //  write result back to global memory
    write(C,cj,ci, Csum);
 }