The general consensus is that you can ease yourself into GPU programming using three approaches:
<ol>
	<li>Use libraries (like ArrayFire) and get quick turnaround times in your acceleration projects.</li>
	<li>Use compiler directives (like OpenACC) that highlight what part of the code can be run in parallel and let the compiler figure out how to best run the code on the GPU and finally,</li>
	<li>The do-it-yourself approach by redesigning your program to run on GPUs using CUDA or OpenCL.</li>
</ol>
In this blog, I will illustrate approach#3 using a C++ application for exploring the Mandelbrot set as a test case. This is a canonical problem that many people do to show enormous speedups.

The original code is given below. We'll call this the CPU version because it is designed to run on the CPU.
<pre class="lang:default decode:true" title="Mandelbrot.cpp">/*

Calculating which points in the specified window of the complex plane
lie inside the Mandelbrot set

CPU version

*/

#include &lt;iostream&gt;
#include &lt;complex&gt;
#include &lt;cmath&gt;
#include &lt;cstdlib&gt;

// -----------------------------
// declare some global variables

const int MAX_ITER = 200;    // maximum number of iterations

const int NX = 1000;  // number of points along X inside the window
const int NY = 1000;  // number of points along Y inside the window

// -------------------------------------------------------
// define a class for complex numbers and their operations

class dcmplx
{
public:
    double re;   // real component
    double im;   // imaginary component

// function to calculate the magnitude or absolute value of the complex number

double magnitude()
{
    return pow((re*re + im*im),0.5);
}

};

// ----------------------------------------------------------------------------------
// function to check all points inside the specified window for membership in the set

void Mandelbrot(double xmin, double xmax, double ymin, double ymax, int *color)
{
    double dx = (xmax - xmin)/NX; // grid spacing along X
    double dy = (ymax - ymin)/NY; // grid spacing along Y

    for(int i=0;i&lt;NX;i++)
    {
        for(int j=0;j&lt;NY;j++)
        {
            double x = xmin + (double) i*dx;   // actual x coordinate (real component)
            double y = ymin + (double) j*dy;   // actual y coordinate (imaginary component)

            dcmplx c;
            c.re = x;
            c.im = y;

            dcmplx z;
            z.re = 0.0;
            z.im = 0.0;

            // ---------------
            // z &lt;---- z*z + c

            int iter = 0;

            while(iter&lt;MAX_ITER)
            {
                iter++;
                dcmplx temp = z;
                z.re = temp.re*temp.re - temp.im*temp.im  +  c.re;
                z.im = 2.0*temp.re*temp.im                +  c.im;

                if (z.magnitude() &gt; 2.0) break;
            }

            // the 2D array "dev_color" stores how many iterations were required for divergence
            // for points outside the Mandelbrot set, this is typically a small number
            // points inside the set do not diverge and thus iter is a large number for such points

            color[i+j*NX] = iter;

            // GRAPHICS: the value of iter can be used to pick a suitable "color" when drawing the set
        }
    }
}

// ------------
// main program

int main(int argc, char* argv[])
{
    // allocate memory for a 2D array
    int *color = new int[NX*NY];

    // fill the 2D array with zeroes
    for(int i=0;i&lt;NX;i++)
    {
        for(int j=0;j&lt;NY;j++)
        {
            color[i+j*NX] = 0;
        }
    }

    // initialize parameters on the host
    double xmin = -2, xmax = 1, ymin = -1.5, ymax = 1.5;

    // check each pixel for membership and fill the 2D array accordingly
    Mandelbrot(xmin,xmax,ymin,ymax,color);

    return 0;
}</pre>
<strong>The CUDA approach</strong>

It is obvious that the CPU version is going to check each point in the 2D window for membership one after the other, in a serial fashion. However, all these calculations can be carried out in parallel because they are completely independent (embarrassingly parallel problem). With a modern GPU, one can launch thousands of parallel calculations simultaneously and, in theory, vastly reduce execution time compared to a single CPU.

The modified code using CUDA is provided below. Note how the loops inside the Mandelbrot function have been eliminated and this function (kernel) is launched on a GPU grid using the triple chevron notation in the main program. We also need to allocate space on the GPU for storing the 2D color array and copy this to the CPU once the GPU finishes filling it in.
<pre class="lang:c++ decode:true" title="GPU_Mandelbrot.cpp">/*

Calculating which points in the specified window of the complex plane
lie inside the Mandelbrot set

GPU version

*/

#include &lt;iostream&gt;
#include &lt;complex&gt;
#include &lt;cmath&gt;
#include &lt;cstdlib&gt;

// -----------------------------
// declare some global variables

const int MAX_ITER = 200;    // maximum number of iterations

