+ //execution kernels
+ calcul_contribs_cols<<<grid, bs, smem>>>( *d_img_x, *d_img_x2, H, L, d_contribs_part );
+ somsom_contribs<<<div,1>>>( d_contribs_part, bpc, d_contribs_cols ) ;
+ calcul_contribs_permutations<<<div,div, CFI(div)*sizeof(double2)>>>(d_contribs_cols, bpc, H, L, (uint64*) *d_stats_snake, d_miniblocs) ;
+ cudaMemcpy( h_miniblocs, d_miniblocs, div*sizeof(double2), cudaMemcpyDeviceToHost ) ;
+
+ //verif minimum des blocs
+ double crit_mini = h_miniblocs[0].x;
+ int j1=0, id_mini=0;
+ for (j1=1 ; j1 < div ; j1++){
+ if (h_miniblocs[j1].x < crit_mini) {
+ crit_mini = h_miniblocs[j1].x ;
+ id_mini = j1 ;
+ }
+ }
+ toc(chrono, "\nCALCUL RECTANGLE");
+ j1 = pas * id_mini ;
+ int j2 = (int)(pas*h_miniblocs[ id_mini ].y) ;
+ printf("pas = %d cols -- critere mini =%f , positions j1=%d j2=%d\n", pas, h_miniblocs[ id_mini ].x, j1, j2);
+
+
+ // transfert datas GPU -> CPU
+ cudaMemcpy( h_contribs_cols, d_contribs_cols, div*sizeof(tcontribs), cudaMemcpyDeviceToHost ) ;
+
+ //verif contribs colonnes
+ for (int c=0 ; c < div ; c++){
+ // calcul valeurs de ref en CPU
+ h_contribs_cols_cpu[c].cx = 0 ;
+ h_contribs_cols_cpu[c].cx2= 0 ;
+ for (int ip=0; ip < H; ip++){
+ h_contribs_cols_cpu[c].cx += img_x[ ip*L + c*pas] ;
+ h_contribs_cols_cpu[c].cx2 += img_x2[ ip*L + c*pas] ;
+ }
+ //comparaison avec valeurs GPU
+ if ( (h_contribs_cols_cpu[c].cx != h_contribs_cols[c].cx) || (h_contribs_cols_cpu[c].cx2 != h_contribs_cols[c].cx2) )
+ printf("ERR colonne %d -> CPUx=%lu CPUx2=%lu | GPUx=%lu GPUx2=%lu\n",
+ c*pas, h_contribs_cols_cpu[c].cx, h_contribs_cols_cpu[c].cx2, h_contribs_cols[c].cx, h_contribs_cols[c].cx2 );
+ }
+ cudaFree(d_contribs_part) ;
+ cudaFree(d_contribs_cols) ;
+ cudaFree(d_miniblocs) ;
+ free(h_contribs_cols);
+ free(h_contribs_cols_cpu);
+ free(h_miniblocs) ;
+
+ //realloc pour lignes horizontales
+ bs = 128 ;
+
+ div = (H+bs-1)/bs ;
+ printf("DIV = %d\n", div ) ;
+
+ int divpow2 = nextPow2(div) ;
+ printf("DIVPOW2 = %d\n", divpow2) ;
+
+ grid = dim3(div, 1, 1) ;
+ smem = CFI(bs)*sizeof(tcontribs) ;
+ cudaMalloc((void**) &d_contribs_part, div*sizeof(tcontribs)) ;
+ cudaMalloc((void**) &d_contribs_cols, div*div*sizeof(tcontribs)) ;
+ cudaMalloc((void**) &d_miniblocs, div*sizeof(double2)) ;
+ h_contribs_cols = new tcontribs[div*div] ;
+ tcontribs * h_contribs_part = new tcontribs[div] ;
+ h_miniblocs = new double2[div] ;
+
+ tic(&chrono, NULL);
+ // Appels kernels optim lignes horizontales
+ calcul_contrib_conjuguee_colonnes<<<grid, bs, CFI(bs)*sizeof(tcontribs) >>>( *d_img_x, *d_img_x2, H, L, j1, j2, d_contribs_part) ;
+
+ /*verif CPU
+ int cpt = 0 ;
+ int cpterr = 0 ;
+ tcontribs * h_contribs_part_cpu = new tcontribs[div] ;
+ cudaMemcpy( h_contribs_part, d_contribs_part, div*sizeof(tcontribs), cudaMemcpyDeviceToHost ) ;
+ for (int bloc=0; bloc < div; bloc++){
+ h_contribs_part_cpu[ bloc ].cx = 0 ;
+ h_contribs_part_cpu[ bloc ].cx2 = 0 ;
+ for (int line=0; ((line < bs)&&(bloc*bs+line < H)); line++){
+ h_contribs_part_cpu[bloc].cx += img_x[ (bloc*bs+line)*L + j2] - img_x[ (bloc*bs+line)*L + j1 ];
+ h_contribs_part_cpu[bloc].cx2 += img_x2[ (bloc*bs+line)*L + j2] - img_x2[ (bloc*bs+line)*L + j1 ];
+ }
+ if ( ( h_contribs_part_cpu[bloc].cx != h_contribs_part[bloc].cx ) || ( h_contribs_part_cpu[bloc].cx2 != h_contribs_part[bloc].cx2 ) )
+ {
+ printf("ERREUR bloc %d -> CPUx=%lu CPUx2=%lu | GPUx=%lu GPUx2=%lu\n", bloc,
+ h_contribs_part_cpu[bloc].cx, h_contribs_part_cpu[bloc].cx2, h_contribs_part[bloc].cx, h_contribs_part[bloc].cx2 ) ;
+ cpterr++;
+ }
+ cpt++ ;
+ }
+ printf("VERIF CONTRIB CONJUGUEES BLOCS --> %d ERREURS / %d BLOCS\n", cpterr, cpt) ;
+ fin verif*/
+
+ grid = dim3(div, div, 1) ;
+ calcul_contribs_snake_rectangle<<<grid,divpow2, CFI(divpow2)*sizeof(tcontribs) >>>(d_contribs_part, d_contribs_cols) ;
+
+ /* verif CPU
+ h_contribs_cols_cpu = new tcontribs[div*div] ;
+ cudaMemcpy( h_contribs_cols, d_contribs_cols, div*div*sizeof(tcontribs), cudaMemcpyDeviceToHost ) ;
+ cpt = 0 ;
+ cpterr = 0 ;
+ for (int i1=0; i1 < div ; i1++){
+ for (int i2=0 ; i2 < div ; i2++){
+ if (i2 >= i1){
+ h_contribs_cols_cpu[ i1*div +i2 ].cx = 0 ;
+ h_contribs_cols_cpu[ i1*div +i2 ].cx2= 0 ;
+ for (int d=i1 ; d <= i2 ; d++){
+ h_contribs_cols_cpu[ i1*div +i2 ].cx += h_contribs_part_cpu[ d ].cx ;
+ h_contribs_cols_cpu[ i1*div +i2 ].cx2+= h_contribs_part_cpu[ d ].cx2 ;
+ }
+ } else {
+ h_contribs_cols_cpu[ i1*div +i2 ].cx = 0 ;
+ h_contribs_cols_cpu[ i1*div +i2 ].cx2= 0 ;
+ }
+
+ 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 ) )
+ && (i2 >= i1))
+ {
+ printf("ERREUR combinaison (%d, %d) -> CPUx=%lu CPUx2=%lu | GPUx=%lu GPUx2=%lu\n", i1, i2,
+ 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 ) ;
+ cpterr++;
+ }
+ cpt++ ;
+ }
+ }
+ printf("VERIF COMBINAISONS LIGNES --> %d ERREURS / %d COMBINAISONS\n", cpterr, cpt) ;
+ fin verif */
+
+
+ grid = dim3(div, 1, 1) ;
+ calcul_critere_permutations_verticales<<< grid, divpow2, CFI(divpow2)*sizeof(double2) >>>(d_contribs_cols, bs, j1, j2, H, L, sigX, sigX2, d_miniblocs) ;
+
+ /* verif CPU
+ cpt = 0 ;
+ cpterr = 0 ;
+ double2 * h_miniblocs_cpu = new double2[ div ] ;
+ cudaMemcpy( h_miniblocs, d_miniblocs, div*sizeof(double2), cudaMemcpyDeviceToHost) ;
+ cudaMemcpy( h_stats_snake, *d_stats_snake, 6*sizeof(int64), cudaMemcpyDeviceToHost) ;
+ for (int lb=0 ; lb < div ; lb++){
+ h_miniblocs_cpu[lb].x = DBL_MAX ;
+ for (int lh=lb ; lh < div ; lh++){
+ 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 )
+ {
+ 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) ;
+ h_miniblocs_cpu[lb].y = (double)lh ;
+ }
+ }
+ 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 ) )
+ {
+ printf("ERREUR MINIMUM BLOC LIGNE %d -> CPU=%lf en i2=%d | GPU=%lf en i2=%d\n", lb,
+ h_miniblocs_cpu[ lb ].x, (int)h_miniblocs_cpu[ lb ].y, h_miniblocs[ lb ].x, (int)h_miniblocs[ lb ].y ) ;
+ cpterr++;
+ }
+ cpt++ ;
+ }
+ printf("VERIF MINIMA PAR BLOC --> %d ERREURS / %d BLOCS\n", cpterr, cpt) ;
+ */
+ /* fin verif */
+
+ /*
+ * determination minimum absolu
+ * a conserver sur CPU
+ */
+ cudaMemcpy( h_miniblocs, d_miniblocs, div*sizeof(double2), cudaMemcpyDeviceToHost) ;
+ crit_mini = h_miniblocs[0].x ;
+ int i1=0 ;
+ id_mini=0 ;
+ for (i1=1 ; i1 < div ; i1++){
+ if (h_miniblocs[i1].x < crit_mini) {
+ crit_mini = h_miniblocs[i1].x ;
+ id_mini = i1 ;
+ }
+ }
+
+ i1 = bs * id_mini ;
+ int i2 = (int)(bs*h_miniblocs[ id_mini ].y) ;
+
+ toc(chrono, "CALCUL RECTANGLE");
+
+ printf("pas = %d lignes -- critere mini =%f , positions i1=%d i2=%d\n", bs, h_miniblocs[ id_mini ].x, i1, i2);
+ /*fin test snake rectangle initial optimal*/
+ cudaFree(d_contribs_part) ;
+ cudaFree(d_contribs_cols) ;
+ cudaFree(d_miniblocs) ;
+
+ tic(&chrono, NULL);
+ //genere_snake_rectangle_4nodes_gpu<<< 1, 1>>>(*d_snake, 140, H, L) ;
+ //genere_snake_bande_gpu<<<1,1>>>(*d_snake, pas*id_mini, (int)(pas*h_miniblocs[ id_mini ].y), H);
+ genere_snake_rectangle<<<1,1>>>(*d_snake, i1, i2, j1, j2);
+
+
