/tomo/pyhst

#define PPT 4

// Iterations per line

#define IPT (16/PPT)

//#define PROJ_LINE (2 * IPT * PPT * 3)

#define PROJ_LINE (PPT * BLOCK_SIZE_X * 3 / 2)

 

__global__ static void hst_cuda_mplinear_kernel(int num_proj, int num_bins, float *d_SLICE, float apos_off_x, float apos_off_y, int batch) {

    float res1[PPT][PPT] = {0.f};

    float res2[PPT][PPT] = {0.f};

    

    const int tidx = threadIdx.x;

    const int tidy = threadIdx.y;

 

    const int bidx = PPT * blockIdx.x * BLOCK_SIZE_X;

    const int bidy = batch + PPT * blockIdx.y * BLOCK_SIZE_Y;

 

 

    const float bx = bidx + apos_off_x;

    const float by = bidy + apos_off_y;

 

    __shared__ float f_minh[128];

    __shared__ float2   cache[IPT * PROJ_LINE + 1];

    

    const int tbidy = tidy/(PPT/2);

    const int tbidx = tidy%(PPT/2);

    

    const int wrapid = (tidy&1) * 16 + tidx;

    

    const int sbidx = 8 * ((tidy / 2) % 2);

    const int sbidy = 8 * ((tidy / 2) / 2);

    

    const int stidx = tidx % 8;

    const int stidy = 2 * (2 * (tidy&1) + tidx / 8);

 

    const int sidx = (sbidx + stidx);

    const int sidy = (sbidy + stidy);

 

    const int idx = bidx + sidx;

    const int idy = bidy + sidy;

 

    const float x = idx + apos_off_x;

    const float y = idy + apos_off_y;

 

    const float exp23 = exp2(23.f);

    

    const int ttidx = 16 * tidy + tidx;

 

    float projf = tbidy + 0.5f;

    

    const int num_proj_blocks = num_proj / 32;

    

    for (int proj_block=0; proj_block<num_proj_blocks; proj_block += 4) {

        const int proj = proj_block * 32;

        

        float4 all = c_all[proj + ttidx];

        float minh = floor(all.z + bx * all.x - by * all.y + all.w);

        f_minh[ttidx] = minh;

        

        __syncthreads();

 

        float fminh[4];

 

#pragma unroll 4

        for (int i = 0; i < 4; i++) {

            fminh[i] = f_minh[32 * i + wrapid];

        }

 

        int max_proj = min(num_proj_blocks - proj_block, 4);

 

#pragma unroll 1

        for (int subproj32 = 0; subproj32 < max_proj; subproj32++) {

#pragma unroll 1

            for (int subproj = 0; subproj < (32/IPT); subproj++) {

                const int loop_proj = 32 * subproj32 + IPT * subproj;

                const int proje = proj + loop_proj;

 

                float4 all = c_all[proje + tbidy];

 

#ifdef HST_OPTIMIZE_KEPLER

                int minh = __shfl(fminh[subproj32], IPT * subproj + tbidy, 32);

#else //  HST_OPTIMIZE_KEPLER 

                int minh = i_minh[loop_proj + tbidy];

#endif //  HST_OPTIMIZE_KEPLER         

 

#pragma unroll 3

                for (int i = 0; i < 3; i++) {

                    int pos = (i * (PPT/2) + tbidx) * BLOCK_SIZE_X + tidx;

                

                    cache[PROJ_LINE * tbidy + pos].x = tex2D(tex_projes, minh + pos, proje + projf);

                }

            

            // we may use fence instead

            //__syncthreads();

 

#pragma unroll 3

                for (int i = 0; i < 3; i++) {

                    int pos = (i * (PPT/2) + tbidx) * BLOCK_SIZE_X + tidx;

                    cache[tbidy * PROJ_LINE + pos].y = cache[tbidy * PROJ_LINE + pos + 1].x - cache[tbidy * PROJ_LINE + pos].x;

                }

 

                __syncthreads();

 

#pragma unroll 4 // IPT

                for (int p = 0; p < IPT; p++) {

                    float4 all = c_all[proje + p];

 

#ifdef HST_OPTIMIZE_KEPLER

                    float minh = __shfl(fminh[subproj32], IPT * subproj + p, 32);

#else // HST_OPTIMIZE_KEPLER

                    float minh = f_minh[loop_proj + p];

#endif // HST_OPTIMIZE_KEPLER

 

                    float h = all.z + x * all.x - y * all.y - minh;

 

 

                    char *cache_row = ((char*)(&cache)) + p * PROJ_LINE * sizeof(float2);

 

 

#pragma unroll 4 // PPT

                    for (int i = 0; i < PPT; i++) {

                        float subh = h;

                        h -= 16.f * all.y;

#pragma unroll 4 // PPT

                        for (int j = 0; j < PPT; j++) {

                            float subh1 = subh;

                            float subh2 = subh - all.y;

 

 

                            float fsubh1 = subh1 + exp23;

                            int idx1 = (*(int*)(&fsubh1)) - 0x4B000000;

                            fsubh1 = subh1 - (fsubh1 - exp23);

 

                            float fsubh2 = subh2 + exp23;

                            int idx2 = (*(int*)(&fsubh2)) - 0x4B000000;

                            fsubh2 = subh2 - (fsubh2 - exp23);

 

 

                            float2 c1 = cache[p * PROJ_LINE + idx1];

                            res1[i][j] += c1.x + fsubh1 * c1.y;

 

                            float2 c2 = cache[p * PROJ_LINE + idx2];

                            res2[i][j] += c2.x + fsubh2 * c2.y;

 

                            subh += 16.f * all.x;

                        }

                    }

 

                }

            }

            __syncthreads();

        }

    }

 

#pragma unroll 4

    for (int i = 0; i < PPT; i++) {

#pragma unroll 4

        for (int j = 0; j < PPT; j++) {

            d_SLICE[BLOCK_SIZE_X * PPT * gridDim.x * (idy + i * BLOCK_SIZE_Y) + idx + j * BLOCK_SIZE_X] = res1[i][j];

            d_SLICE[BLOCK_SIZE_X * PPT * gridDim.x * (idy + i * BLOCK_SIZE_Y + 1) + idx + j * BLOCK_SIZE_X] = res2[i][j];

        }

    }

}