0

我目前正在开展一个项目，在该项目中我正在展开减少的最后一个扭曲。我已经完成了上面的代码；但是，一些修改是通过猜测完成的，我想解释一下原因。我写的代码只有函数kernel4

// in is input array, out is where to store result, n is number of elements from in
// T is a float (32bit)
__global__ void kernel4(T *in, T *out, unsigned int n)

这是一个缩减算法，其余代码已经提供。

代码：

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

#include "timer.h"
#include "cuda_utils.h"

typedef float T;

#define N_ (8 * 1024 * 1024)
#define MAX_THREADS 256
#define MAX_BLOCKS 64

#define MIN(x,y) ((x < y) ? x : y)
#define tid threadIdx.x 
#define bid blockIdx.x 
#define bdim blockDim.x
#define warp_size 32

unsigned int nextPow2( unsigned int x ) {
    --x;
    x |= x >> 1;
    x |= x >> 2;
    x |= x >> 4;
    x |= x >> 8;
    x |= x >> 16;
    return ++x;
}

void getNumBlocksAndThreads(int whichKernel, int n, int maxBlocks, int maxThreads, int &blocks, int &threads)
{
    if (whichKernel < 3) {
        threads = (n < maxThreads) ? nextPow2(n) : maxThreads;
        blocks = (n + threads - 1) / threads;
    } else {
        threads = (n < maxThreads*2) ? nextPow2((n + 1)/ 2) : maxThreads;
        blocks = (n + (threads * 2 - 1)) / (threads * 2);
    }
    if (whichKernel == 5)
        blocks = MIN(maxBlocks, blocks);
}

T reduce_cpu(T *data, int n) {
    T sum = data[0];
    T c = (T) 0.0;
    for (int i = 1; i < n; i++)
    {
        T y = data[i] - c;
        T t = sum + y;
        c = (t - sum) - y;
        sum = t;
    } 
    return sum;
}

__global__ void
kernel4(T *in, T *out, unsigned int n)
{
    __shared__ volatile T d[MAX_THREADS];

    unsigned int i = bid * bdim + tid;

    n >>= 1;
    d[tid] = (i < n) ? in[i] + in[i+n] : 0;
    __syncthreads ();

    for(unsigned int s = bdim >> 1; s > warp_size; s >>= 1) {
        if(tid < s)
            d[tid] += d[tid + s];
        __syncthreads ();
    }

    if (tid < warp_size) {
        if (n > 64) d[tid] += d[tid + 32];
        if (n > 32) d[tid] += d[tid + 16];
        d[tid] += d[tid + 8];
        d[tid] += d[tid + 4];
        d[tid] += d[tid + 2];
        d[tid] += d[tid + 1];
    }

    if(tid == 0)
        out[bid] = d[0];
}

