/tomo/pyhst

/*

 * The PyHST program is Copyright (C) 2002-2011 of the

 * European Synchrotron Radiation Facility (ESRF) and

 * Karlsruhe Institute of Technology (KIT).

 *

 * PyHST is free software: you can redistribute it and/or modify it

 * under the terms of the GNU General Public License as published by the

 * Free Software Foundation, either version 3 of the License, or

 * (at your option) any later version.

 *

 * hst is distributed in the hope that it will be useful, but

 * WITHOUT ANY WARRANTY; without even the implied warranty of

 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.

 * See the GNU General Public License for more details.

 *

 * You should have received a copy of the GNU General Public License along

 * with this program.  If not, see <http://www.gnu.org/licenses/>.

 */

 

 

__constant__ float4  c_all[ MAXNPROJECTIONS  ] ;        //!< (cosinuses, sinuses, axis position, min bin

texture<float, 2, cudaReadModeElementType> tex_projes;

 

 

__global__ static void  hst_cuda_unpack_kernel_fai360(cufftReal *out, int dpitch, cufftReal *data, int spitch, int half, float *params, int batch) {

    float tmp1, tmp2;

 

    int  idx = blockIdx.x * BLOCK_SIZE + threadIdx.x;

    int  idy = blockIdx.y * BLOCK_SIZE + threadIdx.y;

 

    int dest_vector = (idx >= half)?1:0;

    int dest_pos = 2 * (idx - dest_vector * half);

 

    int src_pos = 2 * (idy * spitch + dest_pos) + dest_vector;

    int src_vector = dest_vector;//(src_pos >= 0);

 

    tmp1 = data[src_pos];

    tmp2 = data[src_pos + 2];

 

    float axis_position_corr = c_all[batch + 2*idy + src_vector].z;

    float flat_zone = params[idy * 4 + 2*src_vector];

    float param = params[idy * 4 + 2*src_vector + 1];

 

    int pos = (2 * idy + dest_vector)*dpitch + dest_pos;

 

    float apc_plus_flat = axis_position_corr + flat_zone;

    float apc_param_flat = apc_plus_flat + param - dest_pos;

 

    int flag1 = dest_pos < apc_plus_flat;

    float multiplier1 = __fdividef((flag1?(apc_param_flat - 2 * flat_zone):apc_param_flat),param);

 

    int flag2 =(dest_pos + 1) < apc_plus_flat;

    float multiplier2 = __fdividef((flag2?(apc_param_flat - 2 * flat_zone):apc_param_flat) - 1,param);

 

    if (param > 0) {

        out[pos] = tmp1 * min(2., max(flag1?1.:0., multiplier1));

        out[pos+1] = tmp2 * min(2., max(flag2?1.:0., multiplier2));

    } else {

        out[pos] = tmp1 * max(0., min(flag1?1.:2., multiplier1));

        out[pos+1] = tmp2 * max(0., min(flag2?1.:2., multiplier2));

    }

}

 

 

 

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

    const int idx = blockIdx.x * BLOCK_SIZE_X + threadIdx.x;

    const int idy = blockIdx.y * BLOCK_SIZE_Y + threadIdx.y + batch;

 

    float h;

    float res=0;

 

    const float bx = idx + apos_off_x;

    const float by = idy + apos_off_y;

 

#pragma unroll 8

    for (int proje=0; proje<num_proj; proje++) {

        const float4 all = c_all[proje];

        h = all.z + bx * all.x - by * all.y;

 

        res=res+tex2D(tex_projes, h, proje + 0.5f) ;

    }

 

    d_SLICE[ BLOCK_SIZE_X*gridDim.x*idy + idx ] = res;

}

 

 

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

    float h;

    float res[4][4] = {0};

 

#ifdef HST_OPTIMIZE_KEPLER

    __shared__ float buf[16][18]; // 64b for kepler

    __shared__ float fill[48];

    __shared__ float fin[16][18];

#else // HST_OPTIMIZE_KEPLER

    __shared__ float buf[16][17]; // 32b for Fermi & GT200

    __shared__ float fill[56];

    __shared__ float fin[16][17];

#endif // HST_OPTIMIZE_KEPLER

 

    const int tidx = threadIdx.x;

    const int tidy = threadIdx.y;

 

    const int bidx = blockIdx.x * BLOCK_SIZE_X;

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

 

    const int sidx = tidx % 4;

    const int sidy = tidx / 4;

 

    const float x = bidx + sidx + apos_off_x;

    const float y = bidy + sidy + apos_off_y;

 

    const float projf = tidy + 0.5f;

 

    for (int proje=0; proje<num_proj; proje+=16) {

        const float4 all = c_all[proje+tidy];

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

 

#pragma unroll 4

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

#pragma unroll 4

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

                float subh = h + 4 * j * all.x - 4 * i * all.y;

                res[i][j] += tex2D(tex_projes, subh, proje + projf);

            }

        }

 

    }

 

 

#pragma unroll 4

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

#pragma unroll 4

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

            buf[tidx][tidy] = res[i][j];

 

            __syncthreads();

 

 

