Does cub use global memory to do prefix scan?

[cyk2018]

I have read an article and there are three different methods to get the prefix sum of an array, 1 is using global memory, 2 is using shared memory, 3 is using thrust::inclusive_scan, every way with 20 similar test.
I find that the first test will be always too long time than other 19, such as almost 3 or 4 times(only in method 1 and method 3). But in method 2(shared memory), the time are almost same in 20 tests.
I know why the first test is long than other, maybe cache. But I don't know why the cub are longer than shared memory, seems We use global memory algorithm?
I find code in thrust and there is no memory information, It call function in cub, so I ask in this field, waiting for your reply sincerely.
Thanks

[pauleonix]

It might be helpful to share said article so we know what you mean. Assuming enough data I'm sure that CUB's implementation uses both shared memory for local scans and global memory for communication between blocks to update the local scans with the final term from previous blocks. As far as I know, Single-pass Parallel Prefix Scan with Decoupled Look-back describes what CUB is doing. Could you also provide more information on that slow test case? How many elements are scanned? What GPU are you using?

[cyk2018]

Sure, firstly thanks for your reply.
This article is written by Chinese, and the link is here. The input data is an array which length is 1e6, and the elements are all 1.23. the GPU I used is A800. And I can give you the code.

#include <stdio.h>
#include "error.cuh"
#include <thrust/scan.h>
#include <thrust/execution_policy.h>

#ifdef USE_DP
    typedef double real;
#else
    typedef float real;
#endif


const int N = 1000000;
const int M = sizeof(int) * N;
const int NUM_REAPEATS = 20;
const int BLOCK_SIZE = 1024;
const int GRID_SIZE = (N - 1) / BLOCK_SIZE + 1;


__global__ void globalMemScan(real *d_x, real *d_y) {
    real *x = d_x + blockDim.x * blockIdx.x;
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    real y = 0.0;
    if (n < N) {
        for (int offset = 1; offset < blockDim.x; offset <<= 1) {
            if (threadIdx.x >= offset) y = x[threadIdx.x] + x[threadIdx.x - offset];
            __syncthreads();
            if (threadIdx.x >= offset) x[threadIdx.x] = y;
        }
        if (threadIdx.x == blockDim.x - 1) d_y[blockIdx.x] = x[threadIdx.x];
    } 
}

__global__ void addBaseValue(real *d_x, real *d_y) {
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    real y = blockIdx.x > 0 ? d_y[blockIdx.x - 1] : 0.0;
    if (n < N) {
        d_x[n] += y;
    } 
}

__global__ void sharedMemScan(real *d_x, real *d_y) {
    extern __shared__ real s_x[];
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    s_x[threadIdx.x] = n < N ? d_x[n] : 0.0;
    __syncthreads();
    real y = 0.0;
    if (n < N) {
        for (int offset = 1; offset < blockDim.x; offset <<= 1) {
            if (threadIdx.x >= offset) y = s_x[threadIdx.x] + s_x[threadIdx.x - offset];
            __syncthreads();
            if (threadIdx.x >= offset) s_x[threadIdx.x] = y;
        }
        d_x[n] = s_x[threadIdx.x];
        if (threadIdx.x == blockDim.x - 1) d_y[blockIdx.x] = s_x[threadIdx.x];
    } 
}


void scan(real *d_x, real *d_y, real *d_z, const int method) {
    switch (method)
    {
    case 0:
        globalMemScan<<<GRID_SIZE, BLOCK_SIZE>>>(d_x, d_y);
        globalMemScan<<<1, GRID_SIZE>>>(d_y, d_z);
        addBaseValue<<<GRID_SIZE, BLOCK_SIZE>>>(d_x, d_y);
        break;
    case 1:
        sharedMemScan<<<GRID_SIZE, BLOCK_SIZE, sizeof(real) * BLOCK_SIZE>>>(d_x, d_y);
        sharedMemScan<<<1, GRID_SIZE, sizeof(real) * GRID_SIZE>>>(d_y, d_z);
        addBaseValue<<<GRID_SIZE, BLOCK_SIZE>>>(d_x, d_y);
        break;
    case 2:
        thrust::inclusive_scan(thrust::device, d_x, d_x + N, d_x);
        break;
    
    default:
        break;
    }
}