const int NX = 1000;  // number of points along X inside the window
const int NY = 1000;  // number of points along Y inside the window

// -------------------------------------------------------
// define a class for complex numbers and their operations

class dcmplx
{
public:
    double re;   // real component
    double im;   // imaginary component

// function to calculate the magnitude or absolute value of the complex number
// this function is called from and executes on the device (GPU) 

__device__
double magnitude()
{
    return pow((re*re + im*im),0.5);
}

};

// ----------------------------------------------------------------------------------
// function to check all points inside the specified window for membership in the set
// this function is called from the host (CPU) but executes on the device (GPU)

__global__
void Mandelbrot(double xmin, double xmax, double ymin, double ymax, <strong>int *dev_color</strong>)
{
    double dx = (xmax - xmin)/NX; // grid spacing along X
    double dy = (ymax - ymin)/NY; // grid spacing along Y

    int i = blockIdx.x;
    int j = blockIdx.y;

    double x = xmin + (double) i*dx;   // actual x coordinate (real component)
    double y = ymin + (double) j*dy;   // actual y coordinate (imaginary component)

    dcmplx c;
    c.re = x;
    c.im = y;

    dcmplx z;
    z.re = 0.0;
    z.im = 0.0;

    // ---------------
    // z &lt;---- z*z + c

    int iter = 0;

    while(iter&lt;MAX_ITER)
    {
        iter++;
        dcmplx temp = z;
        z.re = temp.re*temp.re - temp.im*temp.im  +  c.re;
        z.im = 2.0*temp.re*temp.im                +  c.im;

        if (z.magnitude() &gt; 2.0) break;
    }

    // the 2D array "dev_color" stores how many iterations were required for divergence
    // for points outside the Mandelbrot set, this is typically a small number
    // points inside the set do not diverge and thus iter is a large number for such points

    dev_color[i+j*NX] = iter;

    // GRAPHICS: the value of iter can be used to pick a suitable "color" when drawing the set
}

// ---------------------------------------
// main program executes on the host (CPU)

int main(int argc, char* argv[])
{
    // allocate memory for a 2D array on the host (CPU)
    int *color = new int[NX*NY];

    // fill the 2D array with zeroes
    for(int i=0;i&lt;NX;i++)
    {
        for(int j=0;j&lt;NY;j++)
        {
            color[i+j*NX] = 0;
        }
    }

    // allocate memory on the device (GPU)
    int *dev_color;
    cudaMalloc((void **) &amp;dev_color, NX*NY*sizeof(int));

    // initialize parameters on the host
    double xmin = -2, xmax = 1, ymin = -1.5, ymax = 1.5;

    // check each pixel for membership and fill the 2D array accordingly
    // this kernel is executed in parallel using NX*NY blocks on the GPU

    dim3 grid(NX,NY);

    Mandelbrot&lt;&lt;&lt;grid,1&gt;&gt;&gt;(xmin,xmax,ymin,ymax,dev_color);

    // copy the 2D array from device (GPU) to host (CPU)
    cudaMemcpy(color, dev_color, NX*NY*sizeof(int), cudaMemcpyDeviceToHost);

    return 0;
}</pre>
<strong>SPEED UP</strong>

This is a tricky beast to define, but let us look at some of the results from running the CPU code and the GPU code. If we only focus on the kernel we ran on the GPU and compare it to the same function running on the CPU, <strong>the GPU version is ~400X faster!</strong>

<a href="http://blog.accelereyes.com/blog/wp-content/uploads/2013/03/Slide1.png"><img class="alignnone size-medium wp-image-3834" alt="Slide1" src="http://blog.accelereyes.com/blog/wp-content/uploads/2013/03/Slide1-300x225.png" width="300" height="225" /></a>

But, as seen in the picture above, the GPU version also has a significant amount of overhead in allocating memory on the device using "cudaMalloc" and in transferring data back to the CPU using "cudaMemcpy". You should keep this in mind while considering how to redesign your C/C++ code. The rule of thumb is to minimize data transfers to/from the CPU to realize the full potential of GPU acceleration.

<strong>CONCLUSIONS</strong>

You can see that porting an existing C++ code to CUDA is not so difficult. But it does involve rewriting portions of the code and moving data to and from the CPU memory to the GPU memory. In addition, there are several options for deploying the kernel (in this case, the Mandelbrot function that runs on the GPU is called a kernel) and making optimal use of the hardware is critical for achieving the maximum performance possible.

A lot of this guesswork can be eliminated by using ArrayFire. This library can be used in conjunction with your existing C/C++/Fortran 90/Python code and the learning curve is much easier compared to writing your own CUDA/OpenCL code. In Part 2 of this blog, I will outline how to use ArrayFire in the original C++ code and make full use of GPU acceleration without writing a single line of CUDA code.