#ifdef HST_OPTIMIZE_KEPLER

            float val = buf[tidy][tidx];

            for (int k=16; k>=1; k/=2)

                val += __shfl_xor(val, k, 16);

#else // HST_OPTIMIZE_KEPLER

            volatile float *ptr = &buf[tidy][0];

            for (int k=8; k>=1; k/=2)

                ptr[tidx] += ptr[tidx+k];

            float val = ptr[0];

#endif // HST_OPTIMIZE_KEPLER

 

            const int rx = 4 * j + tidy%4;

            const int ry = 4 * i + tidy/4;

 

            if (!tidx) {

                fin[ry][rx] = val;

            }

 

            __syncthreads();

        }

    }

 

 

    const int idx = bidx + tidx;

    const int idy = bidy + tidy;

 

    d_SLICE[BLOCK_SIZE_X * gridDim.x * idy + idx] = fin[tidy][tidx];

}

 

 

 

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

    float2 res[2] = {0.f, 0.f, 0.f, 0.f};

    

    const int tidx = threadIdx.x;

    const int tidy = threadIdx.y;

 

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

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

 

    const int idx = bidx + 2 * tidx;

    const int idy = bidy + 2 * tidy;

 

    const float bx = bidx + apos_off_x;

    const float by = bidy + apos_off_y;

 

    const float x = idx + apos_off_x;

    const float y = idy + apos_off_y;

 

 

    __shared__  float   cache_subh[2*BLOCK_SIZE_X][BLOCK_SIZE_Y + 1];

    __shared__  float   cache_sin[BLOCK_SIZE_Y];

    __shared__  float   cache[BLOCK_SIZE_Y][3 * BLOCK_SIZE_X];

    

//#pragma unroll 8

    float projf = tidy + 0.5;

    for (int proje=0; proje<num_proj; proje += BLOCK_SIZE_Y) {

        float4 all = c_all[proje + tidy];

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

 

        float subh = (all.z + x * all.x - minh);

        cache_subh[2*tidx][tidy] = subh;

        cache_subh[2*tidx+1][tidy] = subh + all.x;

        if (tidx == 0) cache_sin[tidy] = all.y;

 

 

#pragma unroll 3

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

            int pos = i * BLOCK_SIZE_X + tidx;

            cache[tidy][pos] = tex2D(tex_projes, minh + pos, proje + projf);

        }

 

        __syncthreads();

 

#pragma unroll 16

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

            int4 minh = {

                cache_subh[2 * tidx][i] - (y) * cache_sin[i],

                cache_subh[2 * tidx + 1][i] - (y) * cache_sin[i],

                cache_subh[2 * tidx][i] - (y + 1) * cache_sin[i],

                cache_subh[2 * tidx + 1][i] - (y + 1) * cache_sin[i],

            };

 

            res[0].x += cache[i][minh.x];

            res[0].y += cache[i][minh.y];

            res[1].x += cache[i][minh.z];

            res[1].y += cache[i][minh.w];

 

        }

 

 

        __syncthreads();

 

    }

 

    *(float2*)(d_SLICE + BLOCK_SIZE_X * 2 * gridDim.x * (idy    ) + idx) = res[0];

    *(float2*)(d_SLICE + BLOCK_SIZE_X * 2 * gridDim.x * (idy + 1) + idx) = res[1];

}

 

 

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

    float4 res = {0.f, 0.f, 0.f, 0.f};

    

    const int tidx = threadIdx.x;

    const int tidy = threadIdx.y;

 

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

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

 

    const int idx = bidx + 2 * tidx;

    const int idy = bidy + 2 * tidy;

 

    const float bx = bidx + apos_off_x;

    const float by = bidy + apos_off_y;

 

    const float x = idx + apos_off_x;

    const float y = idy + apos_off_y;

 

 

    __shared__  float   cache_subh[2*BLOCK_SIZE_X][BLOCK_SIZE_Y + 1];

    __shared__  float   cache_sin[BLOCK_SIZE_Y];

    __shared__  float2   cache[BLOCK_SIZE_Y][3 * BLOCK_SIZE_X + 1];

    

//#pragma unroll 8

    float projf = tidy + 0.5;

    for (int proje=0; proje<num_proj; proje += BLOCK_SIZE_Y) {

        float4 all = c_all[proje + tidy];

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

 

        float subh = (all.z + x * all.x - minh);

        cache_subh[2*tidx][tidy] = subh;

        cache_subh[2*tidx+1][tidy] = subh + all.x;

        if (tidx == 0) cache_sin[tidy] = all.y;

 

 

#pragma unroll 3

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

            int pos = i * BLOCK_SIZE_X + tidx;

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

        }

 

        __syncthreads();

 

#pragma unroll 3

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

            int pos = i * BLOCK_SIZE_X + tidx;

            cache[tidy][pos].y = cache[tidy][pos + 1].x - cache[tidy][pos].x;

        }

 

        __syncthreads();

 