void timing(const real *h_x, real *d_x, real *d_y, real *d_z, real *h_ret, const int method) {

    float tSum = 0.0;
    float t2Sum = 0.0;
    float tSumVersion2 = 0.0;
    float t2SumVersion2 = 0.0;
    for (int i=0; i<=NUM_REAPEATS; ++i) {
        CHECK(cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice));
        cudaEvent_t start, stop;
        CHECK(cudaEventCreate(&start));
        CHECK(cudaEventCreate(&stop));
        CHECK(cudaEventRecord(start));
        cudaEventQuery(start);

        scan(d_x, d_y, d_z, method);

        CHECK(cudaEventRecord(stop));
        CHECK(cudaEventSynchronize(stop));
        float elapsedTime;
        CHECK(cudaEventElapsedTime(&elapsedTime, start, stop));
        if(i == 0) {
            // do nothing
        } else {
            tSumVersion2 += elapsedTime;
            t2SumVersion2 += elapsedTime * elapsedTime;
        }

        if(i == NUM_REAPEATS) {
            // do nothing
        } else {
            tSum += elapsedTime;
            t2Sum += elapsedTime * elapsedTime;
        }
        
        printf("%g\t", elapsedTime);
        CHECK(cudaEventDestroy(start));
        CHECK(cudaEventDestroy(stop));
    }
    printf("\n======version1======\n");
    printf("\n%g\t", tSum);
    float tAVG = tSum / NUM_REAPEATS;
    float tERR = sqrt(t2Sum / NUM_REAPEATS - tAVG * tAVG);
    printf("Time = %g +- %g ms.\n", tAVG, tERR);
    printf("\n======version2======\n");
    printf("\n%g\t", tSumVersion2);
    float tAVGVersion2 = tSumVersion2 / NUM_REAPEATS;
    float tERRVersion2 = sqrt(t2SumVersion2 / NUM_REAPEATS - tAVGVersion2 * tAVGVersion2);
    printf("Time = %g +- %g ms.\n", tAVGVersion2, tERRVersion2);
    CHECK(cudaMemcpy(h_ret, d_x, M, cudaMemcpyDeviceToHost));
}

int main() {
    real *h_x = new real[N];
    real *h_y = new real[N];    
    real *h_ret = new real[N];
    for (int i=0; i<N; i++) h_x[i] = 1.23;
    real *d_x, *d_y, *d_z;
    CHECK(cudaMalloc((void **)&d_x, M));
    CHECK(cudaMalloc((void **)&d_y, sizeof(real) * GRID_SIZE));
    CHECK(cudaMalloc((void **)&d_z, sizeof(real)));
    
    printf("using global mem:\n");
    timing(h_x, d_x, d_y, d_z, h_ret, 0);
    for (int i = N - 10; i < N; i++) printf("%f  ", h_ret[i]);
    printf("\n");
    printf("using shared mem:\n");
    timing(h_x, d_x, d_y, d_z, h_ret, 1);
    for (int i = N - 10; i < N; i++) printf("%f  ", h_ret[i]);
    printf("\n");
    printf("using thrust lib:\n");
    timing(h_x, d_x, d_y, d_z, h_ret, 2);
    for (int i = N - 10; i < N; i++) printf("%f  ", h_ret[i]);
    printf("\n");

    CHECK(cudaFree(d_x));
    CHECK(cudaFree(d_y));
    CHECK(cudaFree(d_z));
    delete[] h_x;
    delete[] h_y;
    delete[] h_ret;
}

The error.cuh is a helper program:

#pragma once
#include <stdio.h>
#include <cublas_v2.h>

static const char *_cudaGetErrorEnum(cublasStatus_t error)
{
    switch (error)
    {
    case CUBLAS_STATUS_SUCCESS:
        return "CUBLAS_STATUS_SUCCESS";

    case CUBLAS_STATUS_NOT_INITIALIZED:
        return "CUBLAS_STATUS_NOT_INITIALIZED";

    case CUBLAS_STATUS_ALLOC_FAILED:
        return "CUBLAS_STATUS_ALLOC_FAILED";

    case CUBLAS_STATUS_INVALID_VALUE:
        return "CUBLAS_STATUS_INVALID_VALUE";

    case CUBLAS_STATUS_ARCH_MISMATCH:
        return "CUBLAS_STATUS_ARCH_MISMATCH";

    case CUBLAS_STATUS_MAPPING_ERROR:
        return "CUBLAS_STATUS_MAPPING_ERROR";

    case CUBLAS_STATUS_EXECUTION_FAILED:
        return "CUBLAS_STATUS_EXECUTION_FAILED";

    case CUBLAS_STATUS_INTERNAL_ERROR:
        return "CUBLAS_STATUS_INTERNAL_ERROR";

    case CUBLAS_STATUS_NOT_SUPPORTED:
        return "CUBLAS_STATUS_NOT_SUPPORTED";

    case CUBLAS_STATUS_LICENSE_ERROR:
        return "CUBLAS_STATUS_LICENSE_ERROR";
    }
    return "<unknown>";
}

#define CHECK_CUDA_ERROR(call)                             \
    do                                                     \
    {                                                      \
        const cudaError_t errorCode = call;                \
        if (errorCode != cudaSuccess)                      \
        {                                                  \
            printf("CUDA Error:\n");                       \
            printf("    File:   %s\n", __FILE__);          \
            printf("    Line:   %d\n", __LINE__);          \
            printf("    Error code:     %d\n", errorCode); \
            printf("    Error text:     %s\n",             \
                   cudaGetErrorString(errorCode));         \
            exit(1);                                       \
        }                                                  \
    } while (0)

#define CHECK_CUBLAS_STATUS(call)                            \
    do                                                       \
    {                                                        \
        const cublasStatus_t statusCode = call;              \
        if (statusCode != CUBLAS_STATUS_SUCCESS)             \
        {                                                    \
            printf("CUDA Error:\n");                         \
            printf("    File:   %s\n", __FILE__);            \
            printf("    Line:   %d\n", __LINE__);            \
            printf("    Status code:     %d\n", statusCode); \
            printf("    Error text:     %s\n",               \
                   _cudaGetErrorEnum(statusCode));           \
            exit(1);                                         \
        }                                                    \
    } while (0)

You can also get the code here

I can give you a result image:

version1 here is the first 20 avarage time; version2 is the 1-21 result. and the upper array is every time they use.

[elstehle]

As far as I know, Single-pass Parallel Prefix Scan with Decoupled Look-back describes what CUB is doing.

@pauleonix is exactly right. CUB is using the single-pass prefix scan to minimize incurred memory traffic. That is, to communicate partial and inclusive prefix scan results of each tile (I am referring to "a tile", as the items that one thread block processes). To compute the results within one tile, you can assume that the CUB implementation is at least as sophisticated as the shared memory variant you were referring to above.

In general, I would strongly advise to use the CUB algorithms. There's a lot of thought that went into the design of these algorithms and you don't have to maintain your own implementations and making sure they will be kept up to date with recent hardware advancements.

Now, to the performance numbers you are seeing: it is worth pointing out that thrust requires allocating some temporary storage to store some intermediary result, similar to the d_y and d_z memory buffers that your global and shared memory implementations are using. For thrust, however, you are including the time for the memory allocation that happens within the call to thrust::inclusive_scan.

