Different result of computing in GPU (CUDA) and CPU

Question

I've wanted to create program that generates fractals on my GPU. First I created a working project in C, after that I tried to convert it into CUDA/C.

Unfortunately, after I did it I saw that there is a difference in results of CPU and GPU.

I spend few hours thinking what I did wrong and it's a mystery to me.

IMO: It seems that there is a difference of calculating values in while loop, therefore it ends earlier than in normal CPU function.

Question: is there any possibility that it is true? And if, what can I do to avoid that kind of computing error?

Here's my entire code:

// C libs
#include <stdint.h>
#include <stdio.h>
#include <iostream>

// Help libs
#include <windows.h>
#include <math.h>

// CUDA libs
#include "cuda.h"
#include "cuda_runtime.h"
#include "device_launch_parameters.h"

__global__ void calulateFractal(unsigned char *a, int N, double c_re, double c_im, int width, int height, double minX, double maxX, double minY, double maxY, double ratioX, double ratioY, int maxLevel)
{
    int i = blockIdx.x * blockDim.x + threadIdx.x;

    if(i < N)
    {
        int x = i % width;
        int y = i / width;

        double p_im = y * ratioY + minY;
        double p_re = x * ratioX + minX;

        double z_re = p_re;
        double z_im = p_im;

        int iteration = 0;

        while ((z_re * z_re + z_im * z_im) < 4 && iteration < maxLevel)
        {
            double tmp_re = z_re * z_re - z_im * z_im + c_re;
            double tmp_im = 2 * z_re * z_im + c_im;
            z_re = tmp_re;
            z_im = tmp_im;
            iteration++;
        }

        a[i] = iteration;
    }
}

void calulateFractalCPU(unsigned char *a, int i, double c_re, double c_im, int width, int height, double minX, double maxX, double minY, double maxY, double ratioX, double ratioY, int maxLevel)
{
    int x = i % width;
    int y = i / width;

    double p_im = y * ratioY + minY;
    double p_re = x * ratioX + minX;

    double z_re = p_re;
    double z_im = p_im;

    int iteration = 0;

    while ((z_re * z_re + z_im * z_im) < 4 && iteration < 99)
    {
        double tmp_re = z_re * z_re - z_im * z_im + c_re;
        double tmp_im = 2 * z_re * z_im + c_im;
        z_re = tmp_re;
        z_im = tmp_im;
        iteration++;
    }

    a[i] = iteration;
}

int saveFractalToBitmap(unsigned char **colorsArray, unsigned char *bitmap, int width, int height, char *filename)
{
    // Bitmap structures to be written to file
    BITMAPFILEHEADER bfh;
    BITMAPINFOHEADER bih;

    // Fill BITMAPFILEHEADER structure
    memcpy((char *)&bfh.bfType, "BM", 2);
    bfh.bfSize = sizeof(bfh) + sizeof(bih) + 3*height*width;
    bfh.bfReserved1 = 0;
    bfh.bfReserved2 = 0;
    bfh.bfOffBits = sizeof(bfh) + sizeof(bih);

    // Fill BITMAPINFOHEADER structure
    bih.biSize = sizeof(bih);
    bih.biWidth = width;
    bih.biHeight = height;
    bih.biPlanes = 1;
    bih.biBitCount = 24;
    bih.biCompression = BI_RGB; // uncompressed 24-bit RGB
    bih.biSizeImage = 0; // can be zero for BI_RGB bitmaps
    bih.biXPelsPerMeter = 3780; // 96dpi equivalent
    bih.biYPelsPerMeter = 3780;
    bih.biClrUsed = 0;
    bih.biClrImportant = 0;

    // Open bitmap file (binary mode)
    FILE *f;
    f = fopen(filename, "wb");

    if(f == NULL) 
        return -1;

    // Write bitmap file header
    fwrite(&bfh, 1, sizeof(bfh), f);
    fwrite(&bih, 1, sizeof(bih), f);

    // Write bitmap pixel data starting with the
    // bottom line of pixels, left hand side
    for (int i = 0; i < width * height ; i++)
    {
        // Write pixel components in BGR order
        fputc(colorsArray[bitmap[i]][2], f);
        fputc(colorsArray[bitmap[i]][1], f);
        fputc(colorsArray[bitmap[i]][0], f);
    }

    // Close bitmap file
    fclose(f);

    return 0;
}

