CUDA: blocks synchronisation

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <cstdlib>

const int REPEAT_TEST           = 10;
const int TRANSF_ITER           = 100;
const int BLOCKS []             = {  2,   4,   8, 16, 32};                  // проверки на соразмерность масивов нет!
const int THREADS_PER_BLOCK []  = {512, 256, 128, 64, 32};
const int THREADS_PER_BLOCK_max = 512;


__device__ volatile unsigned int count; //volatile
__device__ void do_while_count_not_eq(int val)
{
    /*
    // Цикл задержки в виде ассемблерной вставки. Точно будет в коде, не требует "volatile"-переменной count.
    asm("{\n\t"
        "$my_while_label: \n\t"
        "   .reg .u32      r_count;            \n\t"
        "   .reg .pred     p;                  \n\t"
        "   ld.global.u32  r_count, [count];   \n\t"
        "   setp.ne.u32    p, r_count, %0;     \n\t"
        "@p bra            $my_while_label;    \n\t"
        "}\n\t"
        : : "r"(val));
    */

    // Цикл задержки просто на Си. Компилятор nvcc "выкинет" его, если переменная count не определена с использованием "volatile".
    while (count != val);
}
/* Фунция начальной инициализации флага-счетчика: */
__device__ void InitSyncWholeDevice(const int index)
{
    if (index == 0)                                                         // Первый поток в grid`е запишет нулем (начальным значением) флаг-счетчик блоков.
        count = 0;

    if (threadIdx.x == 0)                                                   // Первый поток каждого block`а будет ждать, пока флаг-счетчик дейcтвительно станет нулем.
        do_while_count_not_eq(0);

    __syncthreads();                                                        // Заставляем остальные потоки каждого block`а ждать, пока первые не выйдут из цикла.
}                                                                           // Все, флаг-аккумулятор записан. Все потоки на device более-менее идут вровень.
/* Фунция синхронизации потоков на device: */
__device__ void SyncWholeDevice()
{
    unsigned int oldc;                                                      // Переменная под значение счетчика до инкремента.
    __threadfence();                                                        // Каждый поток ждет, пока записанное им в gmem и smem, станет видно всему grid`у. __syncthreads отвечает только за видимость для block.
    if (threadIdx.x == 0)                                                   // Первые потоки каждого block`а атомарным образом инкрементируют (каждый по разу) флаг-аккумулятор.
    {
        oldc = atomicInc((unsigned int *)&count, gridDim.x-1);              // В oldc кладется значение count до "+1".
        __threadfence();                                                    // Пусть поток подождет, пока его инкремент "дойдет" до ячейки в gmem.
        if (oldc != (gridDim.x-1))                                          // Если это последний блок (остальные уже инкрементировали count и ждут за счет цикла ниже), то и незачем ему считывать count. Так как предварительно убедились, что его инкремент записан в gmem.
            do_while_count_not_eq(0);                                       // Каждый первый поток в каждом block`е зациклен, пока все блоки не закончат работу.
    }
    __syncthreads();                                                        // Заставляем потоки в каждом блоке ждать, пока первые не выйдут из цикла.
}


__device__ void SingleTransform(const float *src, float *dest)
{
    float a = 0;

    /* Цикл вычисления элемента dest[]. В двух вариантах. */
    // Вариант №1: без использования разделяемой памяти:
        for(int j = 0; j < gridDim.x*blockDim.x; ++j)
            a += src[j];
    // Вариант №2: с использованием разделяемой памяти:
        /*__shared__ float smemBuff [THREADS_PER_BLOCK_max];
        for(int jb = 0; jb < gridDim.x; ++jb)
        {
            __syncthreads();                                                // Необходимо дождться, пока smemBuff[] больше не будет нужен никакому потоку из блока.
            smemBuff[threadIdx.x] = src[threadIdx.x + jb*blockDim.x];
            __syncthreads();                                                // Необходимо дождться, пока в smemBuff[] все потоки блока запишут значения.
            for(int jt = 0; jt < blockDim.x; ++jt)
                a += smemBuff[jt];
        }*/

    a /= (gridDim.x*blockDim.x);
    dest[threadIdx.x + blockIdx.x*blockDim.x] = a;
}

__global__ void Kernel_SKL(float *X, float *P)
{
    InitSyncWholeDevice(threadIdx.x + blockIdx.x*blockDim.x);
    for (int i = 0; i < TRANSF_ITER; ++i)
    {
        SingleTransform(X, P);
        SyncWholeDevice();
        SingleTransform(P, X);
        SyncWholeDevice();
    }
}

__global__ void Kernel_MKL(float *src, float *dest)
{
    SingleTransform(src, dest);
}


