src/lib_gpu.cu

   1
   2 #include <stdio.h>
   3 #include <cfloat>
   4
   5 extern "C"{
   6 #include "structures.h"
   7 #include "lib_math.h"
   8 #include "defines.h"
   9 #include "lib_gpu.h"
  10 #include "lib_snake_2_gpu.h"
  11 }
  12 #include "lib_test_gpu.h"
  13 #include "lib_kernels_cumuls.cu"
  14 #include "lib_kernel_snake_2_gpu.cu"
  15
  16 #define DEBUG_IMG_CUMUL 1
  17 bool DISPLAY_ERR_IMG_CUMUL = 1;
  18 //#define DEBUG_POSITIONS
  19 //#define DEBUG_MOVE
  20 //#define DEBUG_CRST
  21 //#define DEBUG_MV
  22 //#define DEBUG_SOMSOM
  23 //#define DEBUG_SOMBLOCS
  24 //#define DEBUG_LISTES
  25 //#define DEBUG_STATS_REF
  26
  27 double critere_gauss_cpu( int n, int ni, uint64 sxi, uint64 sxi2, int64 *stats_img){
  28   int ne = n - ni ;
  29   uint64 sxe = stats_img[4] - sxi ;
  30   uint64 sxe2= stats_img[5] - sxi2;
  31   double sigi2, sige2, critere ;
  32
  33   /* variance des valeurs des niveaux de gris a l'interieur */
  34   sigi2 =
  35         ((double)sxi2/ni) -
  36         ((double)sxi/ni)*((double)sxi/ni) ;
  37
  38   /* variance des valeurs des niveaux de gris a l'exterieur */
  39   sige2 =
  40         ((double)sxe2)/ne -
  41         ((double)sxe/ne)*((double)sxe/ne) ;
  42
  43   if ( (sigi2 > 0) && (sige2 > 0) )
  44         critere =  0.5*((double)ni*log(sigi2) + (double)ne*log(sige2)) ;
  45   else critere = DBL_MAX ;
  46   return critere;
  47 }
  48
  49 void cuda_init_img_cumul(unsigned short ** img_in, int H, int L, int nb_nodes,
  50                                                  unsigned short ** d_img, t_cumul_x ** d_img_x, t_cumul_x2 ** d_img_x2,
  51                                                  int ** d_freemanDiDj, int ** d_codeNoeud,
  52                                                  snake_node_gpu ** d_snake, uint32 ** d_nb_pix_max,
  53                                                  uint4 ** d_positions, uint64 ** d_contribs_segments, uint4 ** d_freemans_centres,
  54                                                  int ** d_codes_segments, int64 ** d_stats_snake,
  55                                                  int64 ** d_stats, int64 ** d_stats_ref, double ** d_vrais, double ** d_vrais_snake,
  56                                                  uint2 ** d_liste_pixels, uint64 ** d_contribs_segments_blocs,
  57                                                  bool ** d_move
  58                                                  )
  59 {
  60   unsigned int taille = H*L;
  61   timeval chrono;
  62
  63
  64   //allocation cumuls en memoire GPU
  65   tic(&chrono, NULL);
  66   /*
  67         MAX_PIX 20000
  68         MAX_NODES 10000
  69         MAX_LISTE_PIX 10000000
  70    */
  71   cudaMalloc( (void**) d_snake, MAX_NODES*sizeof(snake_node_gpu) );
  72
  73   cudaMalloc( (void**) d_img, taille*sizeof(unsigned short) );
  74   cudaMalloc( (void**) d_img_x, taille*sizeof(t_cumul_x) );
  75   cudaMalloc( (void**) d_img_x2, taille*sizeof(t_cumul_x2) );
  76
  77   cudaMalloc( (void**) d_freemanDiDj, 9*sizeof(int) );
  78   cudaMalloc( (void**) d_codeNoeud, 64*sizeof(int) );
  79
  80   cudaMalloc( (void**) d_stats_snake, 6*sizeof(int64)) ;
  81   cudaMalloc( (void**) d_positions, 8*MAX_NODES*sizeof(uint4)) ;
  82   cudaMalloc( (void**) d_contribs_segments, 3*16*MAX_NODES*sizeof(uint64)) ;
  83   cudaMalloc( (void**) d_contribs_segments_blocs, (3*MAX_LISTE_PIX/32)*sizeof(uint64)) ;
  84   cudaMalloc( (void**) d_freemans_centres, 16*MAX_NODES*sizeof(uint4)) ;
  85   cudaMalloc( (void**) d_codes_segments, 16*MAX_NODES*sizeof(int)) ;
  86   cudaMalloc( (void**) d_stats, 3*8*MAX_NODES*sizeof(int64)) ;
  87   cudaMalloc( (void**) d_stats_ref, 3*MAX_NODES*sizeof(int64)) ;
  88   cudaMalloc( (void**) d_vrais, 8*MAX_NODES*sizeof(double)) ;
  89   cudaMalloc( (void**) d_move, MAX_NODES*sizeof(bool)) ;
  90   cudaMalloc( (void**) d_nb_pix_max, sizeof(uint32)) ;
  91   cudaMalloc( (void**) d_vrais_snake, sizeof(double)) ;
  92
  93   cudaMalloc( (void**) d_liste_pixels, 16*5*(MAX_NODES)*sizeof(uint2) );
  94
  95   printf("TOTAL MEM = %ld octets\n",
  96                  (2*MAX_NODES*(sizeof(snake_node_gpu)+(8+16)*sizeof(uint4)+3*16*8+16*4+24*8+3*8+8*sizeof(double)+sizeof(bool))
  97                  +(MAX_LISTE_PIX)*(sizeof(uint2)+1)
  98                  +taille*(8+sizeof(t_cumul_x)+sizeof(t_cumul_x2))
  99                   +9*4+64*4+6*8+4+sizeof(double)) );
 100
 101   int64 * h_stats_snake = new int64[6];
 102
 103   toc(chrono, "temps alloc mem GPU");
 104
 105   /*detection-choix-initialisation de la carte GPU*/
 106   tic(&chrono, NULL) ;
 107   cudaDeviceProp deviceProp;
 108   deviceProp.major = 2;
 109   deviceProp.minor = 0;
 110   int dev;
 111   cudaChooseDevice(&dev, &deviceProp);
 112   cudaGetDeviceProperties(&deviceProp, dev);
 113   if(deviceProp.major >= 2 )
 114         {
 115           printf("Using Device %d: \"%s\"\n", dev, deviceProp.name);
 116           cudaSetDevice(dev);
 117         }
 118   toc(chrono, "temps acces GPU") ;
 119
 120   //copie tables correspondances freeman en mem GPU
 121   tic(&chrono, NULL) ;
 122   cudaMemcpy( *d_freemanDiDj, CORRESPONDANCE_Di_Dj_FREEMAN , 9*sizeof(int), cudaMemcpyHostToDevice);
 123   cudaMemcpy( *d_codeNoeud, TABLE_CODAGE , 64*sizeof(unsigned int), cudaMemcpyHostToDevice);
 124   toc(chrono, "temps transfert tables de codage") ;
 125
 126   /*transfert image en global mem GPU*/
 127   tic(&chrono, NULL);
 128   cudaMemcpy( *d_img, img_in[0], taille*sizeof(unsigned short), cudaMemcpyHostToDevice);
 129   toc(chrono, "transfert image vers GPU");
 130
 131   //calculs images cumulees sur GPU
 132   int blocs_max = 65536 ;
 133   int bs = 256 ; //arbitraire, d'apres les observations c'est souvent l'optimu
 134   unsigned int base = 0 ;
 135   unsigned int bl_l  = (L+bs-1)/bs ;
 136   unsigned int nb_lines =  blocs_max / bl_l ;
 137   unsigned int lines ;
 138   unsigned int tranches = ( 1 + H / nb_lines ) ;
 139   nb_lines = (H +tranches -1)/ tranches ; // equilibre la taille des tranches
 140
 141   dim3 threads(bs,1,1);
 142   int smem = nextPow2(bl_l)*2; //smem pour le prefixscan des sommes de blocs (etape 2)
 143   smem += smem >> DEC;
 144   smem += smem >> DEC;
 145   int smem_size = smem*sizeof(uint64);
 146   uint64 * d_somblocs ; // sommes des cumuls par bloc de calcul
 147
 148
 149   if(DEBUG_IMG_CUMUL)
 150         {
 151           printf("--- CALCULS IMAGES CUMULEES+STATS GPU  ----\n");
 152           printf("\t%d threads par bloc  -- %u blocs par ligne -- %u tranches -- %u lignes par tranche \n",bs, bl_l, tranches,nb_lines);
 153           printf(" Smem totale pour cumuls : %d\n", CFI(bs)*(sizeof(t_cumul_x)+sizeof(t_cumul_x2)) );
 154           tic(&chrono, NULL);
 155         }
 156   //calculs cumuls generiques : necessitent 3 etapes / 3 kernels
 157   cudaMalloc( (void**) &d_somblocs, 2*bl_l*nb_lines*sizeof(uint64) );
 158   cudaFuncSetCacheConfig(calcul_cumuls_gpu, cudaFuncCachePreferShared);
 159   do
 160         {
 161           if  ( H-base < nb_lines ) lines = H - base ; else lines = nb_lines ;
 162           printf("base = ligne %d -- traitement de %d lignes \n", base, lines) ;
 163           dim3 grid(bl_l*lines,1,1) ;
 164           calcul_cumuls_gpu<<<grid, threads, CFI(bs)*sizeof(tcumuls)>>>(*d_img, *d_img_x, *d_img_x2, H, L, d_somblocs, bl_l, base, lines) ;
 165           scan_somblocs<<<2*lines, nextPow2(bl_l)/2, smem_size>>>(d_somblocs, bl_l) ;
 166           add_soms_to_cumuls<<<grid, threads>>>(*d_img_x, *d_img_x2, H, L, d_somblocs, bl_l, base, lines) ;
 167           base += lines ;
 168         }
 169   while (base < H) ;
 170   cudaFree(d_somblocs) ;
 171
 172   //calcul des sommes totales N, sigX et sigX2 sur l'image
 173   calcul_stats_image<<<1, 1>>>( *d_img_x, *d_img_x2, H, L, (uint64*)*d_stats_snake);
 174
 175   cudaThreadSynchronize()   ;
 176   toc(chrono, "\tTemps GPU");
 177   //allocation memoire CPU
 178           t_cumul_x  * img_x = new t_cumul_x [H*L];
 179           t_cumul_x2 *  img_x2 = new t_cumul_x2 [H*L];
 180           uint64 sigX = 0, sigX2 = 0 ;
 181   if(DEBUG_IMG_CUMUL)
 182         {
 183
 184           /*pour test comparaison*/
 185           t_cumul_x * img_xb = new t_cumul_x [H*L];
 186           t_cumul_x2 * img_x2b = new t_cumul_x2 [H*L];
 187
 188           cudaMemcpy( img_xb, *d_img_x, taille*sizeof(t_cumul_x), cudaMemcpyDeviceToHost);
 189           cudaMemcpy( img_x2b, *d_img_x2, taille*sizeof(t_cumul_x2), cudaMemcpyDeviceToHost);
 190
 191           //cumuls : etape 1 CPU
 192           /*
 193                 for (int i=0; i<H; i++)
 194                 {
 195                         for (int b=0; b<bl_l; b++)
 196                         {
 197                                 int offset = b*bs ;
 198                                 img_x[i*L+offset] = img_in[i][offset] ;
 199                                 img_x2[i*L+offset]= img_in[i][offset]*img_in[i][offset] ;
 200                                 for (int p=1; p<bs; p++)
 201                                 {
 202                                         int j = p+offset ;
 203                                         if (j<L)
 204                                         {
 205                                                 img_x[i*L+j] = img_x[i*L+j-1] + img_in[i][j];
 206                                                 img_x2[i*L+j] = img_x2[i*L+j-1] + img_in[i][j]*img_in[i][j] ;
 207                                         }
 208                                 }
 209                         }
 210                 }
 211           */
 212           //cumuls complets CPU
 213
 214           for (int i=0; i<H; i++)
 215                 {
 216                   img_x[i*L+0] = img_in[i][0] ;
 217                   img_x2[i*L+0]= img_in[i][0]*img_in[i][0] ;
 218                   for (int j=1; j<L; j++)
 219                         {
 220                           img_x[i*L+j]  = img_x[i*L+j-1]  + img_in[i][j] ;
 221                           img_x2[i*L+j] = img_x2[i*L+j-1] + img_in[i][j]*img_in[i][j] ;
 222                         }
 223                 }
 224
 225           int cpt = 0;
 226           int cpt_errx=0, cpt_errx2 = 0;
 227           for (int i=0; i< H; i++){
 228                 for (int j=0; j< L; j++){
 229                   if ( (img_x[i*L+j] !=  img_xb[i*L+j]) ) cpt_errx++ ;
 230                   if ( (img_x2[i*L+j] !=  img_x2b[i*L+j]) ) cpt_errx2++ ;
 231                   if ( (img_x[i*L+j] !=  img_xb[i*L+j]) || (img_x2[i*L+j] !=  img_x2b[i*L+j]))
 232                   {
 233                         //printf("(%d,%d)sxCPU=%lu  sxGPU=%lu -- sx2CPU=%lu  sx2GPU=%lu\n",i,j,img_x[i*L+j], img_xb[i*L+j], img_x2[i*L+j], img_x2b[i*L+j]);
 234                   }
 235                   cpt++;
 236                 }
 237           }
 238           printf("%d erreurs sur CX / %d points\n", cpt_errx, cpt );
 239           printf("%d erreurs sur CX2 / %d points\n", cpt_errx2, cpt );
 240
 241           for (int i=0; i<H; i++)
 242                 {
 243                   sigX += img_x[i*L+L-1] ;
 244                   sigX2+= img_x2[i*L+L-1];
 245                 }
 246           printf("STATS IMAGE  N = %d - sigX = %lu - sigX2 = %lu\n",  H*L, sigX, sigX2 );
 247         }
 248
 249   /*
 250    * generation snake en mem GPU
 251    */
 252   int dist = 140 ;
 253
 254   /* Test de determination de la fenetre verticale  optimale*/
 255   int div = 256;//nb de divisions de l'image : cela définit le pas. La valeur max découle du nb max de threads possible ds une grille
 256   tcontribs *d_contribs_part ; // pour les contribs partielles par bloc
 257   tcontribs *d_contribs_cols, *h_contribs_cols, *h_contribs_cols_cpu ; // pour les contrib des colonnes
 258   double2 *d_miniblocs, *h_miniblocs ;       // pour les minima de chaque colonne haute
 259
 260   //parametres execution
 261   bs = div ;
 262   int bpc= (H + bs -1)/bs ;
 263   dim3 grid = dim3(div, bpc, 1);
 264   smem = CFI(bs)*sizeof(tcontribs) ;
 265
 266   //allocations
 267   cudaMalloc((void**) &d_contribs_part, div*bpc*sizeof(tcontribs));
 268   cudaMalloc((void**) &d_contribs_cols, div*sizeof(tcontribs)) ;
 269   cudaMalloc((void**) &d_miniblocs, div*sizeof(double2)) ;
 270   h_contribs_cols = new tcontribs[div] ;
 271   h_contribs_cols_cpu = new tcontribs[div] ;
 272   h_miniblocs = new double2[div] ;
 273
 274   int pas = L / div ;
 275
 276   tic(&chrono, NULL);
 277
 278   //execution kernels
 279   calcul_contribs_cols<<<grid, bs, smem>>>( *d_img_x, *d_img_x2, H, L, d_contribs_part );
 280   somsom_contribs<<<div,1>>>( d_contribs_part, bpc, d_contribs_cols ) ;
 281   calcul_contribs_permutations<<<div,div, CFI(div)*sizeof(double2)>>>(d_contribs_cols, bpc, H, L, (uint64*) *d_stats_snake, d_miniblocs) ;
 282   cudaMemcpy( h_miniblocs, d_miniblocs, div*sizeof(double2), cudaMemcpyDeviceToHost ) ;
 283
 284   //verif minimum des blocs
 285   double crit_mini = h_miniblocs[0].x;
 286   int j1=0, id_mini=0;
 287   for (j1=1 ; j1 < div ; j1++){
 288         if (h_miniblocs[j1].x < crit_mini) {
 289           crit_mini = h_miniblocs[j1].x ;
 290           id_mini = j1 ;
 291         }
 292   }
 293   toc(chrono, "\nCALCUL RECTANGLE");
 294   j1 = pas * id_mini ;
 295   int j2 = (int)(pas*h_miniblocs[ id_mini ].y) ;
 296   printf("pas = %d cols -- critere mini =%f , positions j1=%d j2=%d\n", pas, h_miniblocs[ id_mini ].x, j1, j2);
 297
 298
 299   // transfert datas GPU -> CPU
 300   cudaMemcpy( h_contribs_cols, d_contribs_cols, div*sizeof(tcontribs), cudaMemcpyDeviceToHost ) ;
 301
 302   //verif contribs colonnes
 303   for (int c=0 ; c < div ; c++){
 304         // calcul valeurs de ref en CPU
 305         h_contribs_cols_cpu[c].cx = 0 ;
 306         h_contribs_cols_cpu[c].cx2= 0 ;
 307         for (int ip=0; ip < H; ip++){
 308           h_contribs_cols_cpu[c].cx  += img_x[ ip*L + c*pas] ;
 309           h_contribs_cols_cpu[c].cx2 += img_x2[ ip*L + c*pas] ;
 310         }
 311         //comparaison avec valeurs GPU
 312         if ( (h_contribs_cols_cpu[c].cx != h_contribs_cols[c].cx) || (h_contribs_cols_cpu[c].cx2 != h_contribs_cols[c].cx2) )
 313           printf("ERR colonne %d -> CPUx=%lu CPUx2=%lu | GPUx=%lu GPUx2=%lu\n",
 314                          c*pas, h_contribs_cols_cpu[c].cx, h_contribs_cols_cpu[c].cx2, h_contribs_cols[c].cx, h_contribs_cols[c].cx2 );
 315   }
 316   cudaFree(d_contribs_part) ;
 317   cudaFree(d_contribs_cols) ;
 318   cudaFree(d_miniblocs) ;
 319   free(h_contribs_cols);
 320   free(h_contribs_cols_cpu);
 321   free(h_miniblocs) ;
 322
 323   //realloc pour lignes horizontales
 324   bs = 128 ;
 325
 326   div = (H+bs-1)/bs ;
 327   printf("DIV = %d\n", div ) ;
 328
 329   int divpow2 = nextPow2(div) ;
 330   printf("DIVPOW2 = %d\n", divpow2) ;
 331
 332   grid = dim3(div, 1, 1) ;
 333   smem = CFI(bs)*sizeof(tcontribs) ;
 334   cudaMalloc((void**) &d_contribs_part, div*sizeof(tcontribs)) ;
 335   cudaMalloc((void**) &d_contribs_cols, div*div*sizeof(tcontribs)) ;
 336   cudaMalloc((void**) &d_miniblocs, div*sizeof(double2)) ;
 337   h_contribs_cols = new tcontribs[div*div] ;
 338   tcontribs * h_contribs_part = new tcontribs[div] ;
 339   h_miniblocs = new double2[div] ;
 340
 341   tic(&chrono, NULL);
 342   // Appels kernels optim lignes horizontales
 343   calcul_contrib_conjuguee_colonnes<<<grid, bs, CFI(bs)*sizeof(tcontribs) >>>( *d_img_x, *d_img_x2, H, L, j1, j2, d_contribs_part) ;
 344
 345   /*verif CPU
 346   int cpt = 0 ;
 347   int cpterr = 0 ;
 348   tcontribs * h_contribs_part_cpu = new tcontribs[div] ;
 349   cudaMemcpy( h_contribs_part, d_contribs_part, div*sizeof(tcontribs), cudaMemcpyDeviceToHost ) ;
 350   for (int bloc=0; bloc < div; bloc++){
 351         h_contribs_part_cpu[ bloc ].cx = 0 ;
 352         h_contribs_part_cpu[ bloc ].cx2 = 0 ;
 353           for (int line=0; ((line < bs)&&(bloc*bs+line < H)); line++){
 354                 h_contribs_part_cpu[bloc].cx += img_x[ (bloc*bs+line)*L + j2] - img_x[ (bloc*bs+line)*L + j1 ];
 355                 h_contribs_part_cpu[bloc].cx2 += img_x2[ (bloc*bs+line)*L + j2] - img_x2[ (bloc*bs+line)*L + j1 ];
 356           }
 357           if ( ( h_contribs_part_cpu[bloc].cx != h_contribs_part[bloc].cx ) || ( h_contribs_part_cpu[bloc].cx2 != h_contribs_part[bloc].cx2 ) )
 358                 {
 359                   printf("ERREUR bloc %d -> CPUx=%lu CPUx2=%lu | GPUx=%lu GPUx2=%lu\n", bloc,
 360                                  h_contribs_part_cpu[bloc].cx, h_contribs_part_cpu[bloc].cx2, h_contribs_part[bloc].cx, h_contribs_part[bloc].cx2 ) ;
 361                   cpterr++;
 362                 }
 363           cpt++ ;
 364   }
 365   printf("VERIF CONTRIB CONJUGUEES BLOCS -->  %d ERREURS / %d BLOCS\n", cpterr, cpt) ;
 366   fin verif*/
 367
 368   grid = dim3(div, div, 1) ;
 369   calcul_contribs_snake_rectangle<<<grid,divpow2, CFI(divpow2)*sizeof(tcontribs) >>>(d_contribs_part, d_contribs_cols) ;
 370
 371   /* verif CPU
 372   h_contribs_cols_cpu = new tcontribs[div*div] ;
 373   cudaMemcpy( h_contribs_cols, d_contribs_cols, div*div*sizeof(tcontribs), cudaMemcpyDeviceToHost ) ;
 374   cpt = 0 ;
 375   cpterr = 0 ;
 376   for (int i1=0; i1 < div ; i1++){
 377         for (int i2=0 ; i2 < div ; i2++){
 378           if (i2 >= i1){
 379                 h_contribs_cols_cpu[ i1*div +i2 ].cx = 0 ;
 380                 h_contribs_cols_cpu[ i1*div +i2 ].cx2= 0 ;
 381                 for (int d=i1 ; d <= i2 ; d++){
 382                   h_contribs_cols_cpu[ i1*div +i2 ].cx += h_contribs_part_cpu[ d ].cx ;
 383                   h_contribs_cols_cpu[ i1*div +i2 ].cx2+= h_contribs_part_cpu[ d ].cx2 ;
 384                 }
 385           } else {
 386                 h_contribs_cols_cpu[ i1*div +i2 ].cx = 0 ;
 387                 h_contribs_cols_cpu[ i1*div +i2 ].cx2= 0 ;
 388           }
 389
 390           if (( ( h_contribs_cols_cpu[ i1*div +i2 ].cx != h_contribs_cols[ i1*div +i2 ].cx ) || ( h_contribs_cols_cpu[ i1*div +i2].cx2 != h_contribs_cols[ i1*div +i2].cx2 ) )
 391                 && (i2 >= i1))
 392                 {
 393                   printf("ERREUR combinaison (%d, %d) -> CPUx=%lu CPUx2=%lu | GPUx=%lu GPUx2=%lu\n", i1, i2,
 394                                  h_contribs_cols_cpu[ i1*div +i2].cx, h_contribs_cols_cpu[ i1*div +i2 ].cx2, h_contribs_cols[ i1*div +i2 ].cx, h_contribs_cols[ i1*div +i2 ].cx2 ) ;
 395                   cpterr++;
 396                 }
 397           cpt++ ;
 398         }
 399   }
 400   printf("VERIF COMBINAISONS LIGNES -->  %d ERREURS / %d COMBINAISONS\n", cpterr, cpt) ;
 401   fin verif */
 402
 403
 404   grid = dim3(div, 1, 1) ;
 405   calcul_critere_permutations_verticales<<< grid, divpow2, CFI(divpow2)*sizeof(double2) >>>(d_contribs_cols, bs, j1, j2, H, L, sigX, sigX2, d_miniblocs) ;
 406
 407   /* verif CPU
 408   cpt = 0 ;
 409   cpterr = 0 ;
 410   double2 * h_miniblocs_cpu = new double2[ div ] ;
 411   cudaMemcpy( h_miniblocs, d_miniblocs, div*sizeof(double2), cudaMemcpyDeviceToHost) ;
 412   cudaMemcpy( h_stats_snake, *d_stats_snake, 6*sizeof(int64), cudaMemcpyDeviceToHost) ;
 413   for (int lb=0 ; lb < div ; lb++){
 414         h_miniblocs_cpu[lb].x = DBL_MAX ;
 415         for (int lh=lb ; lh < div ; lh++){
 416           if ( critere_gauss_cpu(H*L, (lh-lb+1)*bs*(j2 - j1), h_contribs_cols_cpu[ lb*div +lh ].cx, h_contribs_cols_cpu[ lb*div + lh ].cx2, h_stats_snake ) < h_miniblocs_cpu[lb].x )
 417                 {
 418                   h_miniblocs_cpu[lb].x = critere_gauss_cpu(H*L, (lh-lb+1)*bs*(j2 - j1), h_contribs_cols_cpu[ lb*div +lh ].cx, h_contribs_cols_cpu[ lb*div + lh ].cx2, h_stats_snake) ;
 419                   h_miniblocs_cpu[lb].y = (double)lh ;
 420                 }
 421         }
 422         if ( ( h_miniblocs_cpu[lb].x > 1.000001*h_miniblocs[lb].x ) || ( h_miniblocs_cpu[lb].x < 0.999999*h_miniblocs[lb].x ) || ( h_miniblocs_cpu[lb].y != h_miniblocs[lb].y ) )
 423           {
 424                  printf("ERREUR MINIMUM BLOC LIGNE %d -> CPU=%lf en i2=%d | GPU=%lf en i2=%d\n", lb,
 425                                 h_miniblocs_cpu[ lb ].x, (int)h_miniblocs_cpu[ lb ].y, h_miniblocs[ lb ].x, (int)h_miniblocs[ lb ].y ) ;
 426                   cpterr++;
 427           }
 428         cpt++ ;
 429   }
 430   printf("VERIF MINIMA PAR BLOC -->  %d ERREURS / %d BLOCS\n", cpterr, cpt) ;
 431   */
 432   /* fin verif */
 433
 434   /*
 435    * determination  minimum absolu
 436    * a conserver sur CPU
 437    */
 438   cudaMemcpy( h_miniblocs, d_miniblocs, div*sizeof(double2), cudaMemcpyDeviceToHost) ;
 439   crit_mini = h_miniblocs[0].x ;
 440   int i1=0  ;
 441   id_mini=0 ;
 442   for (i1=1 ; i1 < div ; i1++){
 443         if (h_miniblocs[i1].x < crit_mini) {
 444           crit_mini = h_miniblocs[i1].x ;
 445           id_mini = i1 ;
 446         }
 447   }
 448
 449   i1 = bs * id_mini ;
 450   int i2 = (int)(bs*h_miniblocs[ id_mini ].y) ;
 451
 452   toc(chrono, "CALCUL RECTANGLE");
 453
 454   printf("pas = %d lignes -- critere mini =%f , positions i1=%d i2=%d\n", bs, h_miniblocs[ id_mini ].x, i1, i2);
 455   /*fin test snake rectangle initial optimal*/
 456
 457   tic(&chrono, NULL);
 458   //genere_snake_rectangle_4nodes_gpu<<< 1, 1>>>(*d_snake, 140, H, L) ;
 459   //genere_snake_bande_gpu<<<1,1>>>(*d_snake, pas*id_mini, (int)(pas*h_miniblocs[ id_mini ].y), H);
 460   genere_snake_rectangle<<<1,1>>>(*d_snake, i1, i2, j1, j2);
 461
 462
 463   int nnodes = nb_nodes ;
 464   snake_node_gpu * h_snake = new snake_node_gpu[nnodes];
 465   snake_node * h_snake_ll = new snake_node[nnodes] ;
 466   uint4 * h_liste_positions = new uint4[nnodes*8];
 467   double * h_vrais_snake = new double ;
 468   //init les stats du snake
 469   uint2 * d_liste_temp  ;
 470   t_sum_x2 * d_sompart  ;
 471   int tpb, bps, npixmax ;
 472
 473   //calcul nb threads par bloc
 474   npixmax = 2*(H+L-4*dist)/(nnodes-1) ;
 475   tpb = nextPow2(npixmax) ;
 476   if (tpb >= 256) tpb = 256 ;//  /!\ le kernel <<< calcul_contrib...>>> ne supporte pas un bs>256 a cause de la shared-mem nécessaire
 477   if (tpb < 32 ) tpb = 32 ;
 478   tpb=128 ;
 479   bps = (npixmax+tpb-1)/tpb ;
 480   printf("PARAMS EXEC INIT : %d pix max, %d threads/bloc, %d blocs/seg, %d blocs/grille\n", npixmax, tpb, bps, nnodes*bps);
 481
 482   //alloc
 483   cudaMalloc((void**) &d_liste_temp, nnodes*bps*tpb*sizeof(uint2));
 484   cudaMalloc((void**) &d_sompart, 3*nnodes*bps*sizeof(t_sum_x2));
 485   cudaMalloc((void**) &d_stats_ref, 3*nnodes*sizeof(int64));
 486
 487   //DEBUG : pour forcer la mise à zero du tableau intermediaire d_stats_ref
 488   /*
 489   int64 h_stats_ref[3*nnodes] ;
 490   for (int a=0; a<3*nnodes ; a++) h_stats_ref[a] = 0 ;
 491   cudaMemcpy( h_stats_ref, d_stats_ref, sizeof(int64), cudaMemcpyHostToDevice) ;
 492   */
 493   //fin forçage a 0
 494
 495   //DEBUG : pour forcer la mise à zero du tableau intermediaire d_sompart
 496   /*
 497      t_sum_x2 h_sompart[ 3*nnodes*bps ] ;
 498      for (int a=0; a<3*nnodes*bps ; a++) h_sompart[a] = 0 ;
 499      cudaMemcpy( h_sompart, d_sompart, sizeof(t_sum_x2), cudaMemcpyHostToDevice) ;
 500   */
 501   //fin forçage a 0
 502
 503   calcul_contribs_segments_snake<<< nnodes*bps, tpb, (CFI(tpb))*(3*sizeof(t_sum_x2))>>>
 504         (*d_snake, nnodes,
 505          *d_img_x, *d_img_x2,
 506          L, d_liste_temp, d_sompart, *d_freemanDiDj );
 507
 508   //TODO
 509   //parametrer pour ne pas appeler qd tpb=1
 510   //oblige a modifier le kernel <<< calcul_contrib...>>> pour ecrire directement ds d_snake
 511   // au lieu de d_sompart
 512   somsom_snake<<< nnodes , 1 >>>(d_sompart, nnodes, bps, *d_snake);
 513
 514
 515   calcul_stats_snake<<< 1 , 1 >>>(*d_snake, nnodes, *d_stats_snake, *d_vrais_snake,
 516                                                                   *d_img_x, *d_img_x2,
 517                                                                   *d_codeNoeud, L
 518                                                                   );
 519   cudaThreadSynchronize() ;
 520   toc(chrono, "\tTemps") ;
 521
 522   /*
 523         verif stats initiales du snake
 524   */
 525   cudaMemcpy( h_vrais_snake, *d_vrais_snake, sizeof(double), cudaMemcpyDeviceToHost) ;
 526   cudaMemcpy( h_stats_snake, *d_stats_snake, 6*sizeof(int64), cudaMemcpyDeviceToHost) ;
 527
 528   printf("STATS SNAKE log vrais=%lf : c1=%lu - cx=%lu - cx2=%lu - N=%lu - SUMX=%lu - SUMX2=%lu\n",
 529                  *h_vrais_snake,
 530                  h_stats_snake[0],  h_stats_snake[1],  h_stats_snake[2],
 531                  h_stats_snake[3],  h_stats_snake[4],  h_stats_snake[5] );
 532
 533   /*
 534         verif stats diminuees des contribs des 2 segments associes a chaque noeud
 535   */
 536 #ifdef DEBUG_STATS_REF
 537   cudaMemcpy( h_stats_ref, d_stats_ref, 3*nnodes*sizeof(int64), cudaMemcpyDeviceToHost) ;
 538   cudaMemcpy( h_snake, *d_snake, nnodes*sizeof(snake_node_gpu), cudaMemcpyDeviceToHost) ;
 539
 540
 541   printf("******* STATS DIMINUEES\n");
 542   for(int n=0; n<nnodes;n++)
 543         {
 544           int i = h_snake[n].posi, j = h_snake[n].posj ;
 545           printf("node %d (%d,%d) : %ld - %ld - %ld - img1= %lu - imgx= %lu - imgx2= %lu \n", n, i, j,
 546                          h_stats_ref[3*n], h_stats_ref[3*n +1], h_stats_ref[3*n +2],
 547                          img_1[i][j], img_x[i][j], img_x2[i][j]);
 548         }
 549 #endif //DEBUG_STATS_REF
 550
 551   //snake2gpu(d_snake, snake, nb_nodes);
 552   //gpu2snake(*d_snake, &h_snake_ll, nnodes);
 553
 554
 555 #ifdef DEBUG_POSITIONS
 556   for (int n=0; n<nnodes; n++)
 557         {
 558           printf("Node %d :\n", n);
 559           for (int pos=0; pos<8; pos++)
 560                 {
 561                   printf("(%d , %d):%d:%d | ", h_liste_positions[8*n + pos].x, h_liste_positions[8*n + pos].y,
 562                                  h_liste_positions[8*n + pos].z, h_liste_positions[8*n + pos].w);
 563                 }
 564           printf("\n");
 565         }
 566 #endif //DEBUG_POSITIONS
 567
 568 //verif liste pixels noeuds pairs/impairs selon
 569
 570 #ifdef DEBUG_LISTES
 571   printf("NOMBRE PIXELS pour LISTE = %d\n", *h_nb_pix_max) ;
 572   printf("bs = %d - grid = %d - nblocs_seg = %d - npix_max = %d - taille mem = %d\n",
 573                  bs, grid.x, nblocs_seg, *h_nb_pix_max, taille_mem);
 574
 575   cudaMemcpy( h_liste_pix, d_liste_pix, taille_mem*sizeof(uint2), cudaMemcpyDeviceToHost ) ;
 576   cudaMemcpy( h_snake, *d_snake, nnodes*sizeof(snake_node_gpu), cudaMemcpyDeviceToHost );
 577   uint32 * h_liste_pixels_segment = new uint32[2*(*h_nb_pix_max)] ;
 578   int idx_n, idx_nprec, idx_nsuiv ;
 579
 580   printf("********* LISTE PIX  ***********\n");
 581   printf("bs = %d - grid = %d - nblocs_seg = %d - npix_max = %d - taille mem = %d\n",
 582                  bs, grid.x, nblocs_seg, *h_nb_pix_max, taille_mem);
 583
 584   for (int n=0; n<(nnodes/2 + (nnodes%2)*pairs); n++)
 585         {
 586           idx_n = 2*n + !pairs ;
 587           if (idx_n == 0) idx_nprec = nnodes - 1;
 588           else idx_nprec = idx_n - 1;
 589           if (idx_n == nnodes-1) idx_nsuiv = 0;
 590           else idx_nsuiv = idx_n + 1 ;
 591
 592           for (int pos=0; pos < 8 ; pos++) //test des segments avant le noeud
 593                 {
 594
 595                   int nb_pix = calcul_liste_pixel_segment(h_snake[idx_nprec].posi,h_snake[idx_nprec].posj,
 596                                                                                                   h_liste_positions[8*idx_n+pos].x, h_liste_positions[8*idx_n+pos].y,
 597                                                                                                   h_liste_pixels_segment, 0);
 598
 599                   for (int pix=0; pix < nb_pix; pix++)
 600                         {
 601
 602                           if ( (h_liste_pix[(16*n + pos)*nblocs_seg*bs + pix].x != h_liste_pixels_segment[2*pix] )
 603                                    || ( h_liste_pix[(16*n + pos)*nblocs_seg*bs + pix].y != h_liste_pixels_segment[2*pix+1] ) )
 604                                 printf("erreur avant n=%d pix %d/%d segment %d noeuds[ %d-%d-%d ] , CPU (%d,%d) - GPU (%d, %d)\n", n, pix, nb_pix, pos,
 605                                            idx_nprec, idx_n, idx_nsuiv,
 606                                            h_liste_pixels_segment[2*pix], h_liste_pixels_segment[2*pix+1],
 607                                            h_liste_pix[(16*n + pos)*nblocs_seg*bs + pix].x,  h_liste_pix[(16*n + pos)*nblocs_seg*bs + pix].y);
 608
 609                         }
 610
 611                 }
 612           for (int pos=0; pos < 8 ; pos++) //test des segments apres le noeud
 613                 {
 614
 615                   int nb_pix = calcul_liste_pixel_segment(h_liste_positions[8*idx_n+pos].x, h_liste_positions[8*idx_n+pos].y,
 616                                                                                                   h_snake[idx_nsuiv].posi,h_snake[idx_nsuiv].posj,
 617                                                                                                   h_liste_pixels_segment, 0);
 618
 619                   for (int pix=0; pix < nb_pix; pix++)
 620                         {
 621
 622                           if ( (h_liste_pix[(16*n + pos + 8)*nblocs_seg*bs + pix].x != h_liste_pixels_segment[2*pix] )
 623                              || ( h_liste_pix[(16*n + pos + 8)*nblocs_seg*bs + pix].y != h_liste_pixels_segment[2*pix+1] ) )
 624                                 printf("erreur apres n=%d pix %d/%d segment %d noeuds[ %d-%d-%d ] , CPU (%d,%d) - GPU (%d, %d)\n", n, pix, nb_pix, pos+8,
 625                                            idx_nprec, idx_n, idx_nsuiv,
 626                                            h_liste_pixels_segment[2*pix], h_liste_pixels_segment[2*pix+1],
 627                                            h_liste_pix[(16*n + pos + 8)*nblocs_seg*bs + pix].x,  h_liste_pix[(16*n + pos + 8)*nblocs_seg*bs + pix].y);
 628
 629                         }
 630
 631                 }
 632
 633                 }
 634
 635 #endif //DEBUG_LISTES
 636
 637   /*
 638
 639         Test du calcul des sommes partielles 'somblocs' faites par le kernel 'calcul_contribs_segments_blocs_full'
 640
 641    */
 642
 643 #ifdef DEBUG_SOMBLOCS
 644   printf("********* SOMMES PARTIELLES  ***********\n");
 645   printf("bs = %d - grid = %d -  intervalles = %d - nblocs_seg = %d - pairs = %d \n", bs, grid.x, n_interval, nblocs_seg, pairs);
 646   for (int n=0; n< n_interval; n++)
 647         {
 648           idx_n = 2*n + !pairs ;
 649           if (idx_n == 0) idx_nprec = nnodes - 1 ;
 650           else idx_nprec = idx_n - 1 ;
 651           if (idx_n == nnodes-1) idx_nsuiv = 0 ;
 652           else idx_nsuiv = idx_n + 1 ;
 653           printf("******** node %d\n", idx_n) ;
 654           for(int s=0; s<8; s++)
 655                 {
 656                   int nb_pix = calcul_liste_pixel_segment(h_snake[idx_nprec].posi,h_snake[idx_nprec].posj,
 657                                                                                                   h_liste_positions[8*idx_n+s].x, h_liste_positions[8*idx_n+s].y,
 658                                                                                                   h_liste_pixels_segment, 0);
 659                   for (int b=0; b<nblocs_seg; b++)
 660                         {
 661                           uint64 c1=0, cx=0, cx2=0 ;
 662                           int i,j;
 663                           for (int p=0; p<bs; p++)
 664                                 {
 665                                   if ( ((b*bs+p) < (nb_pix-1)) && ((b*bs+p)>0) )
 666                                         {
 667                                           //  /!\ penser a oter le test de prise en
 668                                           // compte pour les pix sur la même ligne dans
 669                                           // le kernel, sinon les comparaisons des
 670                                           // sommes par colonne seront fausses
 671                                           i = h_liste_pixels_segment[2*(b*bs + p)] ;
 672                                           j = h_liste_pixels_segment[2*(b*bs + p) + 1] ;
 673                                           c1 += img_1[i][j] ;
 674                                           cx += img_x[i][j] ;
 675                                           cx2+= img_x2[i][j];
 676                                         }
 677                                 }
 678                           if ( ( c1 != h_sombloc[(16*n + s)*nblocs_seg + b ] ) || ( cx != h_sombloc[(16*n + s)*nblocs_seg + b + grid.x] )
 679                                    ||  ( cx2 != h_sombloc[ (16*n + s)*nblocs_seg + b + 2*grid.x] ) )
 680                                 printf("seg %d - %d pix : bloc %d -> CPU : %lu - %lu - %lu \t|| GPU : %lu - %lu - %lu \n", s, nb_pix, b,
 681                                            c1, cx, cx2, h_sombloc[(16*n+s)*nblocs_seg + b], h_sombloc[(16*n+s)*nblocs_seg + b + grid.x],
 682                                            h_sombloc[(16*n+s)*nblocs_seg + b + 2*grid.x]) ;
 683                         }
 684
 685                 }
 686            for(int s=0; s<8; s++)
 687                 {
 688                   int nb_pix = calcul_liste_pixel_segment( h_liste_positions[8*idx_n+s].x, h_liste_positions[8*idx_n+s].y,
 689                                                                                                   h_snake[idx_nsuiv].posi,h_snake[idx_nsuiv].posj,
 690                                                                                                   h_liste_pixels_segment, 0);
 691                   for (int b=0; b<nblocs_seg; b++)
 692                         {
 693                           uint64 c1=0, cx=0, cx2=0 ;
 694                           int i,j;
 695                           for (int p=0; p<bs; p++)
 696                                 {
 697                                   if ( ((b*bs+p) < (nb_pix-1)) && ((b*bs+p)>0) )
 698                                         {
 699                                           //  /!\ penser a oter le test de prise en
 700                                           // compte pour les pix sur la même ligne dans
 701                                           // le kernel, sinon les comparaisons des
 702                                           // sommes par colonne seront fausses
 703                                           i = h_liste_pixels_segment[2*(b*bs + p)] ;
 704                                           j = h_liste_pixels_segment[2*(b*bs + p) + 1] ;
 705                                           c1 += img_1[i][j] ;
 706                                           cx += img_x[i][j] ;
 707                                           cx2+= img_x2[i][j];
 708                                         }
 709                                 }
 710                           if ( ( c1 != h_sombloc[(16*n + s + 8)*nblocs_seg + b ] ) || ( cx != h_sombloc[(16*n + s + 8)*nblocs_seg + b + grid.x] )
 711                                    ||  ( cx2 != h_sombloc[ (16*n + s + 8)*nblocs_seg + b + 2*grid.x] ) )
 712                                 printf("seg %d - %d pix : bloc %d -> CPU : %lu - %lu - %lu \t|| GPU : %lu - %lu - %lu \n", s, nb_pix, b,
 713                                            c1, cx, cx2, h_sombloc[(16*n+s+8)*nblocs_seg + b], h_sombloc[(16*n+s+8)*nblocs_seg + b + grid.x],
 714                                            h_sombloc[(16*n+s+8)*nblocs_seg + b + 2*grid.x]) ;
 715                         }
 716
 717                 }
 718
 719         }
 720 #endif //DEBUG_SOMBLOCS
 721
 722
 723  /*
 724
 725         Test du calcul des sommes totales 'somsom' faites par le kernel 'somsom_full'
 726
 727    */
 728
 729 #ifdef DEBUG_SOMSOM
 730  printf("********* SOMMES TOTALES  ***********\n");
 731   printf("bs = %d - grid = %d -  intervalles = %d - nblocs_seg = %d - pairs = %d \n", bs, grid.x, n_interval, nblocs_seg, pairs);
 732   for (int n=0; n< n_interval; n++)
 733         {
 734           idx_n = 2*n + !pairs ;
 735           if (idx_n == 0) idx_nprec = nnodes - 1 ;
 736           else idx_nprec = idx_n - 1 ;
 737           if (idx_n == nnodes-1) idx_nsuiv = 0 ;
 738           else idx_nsuiv = idx_n + 1 ;
 739           printf("******** node %d\n", idx_n) ;
 740           for(int s=0; s<8; s++)
 741                 {
 742                   int nb_pix = calcul_liste_pixel_segment(h_snake[idx_nprec].posi,h_snake[idx_nprec].posj,
 743                                                                                                   h_liste_positions[8*idx_n+s].x, h_liste_positions[8*idx_n+s].y,
 744                                                                                                   h_liste_pixels_segment, 0);
 745                   uint64 c1=0, cx=0, cx2=0 ;
 746                   for (int b=0; b<nblocs_seg; b++)
 747                         {
 748                           int i,j;
 749                           for (int p=0; p<bs; p++)
 750                                 {
 751                                   if ( ((b*bs+p) < (nb_pix-1)) && ((b*bs+p)>0) )
 752                                         {
 753                                           //  /!\ penser a oter le test de prise en
 754                                           // compte pour les pix sur la même ligne dans
 755                                           // le kernel, sinon les comparaisons des
 756                                           // sommes par colonne seront fausses
 757                                           i = h_liste_pixels_segment[2*(b*bs + p)] ;
 758                                           j = h_liste_pixels_segment[2*(b*bs + p) + 1] ;
 759                                           c1 += img_1[i][j] ;
 760                                           cx += img_x[i][j] ;
 761                                           cx2+= img_x2[i][j];
 762                                         }
 763                                 }
 764                         }
 765                   if ( ( c1 != h_somsom[3*(16*n + s)] ) || ( cx != h_somsom[3*(16*n + s) + 1] )
 766                            ||  ( cx2 != h_somsom[3*(16*n + s) + 2] ) )
 767                                 printf("seg %d - %d pix -> CPU : %lu - %lu - %lu \t|| GPU : %lu - %lu - %lu \n", s, nb_pix,
 768                                            c1, cx, cx2, h_somsom[3*(16*n + s)], h_somsom[3*(16*n + s) + 1],
 769                                            h_somsom[3*(16*n + s) + 2]) ;
 770
 771                 }
 772
 773            for(int s=0; s<8; s++)
 774                 {
 775                   int nb_pix = calcul_liste_pixel_segment( h_liste_positions[8*idx_n+s].x, h_liste_positions[8*idx_n+s].y,
 776                                                                                                   h_snake[idx_nsuiv].posi,h_snake[idx_nsuiv].posj,
 777                                                                                                    h_liste_pixels_segment, 0);
 778                   uint64 c1=0, cx=0, cx2=0 ;
 779                   for (int b=0; b<nblocs_seg; b++)
 780                         {
 781
 782                           int i,j;
 783                           for (int p=0; p<bs; p++)
 784                                 {
 785                                   if ( ((b*bs+p) < (nb_pix-1)) && ((b*bs+p)>0) )
 786                                         {
 787                                           //  /!\ penser a oter le test de prise en
 788                                           // compte pour les pix sur la même ligne dans
 789                                           // le kernel, sinon les comparaisons des
 790                                           // sommes par colonne seront fausses
 791                                           i = h_liste_pixels_segment[2*(b*bs + p)] ;
 792                                           j = h_liste_pixels_segment[2*(b*bs + p) + 1] ;
 793                                           c1 += img_1[i][j] ;
 794                                           cx += img_x[i][j] ;
 795                                           cx2+= img_x2[i][j];
 796                                         }
 797                                 }
 798                         }
 799                   if ( ( c1 != h_somsom[3*(16*n + s + 8)]  ) || ( cx != h_somsom[3*(16*n + s + 8) + 1] )
 800                            ||  ( cx2 != h_somsom[3*(16*n + s + 8) + 2] ) )
 801                         printf("seg %d - %d pix -> CPU : %lu - %lu - %lu \t|| GPU : %lu - %lu - %lu \n", s, nb_pix,
 802                                    c1, cx, cx2, h_somsom[3*(16*n + s + 8)], h_somsom[3*(16*n + s + 8) + 1],
 803                                    h_somsom[3*(16*n + s + 8)  + 2]) ;
 804
 805                 }
 806
 807         }
 808
 809 #endif
 810
 811
 812 #ifdef DEBUG_MV
 813   printf("**** STATS - REF : %lf \n", *h_vrais_snake);
 814   for(int n=0; n<n_interval; n++)
 815         {
 816           for(int p=0; p<8; p++)
 817                 {
 818                   printf("test %d du node %d : %lu - %lu - %lu - - log_vrais = %lf\n", p, (2*n + !pairs),
 819                                  h_stats[3*(8*n+p)], h_stats[3*(8*n+p)+1], h_stats[3*(8*n+p)+2], h_vrais[8*n+p]);
 820                 }
 821         }
 822 #endif //DEBUG_MV
 823
 824
 825 #ifdef DEBUG_CRST
 826   printf("**** CROISEMENTS \n");
 827   for(int n=0; n<nnodes; n++)
 828         {
 829           printf("test du seg %d : ",  n);
 830           if ( h_croist[n] ) printf("CROISEMENT\n"); else printf("\n");
 831         }
 832 #endif //DEBUG_CRST
 833
 834
 835 #ifdef DEBUG_MOVE
 836   printf("**** MOUVEMENTS \n");
 837   for(int n=0; n<nnodes; n++)
 838         {
 839           printf("Node %d : (%s) ",n, (h_move[n])? "yes":"no");
 840         }
 841 #endif //DEBUG_MOVE
 842
 843   delete h_liste_positions ;
 844   delete h_snake;
 845
 846   /*
 847    * fin generation snake GPU
 848    */
 849 }
 850
 851