#pragma unroll 16

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

            float2 c;

            float4 v1, v2;

            float4 minh = {

                cache_subh[2 * tidx][i] - (y) * cache_sin[i],

                cache_subh[2 * tidx + 1][i] - (y) * cache_sin[i],

                cache_subh[2 * tidx][i] - (y + 1) * cache_sin[i],

                cache_subh[2 * tidx + 1][i] - (y + 1) * cache_sin[i],

            };

 

            c = cache[i][(int)minh.x];

            v1.x = c.x; v2.x = c.y;

            c = cache[i][(int)minh.y];

            v1.y = c.x; v2.y = c.y;

            c = cache[i][(int)minh.z];

            v1.z = c.x; v2.z = c.y;

            c = cache[i][(int)minh.w];

            v1.w = c.x; v2.w = c.y;

            

            res += v1 + (minh - floorf(minh)) * v2;

        }

 

        __syncthreads();

 

    }

 

    *(float2*)(d_SLICE + BLOCK_SIZE_X * 2 * gridDim.x * (idy    ) + idx) = make_float2(res.x, res.y);

    *(float2*)(d_SLICE + BLOCK_SIZE_X * 2 * gridDim.x * (idy + 1) + idx) = make_float2(res.z, res.w);

}

 

 

 

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

    float2 res[2] = {0.f, 0.f, 0.f, 0.f};

    

    const int tidx = threadIdx.x;

    const int tidy = threadIdx.y;

 

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

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

 

    const int idx = bidx + 2 * tidx;

    const int idy = bidy + 2 * tidy;

 

    const float bx = bidx + apos_off_x;

    const float by = bidy + apos_off_y;

 

    const float x = idx + apos_off_x;

    const float y = idy + apos_off_y;

 

 

    __shared__  float   cache_subh[2*BLOCK_SIZE_X][BLOCK_SIZE_Y + 1];

    __shared__  float   cache_sin[BLOCK_SIZE_Y];

    __shared__  float   cache[BLOCK_SIZE_Y][12 * BLOCK_SIZE_X];

    

//#pragma unroll 8

    float projf = tidy + 0.5;

    for (int proje=0; proje<num_proj; proje += BLOCK_SIZE_Y) {

        float4 all = c_all[proje + tidy];

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

 

        float subh = 4 * (all.z + x * all.x - minh);

        cache_subh[2*tidx][tidy] = subh;

        cache_subh[2*tidx+1][tidy] = subh + 4 * all.x;

        if (tidx == 0) cache_sin[tidy] = 4 * all.y;

 

 

#pragma unroll 12

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

            int pos = i * BLOCK_SIZE_X + tidx;

            cache[tidy][pos] = tex2D(tex_projes, minh + ((float)pos)/4, proje + projf);

        }

 

        __syncthreads();

 

#pragma unroll 16

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

            int4 minh = {

                cache_subh[2 * tidx][i] - (y) * cache_sin[i],

                cache_subh[2 * tidx + 1][i] - (y) * cache_sin[i],

                cache_subh[2 * tidx][i] - (y + 1) * cache_sin[i],

                cache_subh[2 * tidx + 1][i] - (y + 1) * cache_sin[i],

            };

 

            res[0].x += cache[i][minh.x];

            res[0].y += cache[i][minh.y];

            res[1].x += cache[i][minh.z];

            res[1].y += cache[i][minh.w];

 

        }

 

 

        __syncthreads();

 

    }

 

    *(float2*)(d_SLICE + BLOCK_SIZE_X * 2 * gridDim.x * (idy    ) + idx) = res[0];

    *(float2*)(d_SLICE + BLOCK_SIZE_X * 2 * gridDim.x * (idy + 1) + idx) = res[1];

}

 

 

// Both kernels have problems with pre sm2 builds. Can't tell.

 

#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 = f_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];

        }

    }

}

*/

 

#undef PROJ_LINE

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

 

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

    float res[PPT][PPT];

    

    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[256];

    __shared__ float cache[IPT * PROJ_LINE];

    

    const int tbidy = tidy/(16 / IPT);

    const int tbidx = tidy%(16 / IPT);

    

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

    

    const int sbidx = tidy % 4;

    const int sbidy = tidy / 4;

    

    const int stidx = tidx % 4;

    const int stidy = tidx / 4;

 

    const int sidx = (4 * sbidx + stidx);

    const int sidy = (4 * 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 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 += 8) {

        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();

 

#ifdef HST_OPTIMIZE_KEPLER

        float fminh[8];

 

#pragma unroll 8

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

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

        }

#endif //  HST_OPTIMIZE_KEPLER         

 

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

 

#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 = f_minh[loop_proj + tbidy];

#endif //  HST_OPTIMIZE_KEPLER         

 

 

#pragma unroll 6

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

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

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

                }

 

 

                __syncthreads();

 

#pragma unroll 4 // IPT

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

                    float4 all = 4 * c_all[proje + p];

                    all.y = -all.y;

 

#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 - 4 * minh;

 

                    all.x *= 16.f;

                    all.y *= 16.f;

 

#pragma unroll 4 // PPT

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

#pragma unroll 4 // PPT

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

                            float subh = h + i * all.y + j * all.x;

                            res[i][j] += cache[p*PROJ_LINE + (int)subh];//cache[p * PROJ_LINE + (int)subh];

                        }

                    }

 

 

                }

            }

            __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] = res[i][j];

        }

    }

}