int main(void)
{
    int NUM_OF_ELEMS;
    float *X_h;
    float *X_d, *P_d;
    float *res[REPEAT_TEST+1];

    float elapsedTime_MKL[2];
    float elapsedTime_SKL[2];

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);

    FILE *fo = fopen("out.txt", "w");

    for(int jnb = 0; jnb < sizeof(BLOCKS)/sizeof(*BLOCKS); ++jnb)
    for(int jnt = 0; jnt < sizeof(THREADS_PER_BLOCK)/sizeof(*THREADS_PER_BLOCK); ++jnt)
    {
        printf("BLOCKS = %3i, THREADS_PER_BLOCK = %3i:\n", BLOCKS[jnb], THREADS_PER_BLOCK[jnt]);
        NUM_OF_ELEMS = BLOCKS[jnb] * THREADS_PER_BLOCK[jnt];

        X_h = new float [NUM_OF_ELEMS];
        cudaMalloc(&X_d, NUM_OF_ELEMS*sizeof(float));
        cudaMalloc(&P_d, NUM_OF_ELEMS*sizeof(float));

        srand(time(NULL));
        for(int j = 0; j < NUM_OF_ELEMS; ++j)
            X_h[j] = float(rand()) / RAND_MAX; // как вариант: float(j+1)/NUM_OF_ELEMS;

        for(int r = 0; r <= REPEAT_TEST; ++r)
            res[r] = new float [NUM_OF_ELEMS];

        elapsedTime_MKL[0] = elapsedTime_SKL[0] = 0;
        for(int r = 1; r <= REPEAT_TEST; ++r)
        {
            printf("  r = %i/%i\r", r, REPEAT_TEST);
            fflush(stdout);

            cudaMemcpy(X_d, X_h, NUM_OF_ELEMS*sizeof(float), cudaMemcpyHostToDevice);
            cudaEventRecord(start, 0);
            for (int i = 0; i < TRANSF_ITER; ++i)
            {
                Kernel_MKL <<<BLOCKS[jnb], THREADS_PER_BLOCK[jnt]>>> (X_d, P_d);
                Kernel_MKL <<<BLOCKS[jnb], THREADS_PER_BLOCK[jnt]>>> (P_d, X_d);
            }
            cudaEventRecord(stop, 0);
            cudaEventSynchronize(stop);
            cudaEventElapsedTime(&elapsedTime_MKL[1], start, stop);
            elapsedTime_MKL[0] += elapsedTime_MKL[1];
            cudaMemcpy(res[0], X_d, NUM_OF_ELEMS*sizeof(float), cudaMemcpyDeviceToHost);    // Дабы не усложнять пример, лишние несколько раз скопируем одно и то же в res[0].

            cudaMemcpy(X_d, X_h, NUM_OF_ELEMS*sizeof(float), cudaMemcpyHostToDevice);
            cudaEventRecord(start, 0);
            Kernel_SKL <<<BLOCKS[jnb], THREADS_PER_BLOCK[jnt]>>> (X_d, P_d);
            cudaEventRecord(stop, 0);
            cudaEventSynchronize(stop);
            cudaEventElapsedTime(&elapsedTime_SKL[1], start, stop);
            elapsedTime_SKL[0] += elapsedTime_SKL[1];
            cudaMemcpy(res[r], X_d, NUM_OF_ELEMS*sizeof(float), cudaMemcpyDeviceToHost);
        }
        printf("  Elapsed Time\n");
        printf("    multi kerenel launch:  %20.14f\n", elapsedTime_MKL[0]/REPEAT_TEST);
        printf("    single kerenel launch: %20.14f\n", elapsedTime_SKL[0]/REPEAT_TEST);
        fprintf(fo, "%i  %i  %20.14f   %20.14f\n", jnb, jnt, elapsedTime_MKL[0]/REPEAT_TEST, elapsedTime_SKL[0]/REPEAT_TEST);

        // Сравним между собой все j-е элементы из каждой r-ой реализации. Выявит численную ошибку за счет использования неактуальных значений элементов массивов.
        for(int j = 0; j < NUM_OF_ELEMS; ++j)
        {
            int r;
            float e = res[0][j];
            for(r = 1; r < REPEAT_TEST; ++r)
                if (e != res[r][j])
                {
                    printf("    Numeric error of SKL-method has been detected!\n");
                    break;
                }
            if (e != res[r][j])
                break;
        }

        cudaFree(X_d);
        cudaFree(P_d);
        delete[] X_h;
        for(int r = 0; r < (REPEAT_TEST+1); ++r)
            delete[] res[r];
    }

    fclose(fo);

    return 0;
}