         d_PartialSuperaccs[gid * BIN_COUNT * MERGE_SIZE + lid] = sum;
     }
 
-    __syncthreads();
+    __syncthreads(); //syncs all threads in a block (but not across blocks)
     if (lid == 0) { //every block normalize its summed superacc
         int imin = IMIN, imax = IMAX;
         Normalize(&d_PartialSuperaccs[gid * BIN_COUNT * MERGE_SIZE], imin, imax);
     }
-
-    //MW: don't we need a global synchronization here??
-    __syncthreads();
+}
+//MW: we need a global synchronization here!!
+//one block of threads with at least 39 threads
+template<uint MERGE_SIZE>
+__global__
+void ExDOTCompleteFinal(
+     int64_t *d_PartialSuperaccs,
+     int64_t *d_superacc
+) {
+    uint lid = threadIdx.x;
+    uint gid = blockIdx.x;
+    uint blocks = gpu::PARTIAL_SUPERACCS_COUNT/gpu::MERGE_SUPERACCS_SIZE;
     if ((lid < BIN_COUNT) && (gid == 0)) {
         int64_t sum = 0;
 
-        for(uint i = 0; i < gridDim.x; i++)
+        for(uint i = 0; i < blocks; i++)
             sum += d_PartialSuperaccs[i * BIN_COUNT * MERGE_SIZE + lid];
 
         d_superacc[lid] = sum;

    //Initialize superaccs
    for (uint i = 0; i < BIN_COUNT; i++)
        l_workingBase[i * WARP_COUNT] = 0;
    __syncthreads();

    //Read data from global memory and scatter it to sub-superaccs
    double a[NBFPE] = {0.0};

    bool *d_error = thrust::raw_pointer_cast( d_errorV.data());
    gpu::ExDOT<NBFPE, gpu::WARP_COUNT><<<gpu::PARTIAL_SUPERACCS_COUNT, gpu::WORKGROUP_SIZE>>>( d_PartialSuperaccs, x1_ptr, x2_ptr, x3_ptr,size,d_error);
    gpu::ExDOTComplete<gpu::MERGE_SUPERACCS_SIZE><<<gpu::PARTIAL_SUPERACCS_COUNT/gpu::MERGE_SUPERACCS_SIZE, gpu::MERGE_WORKGROUP_SIZE>>>( d_PartialSuperaccs, d_superacc );
    //# blocks, # threads per block
    gpu::ExDOTCompleteFinal<gpu::MERGE_SUPERACCS_SIZE><<<1, 64>>>( d_PartialSuperaccs, d_superacc );
    *status = 0;
    if( d_errorV[0] ) *status = 1;
}