I have attached a code snippet based on the code you provided that uses cub::DeviceScan. The snippet pre-allocates the temporary memory required by the cub algorithm, similar to your other benchmarks.

Also make sure to compile your code targeting the GPU architecture(s) you'll be running on. E.g., for V100 / SM 70

nvcc -lineinfo -gencode arch=compute_70,code=sm_70  my_scan_tests.cu

Code snippet using cub::DeviceScan

#include <stdio.h>
#include <cublas_v2.h>

static const char *_cudaGetErrorEnum(cublasStatus_t error)
{
    switch (error)
    {
    case CUBLAS_STATUS_SUCCESS:
        return "CUBLAS_STATUS_SUCCESS";

    case CUBLAS_STATUS_NOT_INITIALIZED:
        return "CUBLAS_STATUS_NOT_INITIALIZED";

    case CUBLAS_STATUS_ALLOC_FAILED:
        return "CUBLAS_STATUS_ALLOC_FAILED";

    case CUBLAS_STATUS_INVALID_VALUE:
        return "CUBLAS_STATUS_INVALID_VALUE";

    case CUBLAS_STATUS_ARCH_MISMATCH:
        return "CUBLAS_STATUS_ARCH_MISMATCH";

    case CUBLAS_STATUS_MAPPING_ERROR:
        return "CUBLAS_STATUS_MAPPING_ERROR";

    case CUBLAS_STATUS_EXECUTION_FAILED:
        return "CUBLAS_STATUS_EXECUTION_FAILED";

    case CUBLAS_STATUS_INTERNAL_ERROR:
        return "CUBLAS_STATUS_INTERNAL_ERROR";

    case CUBLAS_STATUS_NOT_SUPPORTED:
        return "CUBLAS_STATUS_NOT_SUPPORTED";

    case CUBLAS_STATUS_LICENSE_ERROR:
        return "CUBLAS_STATUS_LICENSE_ERROR";
    }
    return "<unknown>";
}

#define CHECK_CUDA_ERROR(call)                             \
    do                                                     \
    {                                                      \
        const cudaError_t errorCode = call;                \
        if (errorCode != cudaSuccess)                      \
        {                                                  \
            printf("CUDA Error:\n");                       \
            printf("    File:   %s\n", __FILE__);          \
            printf("    Line:   %d\n", __LINE__);          \
            printf("    Error code:     %d\n", errorCode); \
            printf("    Error text:     %s\n",             \
                   cudaGetErrorString(errorCode));         \
            exit(1);                                       \
        }                                                  \
    } while (0)

#define CHECK_CUBLAS_STATUS(call)                            \
    do                                                       \
    {                                                        \
        const cublasStatus_t statusCode = call;              \
        if (statusCode != CUBLAS_STATUS_SUCCESS)             \
        {                                                    \
            printf("CUDA Error:\n");                         \
            printf("    File:   %s\n", __FILE__);            \
            printf("    Line:   %d\n", __LINE__);            \
            printf("    Status code:     %d\n", statusCode); \
            printf("    Error text:     %s\n",               \
                   _cudaGetErrorEnum(statusCode));           \
            exit(1);                                         \
        }                                                    \
    } while (0)

#include <stdio.h>
#include <cub/cub.cuh>
#include <thrust/scan.h>
#include <thrust/execution_policy.h>

#ifdef USE_DP
    typedef double real;
#else
    typedef float real;
#endif


const int N = 1000000;
const int M = sizeof(int) * N;
const int NUM_REAPEATS = 20;
const int BLOCK_SIZE = 1024;
const int GRID_SIZE = (N - 1) / BLOCK_SIZE + 1;


