openmp - 使用 OpenMP 并行化高斯模糊算法

Question

我试图使用 OpenMP 并行化高斯模糊函数， 但我是 OpenMP 的新手，当我尝试并行化两个 for 循环时(我不认为有任何变量需要为每个线程私有(private))，它最终 运行速度比以前更慢，并且输出也不同。那么我做错了什么吗？我应该怎样做才能让它运行得更快？

void gaussian_blur(float **src, float **dst, int w, int h, float sigma)
{
    int x, y, i;
    int ksize = (int)(sigma * 2.f * 4.f + 1) | 1;
    int halfk = ksize / 2;
    float scale = -0.5f/(sigma*sigma);
    float sum = 0.f;
    float *kernel, *ringbuf;
    int xmax = w - halfk;
    int ymax = h - halfk;

    // if sigma too small, just copy src to dst
    if (ksize <= 1)
    {
        for (y = 0; y < h; y++)
            for (x = 0; x < w; x++)
                dst[y][x] = src[y][x];
        return;
    }

    // create Gaussian kernel
    kernel = malloc(ksize * sizeof(float));
    ringbuf = malloc(ksize * sizeof(float));

    #pragma omp parallel for reduction(+ : sum)
    for (i = 0; i < ksize; i++)
    {
        float x = (float)(i - halfk);
        float t = expf(scale * x * x);
        kernel[i] = t;
        sum += t;
    }

    scale = 1.f / sum;
    #pragma omp parallel for
    for (i = 0; i < ksize; i++)
        kernel[i] *= scale;

    // blur each row
    #pragma omp parallel for // this is the for loop I parallelized but ended up with wrong output and running slower
    for (y = 0; y < h; y++)
    {
        int x1;
        int bufi0 = ksize-1;
        float tmp = src[y][0];

        for (x1 = 0; x1 < halfk  ; x1++) ringbuf[x1] = tmp;
        for (; x1 < ksize-1; x1++) ringbuf[x1] = src[y][x1-halfk];

        for (x1 = 0; x1 < w; x1++)
        {
            if(x1 < xmax)
                ringbuf[bufi0++] = src[y][x1+halfk];
            else
                ringbuf[bufi0++] = src[y][w-1];
            if (bufi0 == ksize) bufi0 = 0;
            dst[y][x1] = convolve(kernel, ringbuf, ksize, bufi0);
        }
    }

    // blur each column
    #pragma omp parallel for // this is the for loop I parallelized but ended up with wrong output and running slower
    for (x = 0; x < w; x++)
    {
        int y1;
        int bufi0 = ksize-1;
        float tmp = dst[0][x];
        for (y1 = 0; y1 < halfk  ; y1++) ringbuf[y1] = tmp;
        for (     ; y1 < ksize-1; y1++) ringbuf[y1] = dst[y1-halfk][x];

        for (y1 = 0; y1 < h; y1++)
        {
            if(y1 < ymax)
                ringbuf[bufi0++] = dst[y1+halfk][x];
            else
                ringbuf[bufi0++] = dst[h-1][x];

            if (bufi0 == ksize) bufi0 = 0;
            dst[y1][x] = convolve(kernel, ringbuf, ksize, bufi0);
        }
    }

    // clean up
    free(kernel);
    free(ringbuf);
}

Answer 1

除了需要正确识别私有(private)数据和共享数据之外，您还可以采取一些措施来加速程序。

第一步，您应该删除任何不必要的并发。例如，ksize 平均有多大？如果它少于几百个元素，那么使用 OpenMP 来执行计算内核然后对其进行规范化这样的简单操作是完全没有意义的:

#pragma omp parallel for reduction(+ : sum)
for (i = 0; i < ksize; i++)
{
    float x = (float)(i - halfk);
    float t = expf(scale * x * x);
    kernel[i] = t;
    sum += t;
}

scale = 1.f / sum;
#pragma omp parallel for
for (i = 0; i < ksize; i++)
    kernel[i] *= scale;

在典型的现代 CPU 上，引导并行区域比在单核上计算需要更多的周期。此外，在现代 CPU 上，这些循环可以展开和矢量化，并且您可以在单核上获得高达 8 倍的提升。如果内核太小，那么除了 OpenMP 开销之外，过多的错误共享也会导致速度减慢。您必须确保每个线程获取 16 个元素的精确倍数(64 字节缓存行大小/sizeof(float))来处理，以防止错误共享。

您还必须确保线程不会共享列模糊部分中的缓存行。

// blur each column
#pragma omp parallel for
for (x = 0; x < w; x++)
{
    ...
    for (y1 = 0; y1 < h; y1++)
    {
        ...
        dst[y1][x] = convolve(kernel, ringbuf, ksize, bufi0);
    }
}

由于这里的访问模式，您必须确保每个线程获取的列 block 是 16 的倍数，否则将会出现 16*y1 像素的边框重叠区域由每两个连续线程共享，其中会发生过多的错误共享。如果您不能保证 w 能被 16 整除，那么您可以为每个线程指定一个 y 方向的起始偏移量，例如最内层循环变为:

int tid = omp_get_thread_num();

for (y1 = 2*tid; y1 < h; y1++)
{
   ...
}
for (y1 = 0; y1 < 2*tid; y1++)
{
   ...
}

乘数2是任意的。这个想法是让下一个线程比当前线程提前几行，这样两个线程就不会在任何时刻同时处理同一行。您还可以使用加法和模算术来计算 y1，即

for (y2 = 0; y2 < h; y2++)
{
    y1 = (y2 + 2*tid) % h;
    ...
}

但这通常比将循环分成两部分要慢。

还要注意您的数据大小。最后一级缓存(LLC)具有非常高但仍然有限的带宽。如果数据无法放入每个内核的私有(private)缓存中，那么循环向量化等编译器优化可能会给 LLC 带来非常高的压力。如果数据不适合 LLC，事情会变得更加丑陋，因此必须访问主存储器。

如果您不知道什么是虚假共享，Dr.Dobb 中有一篇文章对此进行了解释 here .