int main(int argc, char** argv)
{
    T *h_idata, h_odata, h_cpu;
    T *d_idata, *d_odata;   

    struct stopwatch_t* timer = NULL;
    long double t_kernel_4, t_cpu;

    int whichKernel = 4, threads, blocks, N, i;
    if(argc > 1) {
        N = atoi (argv[1]);
        printf("N: %d\n", N);
    } else {
        N = N_;
        printf("N: %d\n", N);
    }

    getNumBlocksAndThreads (whichKernel, N, MAX_BLOCKS, MAX_THREADS, blocks, threads);

    stopwatch_init ();
    timer = stopwatch_create ();

    h_idata = (T*) malloc (N * sizeof (T));
    CUDA_CHECK_ERROR (cudaMalloc (&d_idata, N * sizeof (T)));
    CUDA_CHECK_ERROR (cudaMalloc (&d_odata, blocks * sizeof (T)));

    srand48(time(NULL));
    for(i = 0; i < N; i++)
        h_idata[i] = drand48() / 100000;
    CUDA_CHECK_ERROR (cudaMemcpy (d_idata, h_idata, N * sizeof (T), cudaMemcpyHostToDevice));

    dim3 gb(blocks, 1, 1);
    dim3 tb(threads, 1, 1);

    kernel4 <<<gb, tb>>> (d_idata, d_odata, N);
    cudaThreadSynchronize ();

    stopwatch_start (timer);

    kernel4 <<<gb, tb>>> (d_idata, d_odata, N);
    int s = blocks;
    while(s > 1) {
        threads = 0;
        blocks = 0;
        getNumBlocksAndThreads (whichKernel, s, MAX_BLOCKS, MAX_THREADS, blocks, threads);

        dim3 gb(blocks, 1, 1);
        dim3 tb(threads, 1, 1);

        kernel4 <<<gb, tb>>> (d_odata, d_odata, s);
        s = (s + threads * 2 - 1) / (threads * 2);
    }
    cudaThreadSynchronize ();

    t_kernel_4 = stopwatch_stop (timer);
    fprintf (stdout, "Time to execute unrolled GPU reduction kernel: %Lg secs\n", t_kernel_4);

    double bw = (N * sizeof(T)) / (t_kernel_4 * 1e9);   // total bits / time
    fprintf (stdout, "Effective bandwidth: %.2lf GB/s\n", bw);
    CUDA_CHECK_ERROR (cudaMemcpy (&h_odata, d_odata, sizeof (T), cudaMemcpyDeviceToHost));

    stopwatch_start (timer);
    h_cpu = reduce_cpu (h_idata, N);
    t_cpu = stopwatch_stop (timer);
    fprintf (stdout, "Time to execute naive CPU reduction: %Lg secs\n", t_cpu);

    if(abs (h_odata - h_cpu) > 1e-5)
        fprintf(stderr, "FAILURE: GPU: %f  CPU: %f\n", h_odata, h_cpu);
    else
        printf("SUCCESS: GPU: %f  CPU: %f\n", h_odata, h_cpu);
    return 0;
}

我的第一个问题是：在声明时

__shared__ volatile T d[MAX_THREADS];

我想验证我对 volatile 的理解。Volatile 防止编译器错误地优化我的代码，并承诺加载/存储是通过缓存完成的，而不仅仅是寄存器（如果有错误请纠正我）。对于归约，如果部分归约和仍存储在寄存器中，为什么会出现问题？

我的第二个问题是：在进行实际的翘曲减少时

    if (tid < warp_size) { // Final log2(32) = 5 strides
        if (n > 64) d[tid] += d[tid + 32];
        if (n > 32) d[tid] += d[tid + 16];
        d[tid] += d[tid + 8];
        d[tid] += d[tid + 4];
        d[tid] += d[tid + 2];
        d[tid] += d[tid + 1];
    }

在没有 (n > 64) 和 (n > 32) 条件的情况下，归约和将产生不正确的结果。我得到的结果是：

FAILURE: GPU: 41.966557  CPU: 41.946209

经过 5 次试验，GPU 缩减始终产生 0.0204 的误差。我很谨慎地认为这是一个浮点运算错误。

老实说，我的老师的助手建议进行此更改以添加 (n > 64) 和 (n > 32) 条件，但没有解释为什么它会修复代码。

由于我的试验中的 n 超过 64，为什么这个条件会改变结果。我很难追溯问题，因为我不能像在 CPU 中那样使用打印功能。

4

2 回答 2

4

在我们解决您的两个问题之前，让我们从一些前言评论开始：

  • 我鼓励你阅读 NVIDIA 的规范化简教程
  • 像这样写的缩减做了几个假设，其中一个是块大小是 2 的幂（为了“正确性”）。
  • 您的代码在最终缩减阶段使用经纱同步编程。你似乎知道你在做什么，所以我不会提供详细的描述，但这里肯定与理解相关。如果需要，您可以 google 并获取说明。它与下面的讨论相关，但我不会在每种情况下都指出它的相关性。

好的，现在你的问题：

我想验证我对 volatile 的理解。Volatile 防止编译器错误地优化我的代码，并承诺加载/存储是通过缓存完成的，而不仅仅是寄存器（如果有错误请纠正我）。对于归约，如果部分归约和仍存储在寄存器中，为什么会出现问题？

关于 的定义volatile，我建议您参考CUDA 编程指南。我已经看到了将其称为防止寄存器优化或防止加载和存储重新排序的摘要描述。我更喜欢前者，并将其用作工作定义。

