A kernel több észre tagolódik. A központi részben egy a hullámok rajzolásához hasonló kórdész áll.

__global__ void juliaKernel(unsigned char *ptr, int w, int h, float scale) {

    int i_x = blockIdx.x*blockDim.x + threadIdx.x;
    int i_y = blockIdx.y*blockDim.y + threadIdx.y;

    float x = scale * (w / 2.0f - i_x) / (w / 2.0f);
    float y = scale * (h / 2.0f - i_y) / (h / 2.0f);

    int tid = (i_y*w + i_x) * 4;

    if (x < w && y < h) {
        int value = onePixel(x, y);
        ptr[tid + 0] = 255 * value;
        ptr[tid + 1] = 255 * value;
        ptr[tid + 2] = 0;
        ptr[tid + 3] = 255;
    }

    return;
}

A különbség főleg abban áll, hogy a pixel koordinátáját át transzformáljuk, hogy a kép közepéhez viszonytva legyen megadva.Ez után meghívjuk a onePixel() függvényt, ami egy logikai értéket ad vissza, ami alapján eldöntjuk a pixel színét.

Korábban írtuk, hogy GPU kódnak két fajtája van. A __device__ előtagú kódok például GPU-ról GPU-ra hívható kódok.

A onepixel() függvény a julia halmaz közeléítését adja. Ennek lényege, hogy egy iteratív számítás eredménye mondja meg, hogy a pixelt kiszínezzük-e.

• A kapott x, y koordinátáknból egy c = x + y*i formájó komplex számot készítünk.
• Vesszük a kapott c számot és egy konstans "k" komplex számot, és a végtelenségig iteráljuk a "c = c*c + k" számítást.
  • Természetesen végtelenségig nem tudunk számolni, de 200 iteráció is megfelelő.
• Ha az iteráció konvergens, akkor igaz értéket adunk vissza, ha nem akkor hamisat.
  • A sorozat 0-hoz konvergál, ezért ennek eldöntésére megvizsgálhatjuk, hogy az eredményt bagyobb-e mint eg elég nagy konstans.

Ez főleg matematikai számítás volt. Ami a programozástechnikai érdekességet adja, az az, hogy a komplex számokat külön megvalósíthatjuk egy struktúrában.

struct gpuComplex {
    float   r;
    float   i;
    __device__ gpuComplex(float _r, float _i) : r(_r), i(_i) {}

    __device__ float magnitudeSqr(void) {
        return r * r + i * i;
    }
    __device__ gpuComplex operator*(const gpuComplex& other) {
        return gpuComplex(r*other.r - i*other.i, i*other.r + r*other.i);
    }
    __device__ gpuComplex operator+(const gpuComplex& other) {
        return gpuComplex(r + other.r, i + other.i);
    }
};

A struktúrhoz viszont tudunk operációkat is adni és a C++ operátor overloading mechanizmusával definiálhatjuk rá az összeadás, és szorzás műveleteket is. Ezzel nagyban megkönnyítjük a dolgunkat. A fenti példában ez látható.

Fontos viszont, hogy a megadott operátorok, a GPU-n futnak így meg kell adnunk hozzájuk a __device__ előtagot.

#include "cuda_runtime.h"
#include "device_launch_parameters.h"
#include "math_constants.h"

#include "ImageIO.h"

#include 
#include 

using namespace std;

#define N           512
#define BLOCKDIM    16
#define SCALE       1.5

#define step_r      -0.8
#define step_i       0.156

void cudaCheckError(char* errInfo);

__global__ void juliaKernel(unsigned char* result, int w, int h, float waveLength);
__device__ int onePixel(float x, float y);

struct gpuComplex {
    float   r;
    float   i;
    __device__ gpuComplex(float _r, float _i) : r(_r), i(_i) {}

    __device__ float magnitudeSqr(void) {
        return r * r + i * i;
    }
    __device__ gpuComplex operator*(const gpuComplex& other) {
        return gpuComplex(r*other.r - i*other.i, i*other.r + r*other.i);
    }
    __device__ gpuComplex operator+(const gpuComplex& other) {
        return gpuComplex(r + other.r, i + other.i);
    }
};

int main(int argc, char** argv)
{
    unsigned char* result = new unsigned char[N*N * 4];
    unsigned char* dev_result;

    cudaMalloc((void**)&dev_result, N*N * 4 * sizeof(unsigned char));
    cudaCheckError("Memory allocation.");

    dim3 blockDim = dim3(BLOCKDIM, BLOCKDIM, 1);
    dim3 gridDim = dim3((N + BLOCKDIM - 1) / BLOCKDIM, (N + BLOCKDIM - 1) / BLOCKDIM, 1);

    juliaKernel << <gridDim, blockDim >> > (dev_result, N, N, SCALE);
    cudaCheckError("Kernel call.");

    cudaMemcpy(result, dev_result, N*N * 4 * sizeof(unsigned char), cudaMemcpyDeviceToHost);
    cudaCheckError("Memcpy: dev -> host");

    writeRGBImageToFile("image.png", result, N, N);

    return 0;
}

__global__ void juliaKernel(unsigned char *ptr, int w, int h, float scale) {

    int i_x = blockIdx.x*blockDim.x + threadIdx.x;
    int i_y = blockIdx.y*blockDim.y + threadIdx.y;

    float x = scale * (w / 2.0f - i_x) / (w / 2.0f);
    float y = scale * (h / 2.0f - i_y) / (h / 2.0f);

    int tid = (i_y*w + i_x) * 4;

    if (x < w && y < h) {
        int value = onePixel(x, y);
        ptr[tid + 0] = 255 * value;
        ptr[tid + 1] = 255 * value;
        ptr[tid + 2] = 0;
        ptr[tid + 3] = 255;
    }

    return;
}

__device__ int onePixel(float x, float y) {

    gpuComplex c(step_r, step_i);
    gpuComplex a(x, y);

    int i = 0;
    for (i = 0; i<200; i++) {
        a = a * a + c;
        if (a.magnitudeSqr() > 1000)
            return 0;
    }

    return 1;
}

void cudaCheckError(char* errInfo)
{
    cudaError_t errCode = cudaGetLastError();

    if (errCode != 0)
    {
        cout << "CUDA error occured: " << endl;
        cout << " - Error description: " << cudaGetErrorString(errCode) << endl;
        cout << " - Error info: " << errInfo << endl;
    }
}