CUDA 曼德布罗套装

Question

这是一个顺序 Mandelbrot Set 实现。

 void mandelbrot(PGMData *I)
{
    float x0,y0,x,y,xtemp;
    int i,j;
    int color;
    int iter;
    int MAX_ITER=1000;  
    for(i=0; i<I->height; i++)
        for(j=0; j<I->width; j++)
        {
            x0 = (float)j/I->width*(float)3.5-(float)2.5; 
            y0 = (float)i/I->height*(float)2.0-(float)1.0;
            x = 0;
            y = 0;
            iter = 0;
            while((x*x-y*y <= 4) && (iter < MAX_ITER))
            { 
                xtemp = x*x-y*y+x0;
                y = 2*x*y+y0;
                x = xtemp;
                iter++;
            }
            color = (int)(iter/(float)MAX_ITER*(float)I->max_gray);
            I->image[i*I->width+j] = I->max_gray-color;
        }
}

我想使用 CUDA 并行化它，但我似乎误解了一些东西，现在我被卡住了。我试过在互联网上搜索，但没有找到真正好的结果。

内核：

__global__ void calc(int *pos)
{
    int row= blockIdx.y * blockDim.y + threadIdx.y;  // WIDTH
    int col = blockIdx.x * blockDim.x + threadIdx.x;  // HEIGHT
    int idx = row * WIDTH + col;

    if(col > WIDTH || row > HEIGHT || idx > N) return;

    float x0 = (float)row/WIDTH*(float)3.5-(float)2.5;
    float y0 = (float)col/HEIGHT*(float)2.0-(float)1.0; 

    int x = 0, y = 0, iter = 0, xtemp = 0;
    while((x*x-y*y <= 4) && (iter < MAX_ITER))
    { 
        xtemp = x*x-y*y+x0;
        y = 2*x*y+y0;
        x = xtemp;
        iter++;
    }
    int color = 255 - (int)(iter/(float)MAX_ITER*(float)255);
    __syncthreads();
    pos[idx] = color;//color;// - color;

}

内核是这样启动的：

dim3 block_size(16, 16);
dim3 grid_size((N)/block_size.x, (int) N / block_size.y);
calc<<<grid_size,block_size>>>(d_pgmData);

这里是常量：

#define HEIGHT 512
#define WIDTH 512   
#define N (HEIGHT*WIDTH)

整个GPU功能

void mandelbrotGPU(PGMData *I)
{
    int *pos = (int *)malloc(HEIGHT*WIDTH*sizeof(int));
    int *d_pgmData;

    cudaMalloc((void **)&d_pgmData, sizeof(int)*WIDTH*HEIGHT);


    cudaMemcpy(d_pgmData, pos ,HEIGHT*WIDTH*sizeof(int) ,cudaMemcpyHostToDevice);

    dim3 block_size(16, 16);
    dim3 grid_size((N)/block_size.x, (int) N / block_size.y);
    calc<<<grid_size,block_size>>>(d_pgmData);

    cudaMemcpy(pos,d_pgmData,HEIGHT*WIDTH*sizeof(int) ,cudaMemcpyDeviceToHost);
    cudaFree(d_pgmData);
    I->image = pos;
}

问题是：要么返回垃圾，要么驱动程序崩溃。我真的很感激一些建议，因为我真的被困住了。

Answer 1

只是一些想法：

1. 不需要__syncthreads()。块中的线程不相互通信。
2. 无需为 I_WIDTH 和 I_HEIGHT 创建设备内存。您只需将它们作为值传递（而不是作为指针或引用）。您确实需要pos 的设备内存。
3. 您需要检查所有 CUDA 函数的返回值（例如，cudaMalloc）并确保一切正常。
4. 当您启动内核时，您的程序可能会在您的 GPU 完成之前返回。在某些情况下，您需要明确等待完成；您可以在启动后调用cudaDeviceSynchronize() 来执行此操作。在您的情况下，您不必这样做，因为您的 CUDA memcpy 将等到内核完成。

Answer 2

这肯定是不正确的：

    dim3 grid_size((N)/block_size.x, (int) N / block_size.y);

这会导致内核中的越界访问。您想总共启动 WIDTH x HEIGHT 线程，一个用于图像中的每个像素。相反，您正在启动 N/16 x N/16 个线程。