基本思想是volatile强制对该变量的任何引用（读取或写入）实际转到内存子系统。我的意思是它将执行读取或写入，并且不会尝试使用先前加载到寄存器中的值。如果没有这个限定符，编译器可以自由地（例如）从实际内存位置加载一次值，然后在寄存器中维护该值（以及对它的任何更新），只要它认为合适。编译器这样做是着眼于性能。（顺便说一句，请注意您在这里使用了“缓存”这个词。我会避免在这里使用这种用法。共享内存在它和处理器加载/存储机制之间没有插入缓存。）

如果没有volatile这种类型的扭曲同步编码，如果我们允许编译器“优化”（即维护）中间值到寄存器中，我们就会遇到问题。这主要是由于线程间通信引起的。为了清楚地了解原因，让我们看一下最终归约的最后两个步骤：

    d[tid] += d[tid + 2];
    d[tid] += d[tid + 1];

tid让我们只考虑值为 0-1的线程。在倒数第二步中，线程 0 将获取该d[2]值并将其添加到该d[0]值中，而线程 1 将获取该d[3]值并将其添加到该d[1]值中。此时，如果我们不使用volatile，编译器没有义务将d[1]线程 1 累积的值写回共享内存。允许将其保存在寄存器中。因此，d[1]在共享内存中看到的值不是“最新的”。

现在让我们进入最后一步。在这一步中，线程 0从共享内存中读取该d[1]值并将其添加到该值中。但是没有，我们在上一步中看到 的共享内存内容不再准确。OTOH，如果我们使用，那么上一步中对共享内存的写入实际上会发生，并且在最后一步中，线程 0 将在读取时获取正确的值。CUDA 线程是一个独立的模型。我的意思是，一个线程不能直接访问属于另一个线程的寄存器中包含的值。因此，warp 级别的线程间通信通常将通过共享内存或通过 warp-shuffle 操作来完成。d[0]volatiled[1]volatiled[1]

__syncthreads()具有类似的行为：它强制将所有像这样的寄存器优化值写入内存，以便它们对块中的其他线程“可见”。因此，更复杂的优化将是仅volatile当归约从基于循环驱动__syncthreads()的归约切换到最终的扭曲同步归约时才切换到合格的指针。您可以在我在此答案开头链接的教程幻灯片中看到一个示例。

另外，这种扭曲同步编程在 CUDA 9 中（更正式地）被弃用。相反，您应该使用合作组。

在没有 (n > 64) 和 (n > 32) 条件的情况下，归约和将产生不正确的结果。

主要使用这些条件是因为代码被设计为对于任何具有 2 次幂大小的块配置都是“正确的”。如果我们假设块大小（每个块的线程数）是 2 的幂，并且大于 64，例如它必须是 128 或更大。您的n变量以块大小开始，但随后乘以 2：

n >>= 1;

因此，如果我们要保证这行代码的正确性：

d[tid] += d[tid + 32];

那么我们应该只在线程块大小为 64（至少）时应用该操作，这就像说n大于 64：

    if (n > 64) d[tid] += d[tid + 32];

关于这个问题，声称无论if (n > 64)是否包含，发布的代码的行为都会有所不同。这样做的原因是发布的代码包含一个循环，该循环在减少过程中重新计算线程数和块数：

int s = blocks;
while(s > 1) {
    threads = 0;
    blocks = 0;
    getNumBlocksAndThreads (whichKernel, s, MAX_BLOCKS, MAX_THREADS, blocks, threads);

此循环最终导致块大小小于 128，这意味着省略 if 条件会导致损坏。（在此循环期间只需打印出threads变量）。

关于这个：

我很难追溯问题，因为我不能像在 CPU 中那样使用打印功能。

我不确定那里有什么问题。 printf应该在内核代码中工作。

于 2018-03-08T03:29:49.053 回答
0

根据此答案，共享变量不能将初始化作为其声明的一部分。因此，如果 n < 64，我们将一些随机共享内存数组数据添加到总和中，这种情况下会出错。

于 2019-05-21T04:38:26.810 回答