int main()
{
    unsigned char **colorsArray;
    unsigned char *fractalLevelsCPU;
    unsigned char *fractalLevelsGPU;

    double minX = -1.7;
    double maxX = 1.7;
    double minY = -1.5;
    double maxY = 1.5;

    double input_re = -0.79;
    double input_im = 0.1463;

    int width = 10;
    int height = 5;
    int N = width * height;
    int maxLevel = 100;
    size_t levelsArraySize = N * sizeof(unsigned char);

    double ratioX = (maxX - minX) / (double) width;
    double ratioY = (maxY - minY) / (double) height;

    bool gpu = true;

    // Allocate memory
    colorsArray = (unsigned char**) malloc((maxLevel+1) * sizeof(unsigned char*));
    for(int i=0; i<=maxLevel; i++)
    {
        colorsArray[i] = (unsigned char *) malloc(3 * sizeof(unsigned char));
        colorsArray[i][0] = (int) (255.0 * i / maxLevel);
        colorsArray[i][1] = (int) (255.0 * i / maxLevel);
        colorsArray[i][2] = (int) (255.0 * log((double) i) / log((double) maxLevel));
    }

    fractalLevelsCPU = (unsigned char*) malloc(levelsArraySize);
    cudaMalloc((unsigned char **) &fractalLevelsGPU, levelsArraySize);

    cudaMemcpy(fractalLevelsCPU, fractalLevelsGPU, levelsArraySize, cudaMemcpyHostToDevice);

    if(gpu)
    {
        // Run GPU method
        calulateFractal <<< 1, N >>> (fractalLevelsGPU, N, input_re, input_im, width, height, minX, maxX, minY, maxY, ratioX, ratioY, maxLevel);

        // Copy data from GPU to CPU array
        cudaMemcpy(fractalLevelsCPU, fractalLevelsGPU, levelsArraySize, cudaMemcpyDeviceToHost);
    }
    else
    {
        // Iterate every element in array and compute level of fractal
        for(int i=0; i<N; i++)
        {
            calulateFractalCPU(fractalLevelsCPU, i, input_re, input_im, width, height, minX, maxX, minY, maxY, ratioX, ratioY, maxLevel);
        }
    }

    // Show results
    for(int i=0; i<N; i++)
    {
        if((i % width) == 0) 
            printf("\n");

        printf("%d\t", fractalLevelsCPU[i]);    
    }
    //saveFractalToBitmap(colorsArray, fractalLevelsCPU, width, height, "frac.bmp");

    // Free memory
    for(int i=0; i<=maxLevel; i++)
    {
        free(colorsArray[i]);
    }
    free(colorsArray);

    free(fractalLevelsCPU);
    cudaFree(fractalLevelsGPU);

    return 0;
}

Answer 1

I've find solution to my problem.

First of all, number of threads per block should be a power of two number. Also I realized that my GPU has it's limits for number of threads per block and blocks itself. NVIDIA Utils showed me that I can use max 65536 blocks and 512 threads per block.

Solution:

int threadsPerBlock = 512;
int blocksNumber = N/threadsPerBlock + (N % threadsPerBlock == 0 ? 0:1);

if(blocksNumber > 65536) 
    return -1;

calulateFractal <<< blocksNumber, threadsPerBlock >>> (fractalLevelsGPU, N, input_re, input_im, width, height, minX, maxX, minY, maxY, ratioX, ratioY, maxLevel);