__global__ void globalMemScan(real *d_x, real *d_y) {
    real *x = d_x + blockDim.x * blockIdx.x;
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    real y = 0.0;
    if (n < N) {
        for (int offset = 1; offset < blockDim.x; offset <<= 1) {
            if (threadIdx.x >= offset) y = x[threadIdx.x] + x[threadIdx.x - offset];
            __syncthreads();
            if (threadIdx.x >= offset) x[threadIdx.x] = y;
        }
        if (threadIdx.x == blockDim.x - 1) d_y[blockIdx.x] = x[threadIdx.x];
    } 
}

__global__ void addBaseValue(real *d_x, real *d_y) {
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    real y = blockIdx.x > 0 ? d_y[blockIdx.x - 1] : 0.0;
    if (n < N) {
        d_x[n] += y;
    } 
}

__global__ void sharedMemScan(real *d_x, real *d_y) {
    extern __shared__ real s_x[];
    const int n = blockDim.x * blockIdx.x + threadIdx.x;
    s_x[threadIdx.x] = n < N ? d_x[n] : 0.0;
    __syncthreads();
    real y = 0.0;
    if (n < N) {
        for (int offset = 1; offset < blockDim.x; offset <<= 1) {
            if (threadIdx.x >= offset) y = s_x[threadIdx.x] + s_x[threadIdx.x - offset];
            __syncthreads();
            if (threadIdx.x >= offset) s_x[threadIdx.x] = y;
        }
        d_x[n] = s_x[threadIdx.x];
        if (threadIdx.x == blockDim.x - 1) d_y[blockIdx.x] = s_x[threadIdx.x];
    } 
}


void scan(real *d_x, real *d_y, real *d_z, const int method) {
    switch (method)
    {
    case 0:
        globalMemScan<<<GRID_SIZE, BLOCK_SIZE>>>(d_x, d_y);
        globalMemScan<<<1, GRID_SIZE>>>(d_y, d_z);
        addBaseValue<<<GRID_SIZE, BLOCK_SIZE>>>(d_x, d_y);
        break;
    case 1:
        sharedMemScan<<<GRID_SIZE, BLOCK_SIZE, sizeof(real) * BLOCK_SIZE>>>(d_x, d_y);
        sharedMemScan<<<1, GRID_SIZE, sizeof(real) * GRID_SIZE>>>(d_y, d_z);
        addBaseValue<<<GRID_SIZE, BLOCK_SIZE>>>(d_x, d_y);
        break;
    
    default:
        break;
    }
}

void cub_scan(real *d_x, void *temp_storage, std::size_t &temp_storage_bytes) {
    cub::DeviceScan::InclusiveSum(temp_storage, temp_storage_bytes, d_x, d_x, N);
}

