src/lib_gpu.cu

   1
   2 #include <stdio.h>
   3
   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
  28
  29 void cuda_init_img_cumul(unsigned short ** img_in, int H, int L, int nb_nodes,
  30                                                  unsigned short ** d_img, t_cumul_x ** d_img_x, t_cumul_x2 ** d_img_x2,
  31                                                  int ** d_freemanDiDj, int ** d_codeNoeud,
  32                                                  snake_node_gpu ** d_snake, uint32 ** d_nb_pix_max,
  33                                                  uint4 ** d_positions, uint64 ** d_contribs_segments, uint4 ** d_freemans_centres,
  34                                                  int ** d_codes_segments, int64 ** d_stats_snake,
  35                                                  int64 ** d_stats, int64 ** d_stats_ref, double ** d_vrais, double ** d_vrais_snake,
  36                                                  uint2 ** d_liste_pixels, uint64 ** d_contribs_segments_blocs,
  37                                                  bool ** d_move
  38                                                  )
  39 {
  40   unsigned int taille = H*L;
  41   timeval chrono;
  42
  43
  44   //allocation cumuls en memoire GPU
  45   tic(&chrono, NULL);
  46   /*
  47         MAX_PIX 20000
  48         MAX_NODES 10000
  49         MAX_LISTE_PIX 10000000
  50    */
  51   cudaMalloc( (void**) d_snake, MAX_NODES*sizeof(snake_node_gpu) );
  52
  53   cudaMalloc( (void**) d_img, taille*sizeof(unsigned short) );
  54   cudaMalloc( (void**) d_img_x, taille*sizeof(t_cumul_x) );
  55   cudaMalloc( (void**) d_img_x2, taille*sizeof(t_cumul_x2) );
  56
  57   cudaMalloc( (void**) d_freemanDiDj, 9*sizeof(int) );
  58   cudaMalloc( (void**) d_codeNoeud, 64*sizeof(int) );
  59
  60   cudaMalloc( (void**) d_stats_snake, 6*sizeof(int64)) ;
  61   cudaMalloc( (void**) d_positions, 8*MAX_NODES*sizeof(uint4)) ;
  62   cudaMalloc( (void**) d_contribs_segments, 3*16*MAX_NODES*sizeof(uint64)) ;
  63   cudaMalloc( (void**) d_contribs_segments_blocs, (3*MAX_LISTE_PIX/32)*sizeof(uint64)) ;
  64   cudaMalloc( (void**) d_freemans_centres, 16*MAX_NODES*sizeof(uint4)) ;
  65   cudaMalloc( (void**) d_codes_segments, 16*MAX_NODES*sizeof(int)) ;
  66   cudaMalloc( (void**) d_stats, 3*8*MAX_NODES*sizeof(int64)) ;
  67   cudaMalloc( (void**) d_stats_ref, 3*MAX_NODES*sizeof(int64)) ;
  68   cudaMalloc( (void**) d_vrais, 8*MAX_NODES*sizeof(double)) ;
  69   cudaMalloc( (void**) d_move, MAX_NODES*sizeof(bool)) ;
  70   cudaMalloc( (void**) d_nb_pix_max, sizeof(uint32)) ;
  71   cudaMalloc( (void**) d_vrais_snake, sizeof(double)) ;
  72
  73   cudaMalloc( (void**) d_liste_pixels, 16*5*(MAX_NODES)*sizeof(uint2) );
  74
  75   printf("TOTAL MEM = %ld octets\n",
  76                  (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))
  77                  +(MAX_LISTE_PIX)*(sizeof(uint2)+1)
  78                  +taille*(8+sizeof(t_cumul_x)+sizeof(t_cumul_x2))
  79                   +9*4+64*4+6*8+4+sizeof(double)) );
  80
  81   int64 * h_stats_snake = new int64[6];
  82
  83   toc(chrono, "temps alloc mem GPU");
  84
  85   /*detection-choix-initialisation de la carte GPU*/
  86   tic(&chrono, NULL) ;
  87   cudaDeviceProp deviceProp;
  88   deviceProp.major = 2;
  89   deviceProp.minor = 0;
  90   int dev;
  91   cudaChooseDevice(&dev, &deviceProp);
  92   cudaGetDeviceProperties(&deviceProp, dev);
  93   if(deviceProp.major >= 2 )
  94         {
  95           printf("Using Device %d: \"%s\"\n", dev, deviceProp.name);
  96           cudaSetDevice(dev);
  97         }
  98   toc(chrono, "temps acces GPU") ;
  99
 100   //copie tables correspondances freeman en mem GPU
 101   tic(&chrono, NULL) ;
 102   cudaMemcpy( *d_freemanDiDj, CORRESPONDANCE_Di_Dj_FREEMAN , 9*sizeof(int), cudaMemcpyHostToDevice);
 103   cudaMemcpy( *d_codeNoeud, TABLE_CODAGE , 64*sizeof(unsigned int), cudaMemcpyHostToDevice);
 104   toc(chrono, "temps transfert tables de codage") ;
 105
 106   /*transfert image en global mem GPU*/
 107   tic(&chrono, NULL);
 108   cudaMemcpy( *d_img, img_in[0], taille*sizeof(unsigned short), cudaMemcpyHostToDevice);
 109   toc(chrono, "transfert image vers GPU");
 110
 111   //calculs images cumulees sur GPU
 112   int blocs_max = 65536 ;
 113   int bs = 256 ; //arbitraire, d'apres les observations c'est souvent l'optimu
 114   unsigned int base = 0 ;
 115   unsigned int bl_l  = (L+bs-1)/bs ;
 116   unsigned int nb_lines =  blocs_max / bl_l ;
 117   unsigned int lines ;
 118   unsigned int tranches = ( 1 + H / nb_lines ) ;
 119   nb_lines = (H +tranches -1)/ tranches ; // equilibre la taille des tranches
 120
 121   dim3 threads(bs,1,1);
 122   int smem = nextPow2(bl_l)*2; //smem pour le prefixscan des sommes de blocs (etape 2)
 123   smem += smem >> DEC;
 124   smem += smem >> DEC;
 125   int smem_size = smem*sizeof(uint64);
 126   uint64 * d_somblocs ; // sommes des cumuls par bloc de calcul
 127
 128
 129   if(DEBUG_IMG_CUMUL)
 130         {
 131           printf("--- CALCULS IMAGES CUMULEES+STATS GPU  ----\n");
 132           printf("\t%d threads par bloc  -- %u blocs par ligne -- %u tranches -- %u lignes par tranche \n",bs, bl_l, tranches,nb_lines);
 133           printf(" Smem totale pour cumuls : %d\n", CFI(bs)*(sizeof(t_cumul_x)+sizeof(t_cumul_x2)) );
 134           tic(&chrono, NULL);
 135         }
 136   //calculs cumuls generiques : necessitent 3 etapes / 3 kernels
 137   cudaMalloc( (void**) &d_somblocs, 2*bl_l*nb_lines*sizeof(uint64) );
 138   cudaFuncSetCacheConfig(calcul_cumuls_gpu, cudaFuncCachePreferShared);
 139   do
 140         {
 141           if  ( H-base < nb_lines ) lines = H - base ; else lines = nb_lines ;
 142           printf("base = ligne %d -- traitement de %d lignes \n", base, lines) ;
 143           dim3 grid(bl_l*lines,1,1) ;
 144           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) ;
 145           scan_somblocs<<<2*lines, nextPow2(bl_l)/2, smem_size>>>(d_somblocs, bl_l) ;
 146           add_soms_to_cumuls<<<grid, threads>>>(*d_img_x, *d_img_x2, H, L, d_somblocs, bl_l, base, lines) ;
 147           base += lines ;
 148         }
 149   while (base < H) ;
 150   cudaFree(d_somblocs) ;
 151
 152   //calcul des sommes totales N, sigX et sigX2 sur l'image
 153   calcul_stats_image<<<1, 1>>>( *d_img_x, *d_img_x2, H, L, (uint64*)*d_stats_snake);
 154
 155
 156   cudaThreadSynchronize()   ;
 157   toc(chrono, "\tTemps GPU");
 158   if(DEBUG_IMG_CUMUL)
 159         {
 160
 161           //allocation memoire CPU
 162           t_cumul_x  * img_x = new t_cumul_x [H*L];
 163           t_cumul_x2 *  img_x2 = new t_cumul_x2 [H*L];
 164
 165           /*pour test comparaison*/
 166           t_cumul_x * img_xb = new t_cumul_x [H*L];
 167           t_cumul_x2 * img_x2b = new t_cumul_x2 [H*L];
 168
 169           cudaMemcpy( img_xb, *d_img_x, taille*sizeof(t_cumul_x), cudaMemcpyDeviceToHost);
 170           cudaMemcpy( img_x2b, *d_img_x2, taille*sizeof(t_cumul_x2), cudaMemcpyDeviceToHost);
 171
 172           //cumuls : etape 1 CPU
 173           /*
 174                 for (int i=0; i<H; i++)
 175                 {
 176                         for (int b=0; b<bl_l; b++)
 177                         {
 178                                 int offset = b*bs ;
 179                                 img_x[i*L+offset] = img_in[i][offset] ;
 180                                 img_x2[i*L+offset]= img_in[i][offset]*img_in[i][offset] ;
 181                                 for (int p=1; p<bs; p++)
 182                                 {
 183                                         int j = p+offset ;
 184                                         if (j<L)
 185                                         {
 186                                                 img_x[i*L+j] = img_x[i*L+j-1] + img_in[i][j];
 187                                                 img_x2[i*L+j] = img_x2[i*L+j-1] + img_in[i][j]*img_in[i][j] ;
 188                                         }
 189                                 }
 190                         }
 191                 }
 192           */
 193           //cumuls complets CPU
 194
 195           for (int i=0; i<H; i++)
 196                 {
 197                   img_x[i*L+0] = img_in[i][0] ;
 198                   img_x2[i*L+0]= img_in[i][0]*img_in[i][0] ;
 199                   for (int j=1; j<L; j++)
 200                         {
 201                           img_x[i*L+j]  = img_x[i*L+j-1]  + img_in[i][j] ;
 202                           img_x2[i*L+j] = img_x2[i*L+j-1] + img_in[i][j]*img_in[i][j] ;
 203                         }
 204                 }
 205
 206           int cpt = 0;
 207           int cpt_errx=0, cpt_errx2 = 0;
 208           for (int i=0; i< H; i++){
 209                 for (int j=0; j< L; j++){
 210                   if ( (img_x[i*L+j] !=  img_xb[i*L+j]) ) cpt_errx++ ;
 211                   if ( (img_x2[i*L+j] !=  img_x2b[i*L+j]) ) cpt_errx2++ ;
 212                   if ( (img_x[i*L+j] !=  img_xb[i*L+j]) || (img_x2[i*L+j] !=  img_x2b[i*L+j]))
 213                   {
 214                         //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]);
 215                   }
 216                   cpt++;
 217                 }
 218           }
 219           printf("%d erreurs sur CX / %d points\n", cpt_errx, cpt );
 220           printf("%d erreurs sur CX2 / %d points\n", cpt_errx2, cpt );
 221           uint64 sigX = 0, sigX2 = 0 ;
 222           for (int i=0; i<H; i++)
 223                 {
 224                   sigX += img_x[i*L+L-1] ;
 225                   sigX2+= img_x2[i*L+L-1];
 226                 }
 227           printf("STATS IMAGE  N = %d - sigX = %lu - sigX2 = %lu\n",  H*L, sigX, sigX2 );
 228         }
 229
 230   /*
 231    * generation snake en mem GPU
 232    */
 233   int dist = 140 ;
 234
 235   /* Test de determination du snake rectangle initial optimal*/
 236   int div = 100;//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
 237   int Nperm = div*div*bs;//nb total de rectangles a tester. La distribution est ainsi irrégulière, mais plus simple.
 238   double best_crit ;
 239   int ind_best_crit ;
 240
 241   t_rectangle_snake * d_all_crit, d_best_crit;//tableaux pour les résultats des différents rectangles / le meilleur
 242   t_rectangle_snake * h_all_crit = new t_rectangle_snake[Nperm];//correspondant CPU
 243
 244   //allocations
 245   cudaMalloc((void**) &d_all_crit, Nperm*sizeof(t_rectangle_snake));
 246   cudaMalloc((void**) &d_best_crit, sizeof(t_rectangle_snake));
 247
 248   tic(&chrono, NULL);
 249
 250   //execution kernel
 251   dim3 grid = dim3(H/div, L/div, 1);
 252   calcul_contribs_snake4<<<grid, bs, CFI(bs)*sizeof(tcontribs) >>>(*d_snake, *d_img_x, *d_img_x2, H, L, *d_stats_snake, d_all_crit) ;
 253   cudaThreadSynchronize();
 254   toc(chrono, "\nCALCULS RECTANGLES");
 255
 256   //recup data rectangles
 257   int ret;
 258   ret = cudaMemcpy( h_all_crit, d_all_crit, Nperm*sizeof(t_rectangle_snake), cudaMemcpyDeviceToHost) ;
 259   printf("COPIE DATA = %s\n",(ret==0)?"OK":"ERR");
 260
 261   //optimum sur CPU
 262   best_crit = h_all_crit[0].crit ;
 263   ind_best_crit = 0 ;
 264   for (int k=1; k<100; k++){
 265         if ((h_all_crit[k].crit > 0) && (h_all_crit[k].crit < best_crit)) {
 266           best_crit = h_all_crit[k].crit ;
 267           ind_best_crit = k ;
 268         }
 269         printf("%d -> ( %d, %d )--( %d, %d)  CRITERE = %f\n", k, h_all_crit[k].bpi, h_all_crit[k].bpj,
 270                  h_all_crit[k].opi, h_all_crit[k].opj,  h_all_crit[k].crit );
 271   }
 272
 273   printf("BEST RECTANGLE/%d tests : %d -> ( %d, %d )--( %d, %d)  CRITERE = %f\n", Nperm, ind_best_crit, h_all_crit[ind_best_crit].bpi, h_all_crit[ind_best_crit].bpj,
 274                  h_all_crit[ind_best_crit].opi, h_all_crit[ind_best_crit].opj, best_crit );
 275
 276   exit(0);
 277   /*fin test snake rectangle initial optimal*/
 278
 279   //genere_snake_rectangle_4nodes_gpu<<< 1, 1>>>(*d_snake, 140, H, L) ;
 280
 281
 282   int nnodes = nb_nodes ;
 283   snake_node_gpu * h_snake = new snake_node_gpu[nnodes];
 284   snake_node * h_snake_ll = new snake_node[nnodes] ;
 285   uint4 * h_liste_positions = new uint4[nnodes*8];
 286   double * h_vrais_snake = new double ;
 287   //init les stats du snake
 288   uint2 * d_liste_temp  ;
 289   t_sum_x2 * d_sompart  ;
 290   int tpb, bps, npixmax ;
 291
 292   //calcul nb threads par bloc
 293   npixmax = 2*(H+L-4*dist)/(nnodes-1) ;
 294   tpb = nextPow2(npixmax) ;
 295   if (tpb >= 256) tpb = 256 ;//  /!\ le kernel <<< calcul_contrib...>>> ne supporte pas un bs>256 a cause de la shared-mem nécessaire
 296   if (tpb < 32 ) tpb = 32 ;
 297   tpb=128 ;
 298   bps = (npixmax+tpb-1)/tpb ;
 299   printf("PARAMS EXEC INIT : %d pix max, %d threads/bloc, %d blocs/seg, %d blocs/grille\n", npixmax, tpb, bps, nnodes*bps);
 300
 301   //alloc
 302   cudaMalloc((void**) &d_liste_temp, nnodes*bps*tpb*sizeof(uint2));
 303   cudaMalloc((void**) &d_sompart, 3*nnodes*bps*sizeof(t_sum_x2));
 304   cudaMalloc((void**) &d_stats_ref, 3*nnodes*sizeof(int64));
 305
 306   //DEBUG : pour forcer la mise à zero du tableau intermediaire d_stats_ref
 307   /*
 308   int64 h_stats_ref[3*nnodes] ;
 309   for (int a=0; a<3*nnodes ; a++) h_stats_ref[a] = 0 ;
 310   cudaMemcpy( h_stats_ref, d_stats_ref, sizeof(int64), cudaMemcpyHostToDevice) ;
 311   */
 312   //fin forçage a 0
 313
 314   //DEBUG : pour forcer la mise à zero du tableau intermediaire d_sompart
 315   /*
 316      t_sum_x2 h_sompart[ 3*nnodes*bps ] ;
 317      for (int a=0; a<3*nnodes*bps ; a++) h_sompart[a] = 0 ;
 318      cudaMemcpy( h_sompart, d_sompart, sizeof(t_sum_x2), cudaMemcpyHostToDevice) ;
 319   */
 320   //fin forçage a 0
 321
 322   calcul_contribs_segments_snake<<< nnodes*bps, tpb, (CFI(tpb))*(3*sizeof(t_sum_x2))>>>
 323         (*d_snake, nnodes,
 324          *d_img_x, *d_img_x2,
 325          L, d_liste_temp, d_sompart, *d_freemanDiDj );
 326
 327   //TODO
 328   //parametrer pour ne pas appeler qd tpb=1
 329   //oblige a modifier le kernel <<< calcul_contrib...>>> pour ecrire directement ds d_snake
 330   // au lieu de d_sompart
 331   somsom_snake<<< nnodes , 1 >>>(d_sompart, nnodes, bps, *d_snake);
 332
 333
 334   calcul_stats_snake<<< 1 , 1 >>>(*d_snake, nnodes, *d_stats_snake, *d_vrais_snake,
 335                                                                   *d_img_x, *d_img_x2,
 336                                                                   *d_codeNoeud, L
 337                                                                   );
 338   cudaThreadSynchronize() ;
 339   toc(chrono, "\tTemps") ;
 340
 341   /*
 342         verif stats initiales du snake
 343   */
 344   cudaMemcpy( h_vrais_snake, *d_vrais_snake, sizeof(double), cudaMemcpyDeviceToHost) ;
 345   cudaMemcpy( h_stats_snake, *d_stats_snake, 6*sizeof(int64), cudaMemcpyDeviceToHost) ;
 346
 347   printf("STATS SNAKE log vrais=%lf : c1=%lu - cx=%lu - cx2=%lu - N=%lu - SUMX=%lu - SUMX2=%lu\n",
 348                  *h_vrais_snake,
 349                  h_stats_snake[0],  h_stats_snake[1],  h_stats_snake[2],
 350                  h_stats_snake[3],  h_stats_snake[4],  h_stats_snake[5] );
 351
 352   /*
 353         verif stats diminuees des contribs des 2 segments associes a chaque noeud
 354   */
 355 #ifdef DEBUG_STATS_REF
 356   cudaMemcpy( h_stats_ref, d_stats_ref, 3*nnodes*sizeof(int64), cudaMemcpyDeviceToHost) ;
 357   cudaMemcpy( h_snake, *d_snake, nnodes*sizeof(snake_node_gpu), cudaMemcpyDeviceToHost) ;
 358
 359
 360   printf("******* STATS DIMINUEES\n");
 361   for(int n=0; n<nnodes;n++)
 362         {
 363           int i = h_snake[n].posi, j = h_snake[n].posj ;
 364           printf("node %d (%d,%d) : %ld - %ld - %ld - img1= %lu - imgx= %lu - imgx2= %lu \n", n, i, j,
 365                          h_stats_ref[3*n], h_stats_ref[3*n +1], h_stats_ref[3*n +2],
 366                          img_1[i][j], img_x[i][j], img_x2[i][j]);
 367         }
 368 #endif //DEBUG_STATS_REF
 369
 370   //snake2gpu(d_snake, snake, nb_nodes);
 371   //gpu2snake(*d_snake, &h_snake_ll, nnodes);
 372
 373
 374 #ifdef DEBUG_POSITIONS
 375   for (int n=0; n<nnodes; n++)
 376         {
 377           printf("Node %d :\n", n);
 378           for (int pos=0; pos<8; pos++)
 379                 {
 380                   printf("(%d , %d):%d:%d | ", h_liste_positions[8*n + pos].x, h_liste_positions[8*n + pos].y,
 381                                  h_liste_positions[8*n + pos].z, h_liste_positions[8*n + pos].w);
 382                 }
 383           printf("\n");
 384         }
 385 #endif //DEBUG_POSITIONS
 386
 387 //verif liste pixels noeuds pairs/impairs selon
 388
 389 #ifdef DEBUG_LISTES
 390   printf("NOMBRE PIXELS pour LISTE = %d\n", *h_nb_pix_max) ;
 391   printf("bs = %d - grid = %d - nblocs_seg = %d - npix_max = %d - taille mem = %d\n",
 392                  bs, grid.x, nblocs_seg, *h_nb_pix_max, taille_mem);
 393
 394   cudaMemcpy( h_liste_pix, d_liste_pix, taille_mem*sizeof(uint2), cudaMemcpyDeviceToHost ) ;
 395   cudaMemcpy( h_snake, *d_snake, nnodes*sizeof(snake_node_gpu), cudaMemcpyDeviceToHost );
 396   uint32 * h_liste_pixels_segment = new uint32[2*(*h_nb_pix_max)] ;
 397   int idx_n, idx_nprec, idx_nsuiv ;
 398
 399   printf("********* LISTE PIX  ***********\n");
 400   printf("bs = %d - grid = %d - nblocs_seg = %d - npix_max = %d - taille mem = %d\n",
 401                  bs, grid.x, nblocs_seg, *h_nb_pix_max, taille_mem);
 402
 403   for (int n=0; n<(nnodes/2 + (nnodes%2)*pairs); n++)
 404         {
 405           idx_n = 2*n + !pairs ;
 406           if (idx_n == 0) idx_nprec = nnodes - 1;
 407           else idx_nprec = idx_n - 1;
 408           if (idx_n == nnodes-1) idx_nsuiv = 0;
 409           else idx_nsuiv = idx_n + 1 ;
 410
 411           for (int pos=0; pos < 8 ; pos++) //test des segments avant le noeud
 412                 {
 413
 414                   int nb_pix = calcul_liste_pixel_segment(h_snake[idx_nprec].posi,h_snake[idx_nprec].posj,
 415                                                                                                   h_liste_positions[8*idx_n+pos].x, h_liste_positions[8*idx_n+pos].y,
 416                                                                                                   h_liste_pixels_segment, 0);
 417
 418                   for (int pix=0; pix < nb_pix; pix++)
 419                         {
 420
 421                           if ( (h_liste_pix[(16*n + pos)*nblocs_seg*bs + pix].x != h_liste_pixels_segment[2*pix] )
 422                                    || ( h_liste_pix[(16*n + pos)*nblocs_seg*bs + pix].y != h_liste_pixels_segment[2*pix+1] ) )
 423                                 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,
 424                                            idx_nprec, idx_n, idx_nsuiv,
 425                                            h_liste_pixels_segment[2*pix], h_liste_pixels_segment[2*pix+1],
 426                                            h_liste_pix[(16*n + pos)*nblocs_seg*bs + pix].x,  h_liste_pix[(16*n + pos)*nblocs_seg*bs + pix].y);
 427
 428                         }
 429
 430                 }
 431           for (int pos=0; pos < 8 ; pos++) //test des segments apres le noeud
 432                 {
 433
 434                   int nb_pix = calcul_liste_pixel_segment(h_liste_positions[8*idx_n+pos].x, h_liste_positions[8*idx_n+pos].y,
 435                                                                                                   h_snake[idx_nsuiv].posi,h_snake[idx_nsuiv].posj,
 436                                                                                                   h_liste_pixels_segment, 0);
 437
 438                   for (int pix=0; pix < nb_pix; pix++)
 439                         {
 440
 441                           if ( (h_liste_pix[(16*n + pos + 8)*nblocs_seg*bs + pix].x != h_liste_pixels_segment[2*pix] )
 442                              || ( h_liste_pix[(16*n + pos + 8)*nblocs_seg*bs + pix].y != h_liste_pixels_segment[2*pix+1] ) )
 443                                 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,
 444                                            idx_nprec, idx_n, idx_nsuiv,
 445                                            h_liste_pixels_segment[2*pix], h_liste_pixels_segment[2*pix+1],
 446                                            h_liste_pix[(16*n + pos + 8)*nblocs_seg*bs + pix].x,  h_liste_pix[(16*n + pos + 8)*nblocs_seg*bs + pix].y);
 447
 448                         }
 449
 450                 }
 451
 452                 }
 453
 454 #endif //DEBUG_LISTES
 455
 456   /*
 457
 458         Test du calcul des sommes partielles 'somblocs' faites par le kernel 'calcul_contribs_segments_blocs_full'
 459
 460    */
 461
 462 #ifdef DEBUG_SOMBLOCS
 463   printf("********* SOMMES PARTIELLES  ***********\n");
 464   printf("bs = %d - grid = %d -  intervalles = %d - nblocs_seg = %d - pairs = %d \n", bs, grid.x, n_interval, nblocs_seg, pairs);
 465   for (int n=0; n< n_interval; n++)
 466         {
 467           idx_n = 2*n + !pairs ;
 468           if (idx_n == 0) idx_nprec = nnodes - 1 ;
 469           else idx_nprec = idx_n - 1 ;
 470           if (idx_n == nnodes-1) idx_nsuiv = 0 ;
 471           else idx_nsuiv = idx_n + 1 ;
 472           printf("******** node %d\n", idx_n) ;
 473           for(int s=0; s<8; s++)
 474                 {
 475                   int nb_pix = calcul_liste_pixel_segment(h_snake[idx_nprec].posi,h_snake[idx_nprec].posj,
 476                                                                                                   h_liste_positions[8*idx_n+s].x, h_liste_positions[8*idx_n+s].y,
 477                                                                                                   h_liste_pixels_segment, 0);
 478                   for (int b=0; b<nblocs_seg; b++)
 479                         {
 480                           uint64 c1=0, cx=0, cx2=0 ;
 481                           int i,j;
 482                           for (int p=0; p<bs; p++)
 483                                 {
 484                                   if ( ((b*bs+p) < (nb_pix-1)) && ((b*bs+p)>0) )
 485                                         {
 486                                           //  /!\ penser a oter le test de prise en
 487                                           // compte pour les pix sur la même ligne dans
 488                                           // le kernel, sinon les comparaisons des
 489                                           // sommes par colonne seront fausses
 490                                           i = h_liste_pixels_segment[2*(b*bs + p)] ;
 491                                           j = h_liste_pixels_segment[2*(b*bs + p) + 1] ;
 492                                           c1 += img_1[i][j] ;
 493                                           cx += img_x[i][j] ;
 494                                           cx2+= img_x2[i][j];
 495                                         }
 496                                 }
 497                           if ( ( c1 != h_sombloc[(16*n + s)*nblocs_seg + b ] ) || ( cx != h_sombloc[(16*n + s)*nblocs_seg + b + grid.x] )
 498                                    ||  ( cx2 != h_sombloc[ (16*n + s)*nblocs_seg + b + 2*grid.x] ) )
 499                                 printf("seg %d - %d pix : bloc %d -> CPU : %lu - %lu - %lu \t|| GPU : %lu - %lu - %lu \n", s, nb_pix, b,
 500                                            c1, cx, cx2, h_sombloc[(16*n+s)*nblocs_seg + b], h_sombloc[(16*n+s)*nblocs_seg + b + grid.x],
 501                                            h_sombloc[(16*n+s)*nblocs_seg + b + 2*grid.x]) ;
 502                         }
 503
 504                 }
 505            for(int s=0; s<8; s++)
 506                 {
 507                   int nb_pix = calcul_liste_pixel_segment( h_liste_positions[8*idx_n+s].x, h_liste_positions[8*idx_n+s].y,
 508                                                                                                   h_snake[idx_nsuiv].posi,h_snake[idx_nsuiv].posj,
 509                                                                                                   h_liste_pixels_segment, 0);
 510                   for (int b=0; b<nblocs_seg; b++)
 511                         {
 512                           uint64 c1=0, cx=0, cx2=0 ;
 513                           int i,j;
 514                           for (int p=0; p<bs; p++)
 515                                 {
 516                                   if ( ((b*bs+p) < (nb_pix-1)) && ((b*bs+p)>0) )
 517                                         {
 518                                           //  /!\ penser a oter le test de prise en
 519                                           // compte pour les pix sur la même ligne dans
 520                                           // le kernel, sinon les comparaisons des
 521                                           // sommes par colonne seront fausses
 522                                           i = h_liste_pixels_segment[2*(b*bs + p)] ;
 523                                           j = h_liste_pixels_segment[2*(b*bs + p) + 1] ;
 524                                           c1 += img_1[i][j] ;
 525                                           cx += img_x[i][j] ;
 526                                           cx2+= img_x2[i][j];
 527                                         }
 528                                 }
 529                           if ( ( c1 != h_sombloc[(16*n + s + 8)*nblocs_seg + b ] ) || ( cx != h_sombloc[(16*n + s + 8)*nblocs_seg + b + grid.x] )
 530                                    ||  ( cx2 != h_sombloc[ (16*n + s + 8)*nblocs_seg + b + 2*grid.x] ) )
 531                                 printf("seg %d - %d pix : bloc %d -> CPU : %lu - %lu - %lu \t|| GPU : %lu - %lu - %lu \n", s, nb_pix, b,
 532                                            c1, cx, cx2, h_sombloc[(16*n+s+8)*nblocs_seg + b], h_sombloc[(16*n+s+8)*nblocs_seg + b + grid.x],
 533                                            h_sombloc[(16*n+s+8)*nblocs_seg + b + 2*grid.x]) ;
 534                         }
 535
 536                 }
 537
 538         }
 539 #endif //DEBUG_SOMBLOCS
 540
 541
 542  /*
 543
 544         Test du calcul des sommes totales 'somsom' faites par le kernel 'somsom_full'
 545
 546    */
 547
 548 #ifdef DEBUG_SOMSOM
 549  printf("********* SOMMES TOTALES  ***********\n");
 550   printf("bs = %d - grid = %d -  intervalles = %d - nblocs_seg = %d - pairs = %d \n", bs, grid.x, n_interval, nblocs_seg, pairs);
 551   for (int n=0; n< n_interval; n++)
 552         {
 553           idx_n = 2*n + !pairs ;
 554           if (idx_n == 0) idx_nprec = nnodes - 1 ;
 555           else idx_nprec = idx_n - 1 ;
 556           if (idx_n == nnodes-1) idx_nsuiv = 0 ;
 557           else idx_nsuiv = idx_n + 1 ;
 558           printf("******** node %d\n", idx_n) ;
 559           for(int s=0; s<8; s++)
 560                 {
 561                   int nb_pix = calcul_liste_pixel_segment(h_snake[idx_nprec].posi,h_snake[idx_nprec].posj,
 562                                                                                                   h_liste_positions[8*idx_n+s].x, h_liste_positions[8*idx_n+s].y,
 563                                                                                                   h_liste_pixels_segment, 0);
 564                   uint64 c1=0, cx=0, cx2=0 ;
 565                   for (int b=0; b<nblocs_seg; b++)
 566                         {
 567                           int i,j;
 568                           for (int p=0; p<bs; p++)
 569                                 {
 570                                   if ( ((b*bs+p) < (nb_pix-1)) && ((b*bs+p)>0) )
 571                                         {
 572                                           //  /!\ penser a oter le test de prise en
 573                                           // compte pour les pix sur la même ligne dans
 574                                           // le kernel, sinon les comparaisons des
 575                                           // sommes par colonne seront fausses
 576                                           i = h_liste_pixels_segment[2*(b*bs + p)] ;
 577                                           j = h_liste_pixels_segment[2*(b*bs + p) + 1] ;
 578                                           c1 += img_1[i][j] ;
 579                                           cx += img_x[i][j] ;
 580                                           cx2+= img_x2[i][j];
 581                                         }
 582                                 }
 583                         }
 584                   if ( ( c1 != h_somsom[3*(16*n + s)] ) || ( cx != h_somsom[3*(16*n + s) + 1] )
 585                            ||  ( cx2 != h_somsom[3*(16*n + s) + 2] ) )
 586                                 printf("seg %d - %d pix -> CPU : %lu - %lu - %lu \t|| GPU : %lu - %lu - %lu \n", s, nb_pix,
 587                                            c1, cx, cx2, h_somsom[3*(16*n + s)], h_somsom[3*(16*n + s) + 1],
 588                                            h_somsom[3*(16*n + s) + 2]) ;
 589
 590                 }
 591
 592            for(int s=0; s<8; s++)
 593                 {
 594                   int nb_pix = calcul_liste_pixel_segment( h_liste_positions[8*idx_n+s].x, h_liste_positions[8*idx_n+s].y,
 595                                                                                                   h_snake[idx_nsuiv].posi,h_snake[idx_nsuiv].posj,
 596                                                                                                    h_liste_pixels_segment, 0);
 597                   uint64 c1=0, cx=0, cx2=0 ;
 598                   for (int b=0; b<nblocs_seg; b++)
 599                         {
 600
 601                           int i,j;
 602                           for (int p=0; p<bs; p++)
 603                                 {
 604                                   if ( ((b*bs+p) < (nb_pix-1)) && ((b*bs+p)>0) )
 605                                         {
 606                                           //  /!\ penser a oter le test de prise en
 607                                           // compte pour les pix sur la même ligne dans
 608                                           // le kernel, sinon les comparaisons des
 609                                           // sommes par colonne seront fausses
 610                                           i = h_liste_pixels_segment[2*(b*bs + p)] ;
 611                                           j = h_liste_pixels_segment[2*(b*bs + p) + 1] ;
 612                                           c1 += img_1[i][j] ;
 613                                           cx += img_x[i][j] ;
 614                                           cx2+= img_x2[i][j];
 615                                         }
 616                                 }
 617                         }
 618                   if ( ( c1 != h_somsom[3*(16*n + s + 8)]  ) || ( cx != h_somsom[3*(16*n + s + 8) + 1] )
 619                            ||  ( cx2 != h_somsom[3*(16*n + s + 8) + 2] ) )
 620                         printf("seg %d - %d pix -> CPU : %lu - %lu - %lu \t|| GPU : %lu - %lu - %lu \n", s, nb_pix,
 621                                    c1, cx, cx2, h_somsom[3*(16*n + s + 8)], h_somsom[3*(16*n + s + 8) + 1],
 622                                    h_somsom[3*(16*n + s + 8)  + 2]) ;
 623
 624                 }
 625
 626         }
 627
 628 #endif
 629
 630
 631 #ifdef DEBUG_MV
 632   printf("**** STATS - REF : %lf \n", *h_vrais_snake);
 633   for(int n=0; n<n_interval; n++)
 634         {
 635           for(int p=0; p<8; p++)
 636                 {
 637                   printf("test %d du node %d : %lu - %lu - %lu - - log_vrais = %lf\n", p, (2*n + !pairs),
 638                                  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]);
 639                 }
 640         }
 641 #endif //DEBUG_MV
 642
 643
 644 #ifdef DEBUG_CRST
 645   printf("**** CROISEMENTS \n");
 646   for(int n=0; n<nnodes; n++)
 647         {
 648           printf("test du seg %d : ",  n);
 649           if ( h_croist[n] ) printf("CROISEMENT\n"); else printf("\n");
 650         }
 651 #endif //DEBUG_CRST
 652
 653
 654 #ifdef DEBUG_MOVE
 655   printf("**** MOUVEMENTS \n");
 656   for(int n=0; n<nnodes; n++)
 657         {
 658           printf("Node %d : (%s) ",n, (h_move[n])? "yes":"no");
 659         }
 660 #endif //DEBUG_MOVE
 661
 662   delete h_liste_positions ;
 663   delete h_snake;
 664
 665   /*
 666    * fin generation snake GPU
 667    */
 668 }
 669
 670