CUDA reduction

reduce.cu

// This software contains source code provided by NVIDIA Corporation.
// http://docs.nvidia.com/cuda/eula/index.html#nvidia-cuda-samples-end-user-license-agreement

#include <cuda_runtime.h>

__global__ void reduce0(int *g_idata, int *g_odata, unsigned int n)
{
    // shared memoryは、うまく使えば速度が速い。(ここでは、うまく使ってない。後で直す)
    // 同じブロック内のスレッド間で共有される。
    extern __shared__ int sdata[];

    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    sdata[tid] = (i < n) ? g_idata[i] : 0;
    __syncthreads();

    // 例えばtidが0..7の場合、
    // sdata[0] += sdata[1]; sdata[2] += sdata[3]; sdata[4] += sdata[5]; sdata[6] += sdata[7];
    // sdata[0] += sdata[2]; sdata[4] += sdata[6];
    // sdata[0] += sdata[4];
    // のように、sdata[0]に集められる。
    for (unsigned int s=1; s<blockDim.x; s*=2) {
        if (tid % (2*s) == 0) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}

// reduce0 の
// if (tid % (2*s) == 0) {
//     sdata[tid] += sdata[tid + s];
// }
// 部分を改善したい。
// CUDAでは、物理的には、ワープという単位でスレッドが実行され、各ワープはプログラムカウンタを共有する。
// また、1ワープは32スレッドである。
// つまり、
// if (tid % (2*s) == 0) {
//     sdata[tid] += sdata[tid + s];
// }
// で、ワープ内に tid % (2*s) == 0 を満たすtidを持つスレッドと、満たさないtidを持つスレッドが混在していた場合、
// 満たさないtidを持つスレッドは、ワープ内の他のスレッドがif文の中身を動かしている間、待機することになる。
// (これをdivergent branchと呼ぶ)
// reduce0 では、例えばs=1のとき、
// tid = 0: sdata[0] += sdata[1];
// tid = 1: 待機
// tid = 2: sdata[2] += sdata[3];
// tid = 3: 待機
// ...
// のようになっていたが、今回は、
// tid = 0: sdata[0] += sdata[1];
// tid = 1: sdata[2] += sdata[3];
// ...
// tid = blockDim.x / 2 - 1: sdata[blockDim.x - 2] += sdata[blockDim.x - 1];
// tid = blockDim.x / 2: 待機
// tid = blockDim.x / 2 + 1: 待機
// ...
// のようになる。この対策で、スレッド数(= blockDim.x)が十分に大きいときは、
// divergent branchが起こるワープの数を減らせる。
__global__ void reduce1(int *g_idata, int *g_odata, unsigned int n)
{
    extern __shared__ int sdata[];

    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    sdata[tid] = (i < n) ? g_idata[i] : 0;
    __syncthreads();

    for (unsigned int s=1; s<blockDim.x; s*=2) {
        int index = 2 * s * tid;
        if (index < blockDim.x) {
            sdata[index] += sdata[index + s];
        }
        __syncthreads();
    }
    if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}

// sharedメモリは、うまく使えばとても速い。
// うまく使うとは? => バンクコンフリクトが起こらないように使う。
// sharedメモリにアクセスする際、内部的には「バンク」に対して、欲しいメモリ番地をリクエストしている。
// メモリ番地ごとに、どのバンクにリクエストを出すかが決まっていて、
// ワープ内の複数のスレッドが、同時に、同じバンクに対して、異なるアドレス(*1)にアクセスしたとき、
// バンクの順番待ちが起こる。これをバンクコンフリクトと呼ぶ。
// なお、同じバンクに同時にアクセスしていても、同じアドレス(*1)に対するアクセスであれば、バンクコンフリクトは起こらない。
// (*1) 異なる/同じ「アドレス」というのは厳密には嘘で、NVIDIAのC Programming Guideでは32-bit wordという表現が使われている。
//      つまり、アドレスが異なっていても、同じ32-bit wordに収まっていれば、バンクコンフリクトは起こらない。
// バンクは全部で32個あり、32-bitごとに、バンク0, バンク1, ..., バンク31, バンク0, ...
// に対して、リクエストを出すことになっている。
//
// さて。reduce0, reduce1では。各スレッドは、sharedメモリに、
// 2個おきにアクセス、4個おきにアクセス、8個おきにアクセス... としている。
// いかにもバンクコンフリクトを起こしそうだ。
//
// バンクコンフリクトを起こさないアクセス方法
// => 少なくとも、2個飛ばしとかせずに、シーケンシャルにアクセスする分には大丈夫。
// バンクは32個あって、1ワープあたりのスレッド数が32個なわけだから。コンフリクトしない。
// (けど、double型みたいに64ビット境界とかだと、それでもバンクコンフリクトするかも)
// つまり、
// sdata[0] += sdata[128]; sdata[1] += sdata[129]; sdata[2] += sdata[130]; sdata[3] += sdata[131]; ...
// のようにしてみる。
__global__ void reduce2(int *g_idata, int *g_odata, unsigned int n)
{
    extern __shared__ int sdata[];

    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * blockDim.x + threadIdx.x;
    sdata[tid] = (i < n) ? g_idata[i] : 0;
    __syncthreads();

    for (unsigned int s=blockDim.x/2; s>0; s>>=1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}

// ループの最初から既に、半分のスレッドが寝ている。こいつ、何もしてないんじゃね?
// メモリをsharedに移す段階で、1回足し算させましょう。
__global__ void reduce3(int *g_idata, int *g_odata, unsigned int n)
{
    extern __shared__ int sdata[];

    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
    int mysum = (i < n) ? g_idata[i] : 0;
    if (i + blockDim.x < n) mysum += g_idata[i + blockDim.x];
    sdata[tid] = mysum;
    __syncthreads();

    for (unsigned int s=blockDim.x/2; s>0; s>>=1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    if (tid == 0) g_odata[blockIdx.x] = sdata[0];
}

// Kepler以降では、ワープ内の別スレッドから変数を読む、shuffle operationというものがある。
// 1ワープに収まって以降は、これを使えば速くなる。
// https://devblogs.nvidia.com/parallelforall/faster-parallel-reductions-kepler/ によると、
// sharedメモリを使うよりも命令数が少なくて済むから速いらしい。
// shuffle operationが使えないアーキテクチャではループアンロールでもしとけ、とNVIDIAのサンプルにはある。
// ちょっと古めのNVIDIAの資料には、同じワープ内だとどうせ全部一緒に動くから__syncthreads不要だよ、とあったが、
// サンプルでは、__syncthreadsされている。このへん、何か変わったんだろうか?
// あと、__shfl_downしてるところではループアンロールしなくていいのかなぁ。
__global__ void reduce4(int *g_idata, int *g_odata, unsigned int n)
{
    extern __shared__ int sdata[];

#if !defined(__CUDA_ARCH__) || (__CUDA_ARCH__ >= 300)
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (blockDim.x * 2) + threadIdx.x;
    int mysum = (i < n) ? g_idata[i] : 0;
    if (i + blockDim.x < n) mysum += g_idata[i + blockDim.x];
    sdata[tid] = mysum;
    __syncthreads();

    for (unsigned int s=blockDim.x/2; s>32; s>>=1) {
        if (tid < s) {
            sdata[tid] = mysum = mysum + sdata[tid + s];
        }
        __syncthreads();
    }
    if (tid < 32) {
       if(blockDim.x >= 64) mysum += sdata[tid + 32];
        for (int offset = 32/2; offset>0; offset>>=1) {
            mysum += __shfl_down(mysum, offset);
        }
    }
    if (tid == 0) g_odata[blockIdx.x] = mysum;
#else
#error "__shfl_down requires CUDA arch >= 300."
#endif
}

// 全部アンロールする。blockDim.xは動的な変数だが、C++のテンプレートを使って、
// blockDim.xに相当する値を静的に求める。(呼び出し時にちょっとがんばらないといけない)
template <unsigned int BLOCKDIM>
__global__ void reduce5(int *g_idata, int *g_odata, unsigned int n)
{
    extern __shared__ int sdata[];

#if !defined(__CUDA_ARCH__) || (__CUDA_ARCH__ >= 300)
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (BLOCKDIM * 2) + threadIdx.x;
    int mysum = (i < n) ? g_idata[i] : 0;
    if (i + BLOCKDIM < n) mysum += g_idata[i + BLOCKDIM];
    sdata[tid] = mysum;
    __syncthreads();

    if (BLOCKDIM >= 1024 && tid < 512) sdata[tid] = mysum = mysum + sdata[tid + 512];
    __syncthreads();
    if (BLOCKDIM >= 512 && tid < 256) sdata[tid] = mysum = mysum + sdata[tid + 256];
    __syncthreads();
    if (BLOCKDIM >= 256 && tid < 128) sdata[tid] = mysum = mysum + sdata[tid + 128];
    __syncthreads();
    if (BLOCKDIM >= 128 && tid < 64) sdata[tid] = mysum = mysum + sdata[tid + 64];
    __syncthreads();
    if (tid < 32) {
        if (BLOCKDIM >= 64) mysum += sdata[tid + 32];
        for (int offset = 32/2; offset>0; offset>>=1) {
            mysum += __shfl_down(mysum, offset);
        }
    }
    if (tid == 0) g_odata[blockIdx.x] = mysum;
#else
#error "__shfl_down requires CUDA arch >= 300."
#endif
}

// ここでは、計算時間ではなく、計算コストの観点で、効率化を測りたい。
// 計算コストとは、各スレッド(あるいはプロセッサ)にかかる時間の合計で、これまでの計算では
// O(N) スレッドがそれぞれ O(log N) の計算時間を使うので、
// 計算コストは O(N log N) となる。
// これは、CPUでの計算コスト O(N) と比べて、よくない。
// ところで、計算時間について、 Brent's theorem という定理がある。
// Tp = T∞ + (T1 - T∞) / p
// ただし、ここで、
// Tp: p個のプロセッサで並列計算したときの計算時間
// p:  計算に使うプロセッサ数
// T∞: プロセッサがいくらでも好きなだけ並列に使えると仮定したときの計算時間
// T1: 1個のプロセッサで計算したときの計算時間
// なので、今回の場合、
// Tp = O(log N + (N - log N) / p)
// となる。
// ここで、p = N / log N とすると、計算コストは O(N) となる。それを示そう。
// p * Tp = O(p log N + (N - log N))
//        = O(N + (N - log N))
//        = O(N)
// よって、 p = N / log N のとき、計算コストは O(N) とできる。
// このような p で計算を行うには、各スレッドが、 O(log N) 個の要素を足し合わせればいい。
// ...のような話が、NVIDIAの資料に載っていて、そこから、なぜか「最大ブロック数を64にする」
// という話になっている。
// まぁ、それでも O(N) になるはずだし、CUDAコア数の関係上、ブロック数をもっと増やしたところで、
// どうせ全部同時には動かないだろうし、いいんだろうけど。Brent's theoremの話は何だったんだ感はなくもない。
// なお、テンプレート引数の N_IS_POW2 は、nが2の累乗だと true となるようにしておく。
// これにより、コンパイル時に無駄な分岐が省かれる。
// が。これがtrueのとき、nによってはバグる気がする。例えばn=64のとき。
template <unsigned int BLOCKDIM, bool N_IS_POW2>
__global__ void reduce6(int *g_idata, int *g_odata, unsigned int n)
{
    extern __shared__ int sdata[];

#if !defined(__CUDA_ARCH__) || (__CUDA_ARCH__ >= 300)
    unsigned int tid = threadIdx.x;
    unsigned int i = blockIdx.x * (BLOCKDIM * 2) + threadIdx.x;
    unsigned int gridsize = BLOCKDIM * 2 * gridDim.x;
    int mysum = 0;

    while (i < n) {
        mysum += g_idata[i];
        if (/*N_IS_POW2 ||*/ i + BLOCKDIM < n) mysum += g_idata[i + BLOCKDIM];
        i += gridsize;
    }
    sdata[tid] = mysum;
    __syncthreads();

    if (BLOCKDIM >= 1024 && tid < 512) sdata[tid] = mysum = mysum + sdata[tid + 512];
    __syncthreads();
    if (BLOCKDIM >= 512 && tid < 256) sdata[tid] = mysum = mysum + sdata[tid + 256];
    __syncthreads();
    if (BLOCKDIM >= 256 && tid < 128) sdata[tid] = mysum = mysum + sdata[tid + 128];
    __syncthreads();
    if (BLOCKDIM >= 128 && tid < 64) sdata[tid] = mysum = mysum + sdata[tid + 64];
    __syncthreads();
    if (tid < 32) {
        if (BLOCKDIM >= 64) mysum += sdata[tid + 32];
        for (int offset = 32/2; offset>0; offset>>=1) {
            mysum += __shfl_down(mysum, offset);
        }
    }
    if (tid == 0) g_odata[blockIdx.x] = mysum;
#else
#error "__shfl_down requires CUDA arch >= 300."
#endif
}

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

#ifndef N
#define N ((1 << 28) + 1234)
#endif

template<class T>
bool device_assert_eq(T *pt, T v) {
    T t;
    cudaMemcpy(&t, pt, sizeof(T), cudaMemcpyDeviceToHost);
    assert(t == v);
}

void launch_reduce5(int griddim, int blockdim, int sharedsize, int* g_idata, int* g_odata, unsigned int n)
{
    switch(blockdim) {
    case 32:
        reduce5<32><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
        break;
    case 64:
        reduce5<64><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
        break;
    case 128:
        reduce5<128><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
        break;
    case 256:
        reduce5<256><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
        break;
    case 512:
        reduce5<512><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
        break;
    case 1024:
        reduce5<1024><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
        break;
    default:
        printf("Invalid thread size.\n");
    }
}

bool isPow2(unsigned int x)
{
    return ((x&(x-1))==0);
}

void launch_reduce6(int griddim, int blockdim, int sharedsize, int* g_idata, int* g_odata, unsigned int n)
{
    bool is_pow2 = isPow2(n);
    if (is_pow2) {
        switch(blockdim) {
        case 32:
            reduce6<32, true><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        case 64:
            reduce6<64, true><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        case 128:
            reduce6<128, true><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        case 256:
            reduce6<256, true><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        case 512:
            reduce6<512, true><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        case 1024:
            reduce6<1024, true><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        default:
            printf("Invalid thread size.\n");
        }
    }
    else {
        switch(blockdim) {
        case 32:
            reduce6<32, false><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        case 64:
            reduce6<64, false><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        case 128:
            reduce6<128, false><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        case 256:
            reduce6<256, false><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        case 512:
            reduce6<512, false><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        case 1024:
            reduce6<1024, false><<<griddim, blockdim, sharedsize>>>(g_idata, g_odata, n);
            break;
        default:
            printf("Invalid thread size.\n");
        }
    }
}

int main(void)
{
    printf("N:      %10d\n", N);
    int cpu_result = 0;
    int *arr_host = (int*)malloc(N * sizeof(int));
    for (int i=0; i<N; i++) arr_host[i] = 1;

    clock_t cl = clock();
    for (int i=0; i<N; i++) cpu_result += arr_host[i];
    assert(cpu_result = N);
    double time = (clock() - cl) / (double)CLOCKS_PER_SEC * 1000.0;
    printf("CPU:    %10.4f\n", time);

    int threads[] = { 32, 64, 128, 256, 512, 1024 };

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

    int *arr_dev;
    cudaMalloc((void**)&arr_dev, sizeof(int) * N);
    cudaMemcpy(arr_dev, arr_host, N*sizeof(int), cudaMemcpyHostToDevice);

    printf(
    "threads   %8d  %8d  %8d  %8d  %8d  %8d\n"
    "reduce0:  ", threads[0], threads[1], threads[2], threads[3], threads[4], threads[5]);
    for (int i=0; i<sizeof(threads) / sizeof(threads[0]); i++) {
        int *out1_dev, *out2_dev;
        int th = threads[i];
        int blocks = (N - 1) / th + 1;
        int shared_mem_size = 2 * th * sizeof(int);
        cudaMalloc((void**)&out1_dev, sizeof(int) * blocks);
        cudaMalloc((void**)&out2_dev, sizeof(int) * blocks);

        cudaEventRecord(start);
        int **in = &arr_dev, **out = &out1_dev;
        int n = N;
        while (blocks > 1) {
            reduce0<<<blocks, th, shared_mem_size>>>(*in, *out, n);
            cudaError_t err = cudaGetLastError();
            if (err != cudaSuccess) {
                puts(cudaGetErrorString(err));
            }
            if (*out == out1_dev) {
                out = &out2_dev; in = &out1_dev;
            }
            else {
                out = &out1_dev; in = &out2_dev;
            }
            n = blocks;
            blocks = (blocks - 1) / th + 1;
            cudaThreadSynchronize();
        }
        reduce0<<<blocks, th, shared_mem_size>>>(*in, *out, n);
        int result;
        cudaMemcpy(&result, *out, sizeof(int), cudaMemcpyDeviceToHost);
        //printf("%d\n", result);
        cudaEventRecord(stop);
        cudaEventSynchronize(stop);

        assert(result == N);
        float milliseconds = 0;
        cudaEventElapsedTime(&milliseconds, start, stop);
        printf("%8.4f  ", milliseconds);

        cudaFree(out1_dev);
        cudaFree(out2_dev);
    }
    printf("\n");

    printf("reduce1:  ");
    for (int i=0; i<sizeof(threads) / sizeof(threads[0]); i++) {
        int *out1_dev, *out2_dev;
        int th = threads[i];
        int blocks = (N - 1) / th + 1;
        int shared_mem_size = 2 * th * sizeof(int);
        cudaMalloc((void**)&out1_dev, sizeof(int) * blocks);
        cudaMalloc((void**)&out2_dev, sizeof(int) * blocks);

        cudaEventRecord(start);
        int **in = &arr_dev, **out = &out1_dev;
        int n = N;
        while (blocks > 1) {
            reduce1<<<blocks, th, shared_mem_size>>>(*in, *out, n);
            cudaError_t err = cudaGetLastError();
            if (err != cudaSuccess) {
                puts(cudaGetErrorString(err));
            }
            if (*out == out1_dev) {
                out = &out2_dev; in = &out1_dev;
            }
            else {
                out = &out1_dev; in = &out2_dev;
            }
            n = blocks;
            blocks = (blocks - 1) / th + 1;
            cudaThreadSynchronize();
        }
        reduce1<<<blocks, th, shared_mem_size>>>(*in, *out, n);
        int result;
        cudaMemcpy(&result, *out, sizeof(int), cudaMemcpyDeviceToHost);
        //printf("%d\n", result);
        cudaEventRecord(stop);
        cudaEventSynchronize(stop);

        assert(result == N);
        float milliseconds = 0;
        cudaEventElapsedTime(&milliseconds, start, stop);
        printf("%8.4f  ", milliseconds);

        cudaFree(out1_dev);
        cudaFree(out2_dev);
    }
    printf("\n");

    printf("reduce2:  ");
    for (int i=0; i<sizeof(threads) / sizeof(threads[0]); i++) {
        int *out1_dev, *out2_dev;
        int th = threads[i];
        int blocks = (N - 1) / th + 1;
        int shared_mem_size = 2 * th * sizeof(int);
        cudaMalloc((void**)&out1_dev, sizeof(int) * blocks);
        cudaMalloc((void**)&out2_dev, sizeof(int) * blocks);

        cudaEventRecord(start);
        int **in = &arr_dev, **out = &out1_dev;
        int n = N;
        while (blocks > 1) {
            reduce2<<<blocks, th, shared_mem_size>>>(*in, *out, n);
            cudaError_t err = cudaGetLastError();
            if (err != cudaSuccess) {
                puts(cudaGetErrorString(err));
            }
            if (*out == out1_dev) {
                out = &out2_dev; in = &out1_dev;
            }
            else {
                out = &out1_dev; in = &out2_dev;
            }
            n = blocks;
            blocks = (blocks - 1) / th + 1;
            cudaThreadSynchronize();
        }
        reduce2<<<blocks, th, shared_mem_size>>>(*in, *out, n);
        int result;
        cudaMemcpy(&result, *out, sizeof(int), cudaMemcpyDeviceToHost);
        //printf("%d\n", result);
        cudaEventRecord(stop);
        cudaEventSynchronize(stop);

        assert(result == N);
        float milliseconds = 0;
        cudaEventElapsedTime(&milliseconds, start, stop);
        printf("%8.4f  ", milliseconds);

        cudaFree(out1_dev);
        cudaFree(out2_dev);
    }
    printf("\n");

    printf("reduce3:  ");
    for (int i=0; i<sizeof(threads) / sizeof(threads[0]); i++) {
        int *out1_dev, *out2_dev;
        int th = threads[i];
        int blocks = (N - 1) / (2 * th) + 1;
        int shared_mem_size = 2 * th * sizeof(int);
        cudaMalloc((void**)&out1_dev, sizeof(int) * blocks);
        cudaMalloc((void**)&out2_dev, sizeof(int) * blocks);

        cudaEventRecord(start);
        int **in = &arr_dev, **out = &out1_dev;
        int n = N;
        while (blocks > 1) {
            reduce3<<<blocks, th, shared_mem_size>>>(*in, *out, n);
            cudaError_t err = cudaGetLastError();
            if (err != cudaSuccess) {
                puts(cudaGetErrorString(err));
            }
            if (*out == out1_dev) {
                out = &out2_dev; in = &out1_dev;
            }
            else {
                out = &out1_dev; in = &out2_dev;
            }
            n = blocks;
            blocks = (blocks - 1) / (2 * th) + 1;
            cudaThreadSynchronize();
        }
        reduce3<<<blocks, th, shared_mem_size>>>(*in, *out, n);
        int result;
        cudaMemcpy(&result, *out, sizeof(int), cudaMemcpyDeviceToHost);
        cudaEventRecord(stop);
        cudaEventSynchronize(stop);

        assert(result == N);
        float milliseconds = 0;
        cudaEventElapsedTime(&milliseconds, start, stop);
        printf("%8.4f  ", milliseconds);

        cudaFree(out1_dev);
        cudaFree(out2_dev);
    }
    printf("\n");

    printf("reduce4:  ");
    for (int i=0; i<sizeof(threads) / sizeof(threads[0]); i++) {
        int *out1_dev, *out2_dev;
        int th = threads[i];
        int blocks = (N - 1) / (2 * th) + 1;
        int shared_mem_size = 2 * th * sizeof(int);
        cudaMalloc((void**)&out1_dev, sizeof(int) * blocks);
        cudaMalloc((void**)&out2_dev, sizeof(int) * blocks);

        cudaEventRecord(start);
        int **in = &arr_dev, **out = &out1_dev;
        int n = N;
        while (blocks > 1) {
            reduce4<<<blocks, th, shared_mem_size>>>(*in, *out, n);
            cudaError_t err = cudaGetLastError();
            if (err != cudaSuccess) {
                puts(cudaGetErrorString(err));
            }
            if (*out == out1_dev) {
                out = &out2_dev; in = &out1_dev;
            }
            else {
                out = &out1_dev; in = &out2_dev;
            }
            n = blocks;
            blocks = (blocks - 1) / (2 * th) + 1;
            cudaThreadSynchronize();
        }
        reduce4<<<blocks, th, shared_mem_size>>>(*in, *out, n);
        int result;
        cudaMemcpy(&result, *out, sizeof(int), cudaMemcpyDeviceToHost);
        cudaEventRecord(stop);
        cudaEventSynchronize(stop);

        assert(result == N);
        float milliseconds = 0;
        cudaEventElapsedTime(&milliseconds, start, stop);
        printf("%8.4f  ", milliseconds);

        cudaFree(out1_dev);
        cudaFree(out2_dev);
    }
    printf("\n");

    printf("reduce5:  ");
    for (int i=0; i<sizeof(threads) / sizeof(threads[0]); i++) {
        int *out1_dev, *out2_dev;
        int th = threads[i];
        int blocks = (N - 1) / (2 * th) + 1;
        int shared_mem_size = 2 * th * sizeof(int);
        cudaMalloc((void**)&out1_dev, sizeof(int) * blocks);
        cudaMalloc((void**)&out2_dev, sizeof(int) * blocks);

        cudaEventRecord(start);
        int **in = &arr_dev, **out = &out1_dev;
        int n = N;
        while (blocks > 1) {
            launch_reduce5(blocks, th, shared_mem_size, *in, *out, n);
            cudaError_t err = cudaGetLastError();
            if (err != cudaSuccess) {
                puts(cudaGetErrorString(err));
            }
            if (*out == out1_dev) {
                out = &out2_dev; in = &out1_dev;
            }
            else {
                out = &out1_dev; in = &out2_dev;
            }
            n = blocks;
            blocks = (blocks - 1) / (2 * th) + 1;
            cudaThreadSynchronize();
        }
        launch_reduce5(blocks, th, shared_mem_size, *in, *out, n);
        int result;
        cudaMemcpy(&result, *out, sizeof(int), cudaMemcpyDeviceToHost);
        cudaEventRecord(stop);
        cudaEventSynchronize(stop);

        assert(result == N);
        float milliseconds = 0;
        cudaEventElapsedTime(&milliseconds, start, stop);
        printf("%8.4f  ", milliseconds);

        cudaFree(out1_dev);
        cudaFree(out2_dev);
    }
    printf("\n");

    printf("reduce6:  ");
    for (int i=0; i<sizeof(threads) / sizeof(threads[0]); i++) {
        int *out1_dev, *out2_dev;
        int th = threads[i];
        int blocks = (N - 1) / (2 * th) + 1;
        int shared_mem_size = 2 * th * sizeof(int);
        cudaMalloc((void**)&out1_dev, sizeof(int) * blocks);
        cudaMalloc((void**)&out2_dev, sizeof(int) * blocks);

        cudaEventRecord(start);
        int **in = &arr_dev, **out = &out1_dev;
        int n = N;
        while (blocks > 1) {
            int b = blocks;
            if (b > 64) b = 64;
            launch_reduce6(b, th, shared_mem_size, *in, *out, n);
            cudaError_t err = cudaGetLastError();
            if (err != cudaSuccess) {
                puts(cudaGetErrorString(err));
            }
            if (*out == out1_dev) {
                out = &out2_dev; in = &out1_dev;
            }
            else {
                out = &out1_dev; in = &out2_dev;
            }
            n = b;
            blocks = (b - 1) / (2 * th) + 1;
            cudaThreadSynchronize();
        }
        launch_reduce6(blocks, th, shared_mem_size, *in, *out, n);
        int result;
        cudaMemcpy(&result, *out, sizeof(int), cudaMemcpyDeviceToHost);
        cudaEventRecord(stop);
        cudaEventSynchronize(stop);

        assert(result == N);
        float milliseconds = 0;
        cudaEventElapsedTime(&milliseconds, start, stop);
        printf("%8.4f  ", milliseconds);

        cudaFree(out1_dev);
        cudaFree(out2_dev);
    }
    printf("\n");
}