     //Initialize superaccs
     for (uint i = 0; i < BIN_COUNT; i++)
         l_workingBase[i * WARP_COUNT] = 0;
-    __syncthreads();
+    __syncthreads(); //syncs all threads in a block (but not across blocks)
 
     //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 );
    *status = 0;
    if( d_errorV[0] ) *status = 1;
}

        d_PartialSuperaccs[gid * BIN_COUNT * MERGE_SIZE + lid] = sum;
    }

    __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: 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 < blocks; i++)
            sum += d_PartialSuperaccs[i * BIN_COUNT * MERGE_SIZE + lid];

        d_superacc[lid] = sum;
