src/lib_kernels_cumuls.cu

   1
   2 texture<int, 2, cudaReadModeElementType> tex_img_in ;
   3
   4 __global__ void calcul_cumuls_gpu(unsigned short * img, t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, uint32 h, uint32 l, uint64 * gsomblocs, unsigned int N_BLOCS_LIGNE, unsigned int baseline, unsigned int nb_lines )
   5 {
   6   int bs = blockDim.x ;
   7   int offsetSom = N_BLOCS_LIGNE * nb_lines;         //indice du premier element de cumul_x2 dans le tableau
   8
   9   int tib = threadIdx.x ;                           //indice du thread dans le bloc
  10   int idx = blockIdx.x*bs + tib ;                   //indice global du thread dans la portion d'image
  11   int num_ligne = blockIdx.x / N_BLOCS_LIGNE ;      //indice de la ligne du thread dans la portion d'image
  12   int til = idx - num_ligne*bs*N_BLOCS_LIGNE ;      //indice du thread ds la ligne
  13   int pos = num_ligne*l + til ;                     //indice du pixel dans la 'tranche' courante
  14   int pos_img = pos + baseline*l ;                  //indice correspondant dans l'image originale
  15   int id, stride=1;
  16   unsigned int z ;
  17
  18   extern __shared__ tcumuls sdata[];
  19   //chargement en smem avec complement de 0
  20   id = CFI(tib);
  21   if (til < l){
  22         z = tex2D(tex_img_in, til, num_ligne) ;
  23             sdata[id].x = z ;
  24             sdata[id].x2 = z*z ;
  25          } else {
  26            sdata[id].x = 0 ;
  27            sdata[id].x2 = 0 ;
  28         }
  29         /*
  30          *prefix sum en smem
  31          */
  32         //passe 1 : construction des sommes
  33
  34         for (int d = bs/2; d > 0; d >>= 1)
  35           {
  36                 __syncthreads();
  37
  38                 if (tib < d)
  39                   {
  40                         int i  = __mul24(__mul24(2, stride), tib);
  41                         int ai = i + stride - 1;
  42                         int bi = CFI(ai + stride);
  43
  44                         ai = CFI(ai);
  45                         //bi = CFI(bi);
  46
  47                         sdata[bi].x += sdata[ai].x;
  48                         sdata[bi].x2 += sdata[ai].x2;
  49
  50                   }
  51                 stride *= 2;
  52           }
  53
  54         //stockage somme du bloc en gmem
  55         // et
  56         //mise a zero dernier element du bloc
  57
  58         __syncthreads();
  59
  60         if (tib == bs-1) {
  61           t_cumul_x som_x = sdata[CFI(tib)].x;
  62           t_cumul_x2 som_x2 = sdata[CFI(tib)].x2;
  63           if (til < l){                                 //il ne faut pas ecrire de somme dans le dernier bloc, elle serait hors image
  64              cumul_x[ pos_img ] = som_x;
  65              cumul_x2[pos_img ] = som_x2;
  66              }
  67           gsomblocs[blockIdx.x] = som_x;
  68           gsomblocs[blockIdx.x + offsetSom] = som_x2;
  69           sdata[CFI(tib)].x = 0;
  70           sdata[CFI(tib)].x2 = 0;
  71         }
  72
  73         //Passe 2 : construction des scans
  74         for (int d = 1; d <= (bs/2); d *= 2)
  75           {
  76                 stride >>= 1;
  77                 __syncthreads();
  78
  79                 if (tib < d)
  80                   {
  81                         int i  = __mul24(__mul24(2, stride), tib);
  82                         int ai = i + stride - 1;
  83                         int bi = CFI(ai + stride);
  84
  85                         ai = CFI(ai);
  86                         //bi = CFI(bi);
  87
  88                         t_cumul_x tx = sdata[ai].x;
  89                         t_cumul_x2 tx2 = sdata[ai].x2;
  90
  91                         sdata[ai].x = sdata[bi].x;
  92                         sdata[bi].x += tx;
  93
  94                         sdata[ai].x2 = sdata[bi].x2;
  95                         sdata[bi].x2 += tx2;
  96                   }
  97           }
  98
  99         //transfert en gmem des scans
 100         //il faut décaler à gauche de 1 element
 101         __syncthreads();
 102         if ( (til < l) && ( tib < (bs-1) ) ){
 103                 cumul_x[ pos_img ] = sdata[CFI(tib+1)].x;
 104                 cumul_x2[pos_img ] = sdata[CFI(tib+1)].x2;
 105          }
 106
 107 }
 108
 109
 110
 111 /*
 112   kernel pour faire le prefixsum des sommes de blocs des lignes de l'image
 113   les params d'executions kivonbien sont :
 114   - threads : la moitié de la nextpow2 du nb de blocs par ligne
 115   - grid : 2 x le nombre de blocs de sommes dans l'image , soit le nb_blocs_par_ligne x nb_lignes x 2
 116   - shared mem : 2 x nextPow2(nb_blocs par ligne)
 117 */
 118
 119
 120
 121 __global__ void scan_somblocs(uint64 * g_idata, int n_bl_l){
 122
 123   extern __shared__  uint64 temp[];
 124
 125   int tib = threadIdx.x;
 126   int toff = blockIdx.x*n_bl_l;         // offset des positions des sommes partielles dans le tableau g_idata
 127   int nmax = blockDim.x*2;              //nb max de sommes partielles par ligne
 128
 129
 130   int aai = tib;
 131   int bbi = tib + blockDim.x;
 132
 133
 134   // chargement data en smem
 135   temp[ CFI(aai) ] = (aai < n_bl_l)? g_idata[ toff + aai ] : 0;
 136   temp[ CFI(bbi) ] = (bbi < n_bl_l)? g_idata[ toff + bbi ] : 0;
 137
 138   int offset = 1;
 139
 140   // passe 1 de bas en haut
 141   for (int d = nmax/2; d > 0; d >>= 1)
 142     {
 143           __syncthreads();
 144
 145           if (tib < d)
 146         {
 147                   int ai = CFI(offset*(2*tib+1)-1);
 148                   int bi = CFI(offset*(2*tib+2)-1);
 149
 150                   temp[bi] += temp[ai];
 151         }
 152
 153           offset *= 2;
 154     }
 155
 156   // zero le dernier element
 157   __syncthreads();
 158   if (tib == 0)
 159     {
 160           int index = CFI(nmax - 1);
 161           g_idata[ toff + n_bl_l - 1 ] =  temp[ index ]; //memorisation de la somme du bloc a sa place en gmem
 162           temp[index] = 0;
 163     }
 164
 165   // passe 2 de haut en bas
 166   for (int d = 1; d < nmax; d *= 2)
 167     {
 168           offset /= 2;
 169
 170           __syncthreads();
 171
 172           if (tib < d)
 173         {
 174                   int ai = CFI(offset*(2*tib+1)-1);
 175                   int bi = CFI(offset*(2*tib+2)-1);
 176
 177                   uint64 t  = temp[ai];
 178                   temp[ai] = temp[bi];
 179                   temp[bi] += t; //tester atomicadd 64 bits
 180                   }
 181
 182     }
 183
 184   __syncthreads();
 185
 186   // transfert resultats en gmem decales d'un element vers la gauche
 187   g_idata[ toff + aai ] = temp[ CFI(aai + 1) ]; //pour le demi tableau smem inferieur
 188   if ( bbi < n_bl_l-1 ) g_idata[ toff + bbi ] = temp[CFI( bbi + 1) ]; //demi tableau smem sup. sauf  dernier=somme.
 189
 190 }
 191
 192
 193 __global__ void add_soms_to_cumuls(t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, uint32 h, uint32 l, uint64 * gsomblocs,
 194                                                                    unsigned int  N_BLOCS_LIGNE, unsigned int baseline, unsigned int nb_lines){
 195
 196   int bs = blockDim.x ;
 197   int offsetSom = N_BLOCS_LIGNE * nb_lines;
 198   int tib = threadIdx.x ;                           //indice du thread dans le bloc
 199   int idx = blockIdx.x*blockDim.x + tib ;           //indice global du thread
 200   int num_ligne = blockIdx.x / N_BLOCS_LIGNE ;      //indice de la ligne du thread dans l'image partielle
 201   int til = idx - num_ligne*bs*N_BLOCS_LIGNE ;      //indice du thread ds la ligne
 202   int pos = num_ligne*l + til;
 203   int pos_img = pos + baseline*l ;
 204
 205
 206   //chargement des valeurs a ajouter en smem 0<-->x et 1<-->x2
 207   uint64 __shared__ ajout[2];
 208
 209   if ( til >= bs ){
 210           ajout[0] = gsomblocs[blockIdx.x -1 ];
 211           ajout[1] = gsomblocs[blockIdx.x -1 + offsetSom];
 212         __syncthreads(); //tester les nouveaux __sync__
 213
 214         //addition par bloc
 215         if (til < l) {
 216           cumul_x[pos_img] += ajout[0];
 217           cumul_x2[pos_img]+= ajout[1];
 218         }
 219   }
 220 }
 221
 222 //calcul des SUM_1, SUM_X et SUM_X2 de l'image
 223
 224 __global__ void calcul_stats_image( t_cumul_x * cumul_x, t_cumul_x2 * cumul_x2, uint32 h, uint32 l, uint64 * d_stats_snake){
 225   uint64 sigX = 0, sigX2 = 0 ;
 226   for (int i=l-1; i<h*l; i+=l)
 227         {
 228           sigX += cumul_x[i] ;
 229           sigX2+= cumul_x2[i];
 230         }
 231   d_stats_snake[3] = h*l ;
 232   d_stats_snake[4] = sigX ;
 233   d_stats_snake[5] = sigX2;
 234
 235 }