0

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;
}
localhost
  • 93
  • 8
  • You seem to have forgotten to ask a question. – talonmies Mar 06 '13 at 14:14
  • You're right ;) Sorry, first time here. – localhost Mar 06 '13 at 14:18
  • I have no experience with CUDA, but this may be the source of the difference: http://www.herikstad.net/2009/05/cuda-and-double-precision-floating.html : *all doubles are silently converted into floats inside kernels and any double precision calculations computed are incorrect.* – Martin Tournoij Mar 06 '13 at 14:31
  • 1
    Are you sure you understand how cuda works? You do realize you are running calulateFractal N * N times on the GPU, at the same time (more or less). Also note that only the first N values are executed (where blockIDx is 0, it goes to N).Have a look here: [1]: https://devtalk.nvidia.com/default/topic/417278/cuda-programming-and-performance/difference-between-threadidx-blockidx-statements/ [2]: http://www.nvidia.com/content/cudazone/download/Getting_Started_w_CUDA_Training_NVISION08.pdf – Marius Brendmoe Mar 06 '13 at 14:34
  • 1
    @localhost... seems like a precision problem. – sgarizvi Mar 06 '13 at 14:35
  • Modern cuda (2.0+) supports double precision floating point see: http://stackoverflow.com/questions/2816992/double-precision-floating-point-in-cuda – Marius Brendmoe Mar 06 '13 at 14:38
  • There is no precision problem if you have CC of the GPU greater or equal to 2.0. In CC less than 2.0 double are not supported. Check the CC of your card and see if supports double precision. – KiaMorot Mar 06 '13 at 14:38
  • My GPU supports only 1.1 version. @Marius Brendmoe: I've tried everything and forgot to change it back. – localhost Mar 06 '13 at 14:43
  • 1
    You won't be able to use variables of type double, then. – KiaMorot Mar 06 '13 at 14:45
  • I changed all types from double to float and CPU computed all correctly. – localhost Mar 06 '13 at 14:46
  • 2
    @KiaMorot... Correct your info. [double precision floating point support started from Compute 1.3](http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html#features-technical-specifications.xml). – sgarizvi Mar 06 '13 at 15:29
  • yes, the answer with high probability is that it is a precision problem, also I don't think this program is correctly converted to CUDA. – Soroosh Bateni Mar 06 '13 at 17:37

1 Answers1

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);
localhost
  • 93
  • 8