void timing(const real *h_x, real *d_x, real *d_y, real *d_z, real *h_ret, const int method) {

    float tSum = 0.0;
    float t2Sum = 0.0;
    float tSumVersion2 = 0.0;
    float t2SumVersion2 = 0.0;

    // Query CUB::DeviceScan temporary storage requirements
    void *temp_storage = nullptr;
    std::size_t temp_storage_bytes{};
    cub_scan(d_x, temp_storage, temp_storage_bytes);
    CHECK_CUDA_ERROR(cudaMalloc((void **)&temp_storage, temp_storage_bytes));

    for (int i=0; i<=NUM_REAPEATS; ++i) {
        CHECK_CUDA_ERROR(cudaMemcpy(d_x, h_x, M, cudaMemcpyHostToDevice));
        cudaEvent_t start, stop;
        CHECK_CUDA_ERROR(cudaEventCreate(&start));
        CHECK_CUDA_ERROR(cudaEventCreate(&stop));
        CHECK_CUDA_ERROR(cudaEventRecord(start));
        cudaEventQuery(start);

        if(method == 2){
            cub_scan(d_x, temp_storage, temp_storage_bytes);
        }else{
        scan(d_x, d_y, d_z, method);
        }

        CHECK_CUDA_ERROR(cudaEventRecord(stop));
        CHECK_CUDA_ERROR(cudaEventSynchronize(stop));
        float elapsedTime;
        CHECK_CUDA_ERROR(cudaEventElapsedTime(&elapsedTime, start, stop));
        if(i == 0) {
            // do nothing
        } else {
            tSumVersion2 += elapsedTime;
            t2SumVersion2 += elapsedTime * elapsedTime;
        }

        if(i == NUM_REAPEATS) {
            // do nothing
        } else {
            tSum += elapsedTime;
            t2Sum += elapsedTime * elapsedTime;
        }
        
        printf("%g\t", elapsedTime);
        CHECK_CUDA_ERROR(cudaEventDestroy(start));
        CHECK_CUDA_ERROR(cudaEventDestroy(stop));
    }
    printf("\n======version1======\n");
    printf("\n%g\t", tSum);
    float tAVG = tSum / NUM_REAPEATS;
    float tERR = sqrt(t2Sum / NUM_REAPEATS - tAVG * tAVG);
    printf("Time = %g +- %g ms.\n", tAVG, tERR);
    printf("\n======version2======\n");
    printf("\n%g\t", tSumVersion2);
    float tAVGVersion2 = tSumVersion2 / NUM_REAPEATS;
    float tERRVersion2 = sqrt(t2SumVersion2 / NUM_REAPEATS - tAVGVersion2 * tAVGVersion2);
    printf("Time = %g +- %g ms.\n", tAVGVersion2, tERRVersion2);
    CHECK_CUDA_ERROR(cudaMemcpy(h_ret, d_x, M, cudaMemcpyDeviceToHost));
}

int main() {
    real *h_x = new real[N];
    real *h_y = new real[N];    
    real *h_ret = new real[N];
    for (int i=0; i<N; i++) h_x[i] = 1.23;
    real *d_x, *d_y, *d_z;
    CHECK_CUDA_ERROR(cudaMalloc((void **)&d_x, M));
    CHECK_CUDA_ERROR(cudaMalloc((void **)&d_y, sizeof(real) * GRID_SIZE));
    CHECK_CUDA_ERROR(cudaMalloc((void **)&d_z, sizeof(real)));
    
    printf("using global mem:\n");
    timing(h_x, d_x, d_y, d_z, h_ret, 0);
    for (int i = N - 10; i < N; i++) printf("%f  ", h_ret[i]);
    printf("\n");
    printf("using shared mem:\n");
    timing(h_x, d_x, d_y, d_z, h_ret, 1);
    for (int i = N - 10; i < N; i++) printf("%f  ", h_ret[i]);
    printf("\n");
    printf("using thrust lib:\n");
    timing(h_x, d_x, d_y, d_z, h_ret, 2);
    for (int i = N - 10; i < N; i++) printf("%f  ", h_ret[i]);
    printf("\n");

    CHECK_CUDA_ERROR(cudaFree(d_x));
    CHECK_CUDA_ERROR(cudaFree(d_y));
    CHECK_CUDA_ERROR(cudaFree(d_z));
    delete[] h_x;
    delete[] h_y;
    delete[] h_ret;
}

[cyk2018]

Great! I use your code in my GPU with arch=86. And the result are here:

Thanks for you and @pauleonix 's reply.
I have an other question that I want to test the time of operators in different implementation such as pytorch and tensorflow, which tool or framework I can use. I hope you can give me some great advice about this question.

[jrhemstad]

@cyk2018 we recently gave a talk that provided a deep dive into CUB's decoupled lookback scan algorithm and explicitly compared it to other common scan algorithms. I highly recommend you check this out to get a deeper understanding: https://www.youtube.com/watch?v=VLdm3bV4bKo

[cyk2018]

Thanks I will watch it.

[cyk2018]

I have watch your great video, I found that the definition of sol is wonderful, It give us a destination about how to optimize our algorithm, But I am confused at how can we calculate sol for a big program and are there some framework or tools for this?

I also found there is no your slide in the CUDA_MODE github repo, can you upload that?
Thanks for you two!