您的内核中似乎有一个线程检查行（应该可以防止错误线程的越界访问），但它的表述不正确：

if(col > WIDTH || row > HEIGHT || idx > N) return;

例如，这允许idx = N 通过线程检查，但是当内核的最后一行写入时，这不是有效的内存位置：

pos[idx] = color;

您可以使用以下方法修复此线程检查：

if(col >= WIDTH || row >= HEIGHT || idx >= N) return;

其他几个cmets：

• 你没有做任何事情proper cuda error checking
• 如果您使用cuda-memcheck 运行代码，您可以在代码中看到错误

Answer 3

这是您的代码的工作版本（使用 OpenCV）：

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include <iostream>
#include <opencv2/core/core.hpp>
#include <opencv2/highgui/highgui.hpp>

using namespace cv;
using namespace std;

#define HEIGHT 512 // must be multiple of block_size.y
#define WIDTH 512 // must be multiple of block_size.x
#define MAX_ITER 10000

void mandelbrotGPU(char*);
__global__ void calc(char* image_buffer);

#define cudaAssertSuccess(ans) { _cudaAssertSuccess((ans), __FILE__, __LINE__); }
inline void _cudaAssertSuccess(cudaError_t code, char *file, int line)
{
  if (code != cudaSuccess)  {
    fprintf(stderr,"_cudaAssertSuccess: %s %s %d\n", cudaGetErrorString(code), file, line);
    exit(code);
  }
}

int main(int argc, char** argv)
{
  IplImage* image_output = cvCreateImage(cvSize(WIDTH, HEIGHT), IPL_DEPTH_8U, 1);
  mandelbrotGPU(image_output->imageData);
  cvShowImage("GPU", image_output);
  waitKey(0);
  cvReleaseImage(&image_output);
}

void mandelbrotGPU(char* image_buffer)
{
  char* d_image_buffer;
  cudaAssertSuccess(cudaMalloc(&d_image_buffer, WIDTH * HEIGHT));
  dim3 block_size(16, 16);
  dim3 grid_size(WIDTH / block_size.x, HEIGHT / block_size.y);
  calc<<<grid_size, block_size>>>(d_image_buffer);
  cudaAssertSuccess(cudaPeekAtLastError());
  cudaAssertSuccess(cudaDeviceSynchronize());
  cudaAssertSuccess(cudaMemcpy(image_buffer, d_image_buffer, HEIGHT * WIDTH, cudaMemcpyDeviceToHost));
  cudaAssertSuccess(cudaFree(d_image_buffer));
}

__global__ void calc(char* image_buffer)
{
  int row = blockIdx.y * blockDim.y + threadIdx.y;  // WIDTH
  int col = blockIdx.x * blockDim.x + threadIdx.x;  // HEIGHT
  int idx = row * WIDTH + col;
  if(col >= WIDTH || row >= HEIGHT) return;

  float x0 = ((float)col / WIDTH) * 3.5f - 2.5f;
  float y0 = ((float)row / HEIGHT) * 3.5f - 1.75f;

  float x = 0.0f;
  float y = 0.0f;
  int iter = 0;
  float xtemp;
  while((x * x + y * y <= 4.0f) && (iter < MAX_ITER))
  { 
    xtemp = x * x - y * y + x0;
    y = 2.0f * x * y + y0;
    x = xtemp;
    iter++;
  }

  int color = iter * 5;
  if (color >= 256) color = 0;
  image_buffer[idx] = color;
}

输出：

最重要的变化：

• 删除了__syncthreads();。该算法不使用其他线程生成的数据，因此无需同步线程。
• 删除了将主机缓冲区复制到设备。没有必要，因为 Mandelbrot 算法会写入整个设备缓冲区。
• 修复了不正确的网格大小计算。
• 删除了主机内存的 malloc，因为结果直接复制到 OpenCV 图像缓冲区中。
• 将缓冲区更改为使用字节而不是整数，当您有一个具有 8 位分辨率的灰色通道时，这会更方便。
• 删除了一些不必要的浮动转换。当您在计算中使用整数和浮点数时，整数会自动提升为浮点数。
• 修复了 Mandelbrot 算法中的两个问题：
  • x 和 y 被声明为 ints，而它们应该是 floats。
  • while 循环中的第一个表达式应包含 +，而不是 -。