Laplace's equation is one of the simplest possible partial differential equations to solve numerically. It applies to various physical processes like electric potentials and temperature fields. In this blog, I will highlight a simple numerical scheme to solve this equation referred to as Jacobi iteration. Starting with a simple CPU based C++ code, I will highlight how easy it is to port this code to the GPU using CUDA.

We'll focus on a particular example from heat transfer, where the goal is to solve for the steady-state temperature field T(x,y) and the solution domain is a square of size L in the X-Y plane. The boundary conditions are as follows:
<ul>
	<li>Left edge: x = 0, T(0, y) = 1</li>
	<li>Right edge: x = L, T(L, y) = 0</li>
	<li>Bottom edge: y = 0, T(x,0) = 0</li>
	<li>Top edge: y = L, T(x,L) = 0</li>
</ul>
Without going into too many details, the core algorithm for solving the temperature field inside the domain boils down to this simple equation:

T(i,j) = 0.25 * [T(i-1,j) + T(i+1,j) + T(i,j-1) + T(i,j+1)]       -------   (1)

Here, it is assumed that the square domain is represented using a computational grid of N x N points, equally spaced along X and Y. Equation (1) essentially states that the temperature value at any node point equals the arithmetic average of the 4 neighboring temperature values. The value at the boundary node points are known from the boundary conditions and the objective is to iteratively find the interior T(i,j) values using Eq.(1).

<strong>The CPU approach</strong>

We define two 2D arrays, T(1:N,1:N) and T_new(1:N,1:N) and essentially loop over all interior node points in the domain, updating T_new(i,j) at each interior node point using values from T(i,j) from the previous iteration. The core Jacobi function from the CPU code is copied below for reference.
<pre class="lang:default decode:true">void Jacobi(float* T, float* T_new)
{
    const int MAX_ITER = 1000;
    int iter = 0;

    while(iter &lt; MAX_ITER)
    {
        for(int i=1; i&lt;NX-1; i++)
        {
            for(int j=1; j&lt;NY-1; j++)
            {
                float T_E = T[(i+1) + NX*j];
                float T_W = T[(i-1) + NX*j];
                float T_N = T[i + NX*(j+1)];
                float T_S = T[i + NX*(j-1)];
                T_new[i+NX*j] = 0.25*(T_E + T_W + T_N + T_S);
            }
        }

        for(int i=1; i&lt;NX-1; i++)
        {
            for(int j=1; j&lt;NY-1; j++)
            {
                float T_E = T_new[(i+1) + NX*j];
                float T_W = T_new[(i-1) + NX*j];
                float T_N = T_new[i + NX*(j+1)];
                float T_S = T_new[i + NX*(j-1)];
                T[i+NX*j] = 0.25*(T_E + T_W + T_N + T_S);
            }
        }
        iter+=2;
    }
}</pre>
In this implementation, we update T_new(i,j) using values from T(i,j) and then update T(i,j) using values from T_new(i,j). This essentially completes two Jacobi iterations. We keep going until a set number of maximum iterations is reached (1,000 in the above code snippet). The CPU implementation is slow mainly because we update points in a sequential fashion, one after the other. But this was pretty much the best you could hope for before parallel computing.

<strong>The GPU approach: CUDA</strong>

Notice that each update of T_new(i,j) in the above code is independent of every other update. So what if we could update all node points in parallel at the same time? This is exactly what we do using GPUs. The entire CUDA code for solving this problem is copied below for reference.
<pre class="lang:default decode:true">/*

Solution of the Laplace equation for heat conduction in a square plate

*/

#include &lt;iostream&gt;

// global variables

const int NX = 256;      // mesh size (number of node points along X)
const int NY = 256;      // mesh size (number of node points along Y)

const int MAX_ITER=1000;  // number of Jacobi iterations

// device function to update the array T_new based on the values in array T_old
// note that all locations are updated simultaneously on the GPU 
__global__ void Laplace(float *T_old, float *T_new)
{
    // compute the "i" and "j" location of the node point
    // handled by this thread

    int i = blockIdx.x * blockDim.x + threadIdx.x ;
    int j = blockIdx.y * blockDim.y + threadIdx.y ;

    // get the natural index values of node (i,j) and its neighboring nodes
                                //                         N 
    int P = i + j*NX;           // node (i,j)              |
    int N = i + (j+1)*NX;       // node (i,j+1)            |
    int S = i + (j-1)*NX;       // node (i,j-1)     W ---- P ---- E
    int E = (i+1) + j*NX;       // node (i+1,j)            |
    int W = (i-1) + j*NX;       // node (i-1,j)            |
                                //                         S 

    // only update "interior" node points
    if(i&gt;0 &amp;&amp; i&lt;NX-1 &amp;&amp; j&gt;0 &amp;&amp; j&lt;NY-1) {
        T_new[P] = 0.25*( T_old[E] + T_old[W] + T_old[N] + T_old[S] );
    }
}

// initialization

void Initialize(float *TEMPERATURE)
{
    for(int i=0;i&lt;NX;i++) {
        for(int j=0;j&lt;NY;j++) {
            int index = i + j*NX;
            TEMPERATURE[index]=0.0;
        }
    }

    // set left wall to 1

    for(int j=0;j&lt;NY;j++) {
        int index = j*NX;
        TEMPERATURE[index]=1.0;
    }
}

int main(int argc,char **argv)
{
    float *T;          // pointer to host (CPU) memory
    float *_T1, *_T2;  // pointers to device (GPU) memory

    // allocate a "pre-computation" T array on the host
    float *T = new float [NX*NY];

    // initialize array on the host
    Initialize(T);

    // allocate storage space on the GPU
    cudaMalloc((void **)&amp;_T1,NX*NY*sizeof(float));
    cudaMalloc((void **)&amp;_T2,NX*NY*sizeof(float));

    // copy (initialized) host arrays to the GPU memory from CPU memory
    cudaMemcpy(_T1,T,NX*NY*sizeof(float),cudaMemcpyHostToDevice);
    cudaMemcpy(_T2,T,NX*NY*sizeof(float),cudaMemcpyHostToDevice);

    // assign a 2D distribution of CUDA "threads" within each CUDA "block"    
    int ThreadsPerBlock=16;
    dim3 dimBlock( ThreadsPerBlock, ThreadsPerBlock );

    // calculate number of blocks along X and Y in a 2D CUDA "grid"
    dim3 dimGrid( ceil(float(NX)/float(dimBlock.x)), ceil(float(NY)/float(dimBlock.y)), 1 );

    // begin Jacobi iteration
    int k = 0;
    while(k&lt;MAX_ITER) {
        Laplace&lt;&lt;&lt;dimGrid, dimBlock&gt;&gt;&gt;(_T1,_T2);   // update T1 using data stored in T2
        Laplace&lt;&lt;&lt;dimGrid, dimBlock&gt;&gt;&gt;(_T2,_T1);   // update T2 using data stored in T1
        k+=2;
    }

    // copy final array to the CPU from the GPU 
    cudaMemcpy(T,_T2,NX*NY*sizeof(float),cudaMemcpyDeviceToHost);
    cudaThreadSynchronize();
/*
    // print the results to screen
    for (int j=NY-1;j&gt;=0;j--) {
        for (int i=0;i&lt;NX;i++) {
            int index = i + j*NX;
            std::cout &lt;&lt; T[index] &lt;&lt; " ";
        }
        std::cout &lt;&lt; std::endl;
    }
*/
    // release memory on the host 
    delete T;

    // release memory on the device 
    cudaFree(_T1);
    cudaFree(_T2);

    return 0;
}</pre>
Notice that the basic instruction for updating the node point using the average of its neighbors is identical in the GPU code. However, the loop over i and j has disappeared in the kernel (function running on the GPU). The key point to understand is that this kernel is executed in parallel using several threads on the GPU using the triple chevron notation &lt;&lt;&lt;   &gt;&gt;&gt;.

Each thread figures out which node point it is updating and updates the value in the appropriate memory location. When all threads are done updating their node point, we move on to the next iteration. Obviously, the GPU approach is faster than the CPU approach, but there is also some time spent in allocating resources on the GPU and in copying data from the CPU to the GPU. Thus, the speed up in the actual execution of the Jacobi iteration may not become apparent until the mesh size is large.

<strong>The GPU approach: ArrayFire</strong>

Now that you have a feel of how to write a simple CUDA kernel for implementing Jacobi iteration, the final step in our journey will be to implement the GPU version of this code using the ArrayFire library. At the core of the code to follow is a ArrayFire function called convolution. This essentially alters an array (1D, 2D or 3D) using a pre-defined "convolution kernel". In this particular example, we have a 3 x 3 kernel which tells the code how to calculate the new value at a certain location using neighboring values. If we are only concerned with averaging neighbors along the vertical and horizontal directions, this kernel can be visualized as follows:
<pre class="lang:default decode:true">0      1/4      0

1/4     0      1/4

0      1/4      0</pre>
During the convolution operation, this kernel is used to update the value at the central location, using the weights for the neighboring points as provided above. This is just another way of implementing Eq.(1) above.

The ArrayFire code itself is much simpler than the CUDA version copied above. The complete code for implementing Jacobi iteration is provided below.
<pre class="lang:default decode:true">#include &lt;stdio.h&gt;
#include &lt;arrayfire.h&gt;
using namespace af;

const int N = 4096;

int main(int argc, char** argv)
{
    array T = constant(0,N,N);
    T(seq(1,end-1),0) = 1;

//  print(T);

    float fval[]={0,0.25,0,0.25,0,0.25,0,0.25,0};
    array f(3,3,fval);

    int iter=0;
    while(iter&lt;1000) {
        T= convolve(T, f);
        T(seq(1,end-1),0) = 1;
        T(seq(1,end-1),end) = 0;
        T(0,seq(1,end-1)) = 0;
        T(end,seq(1,end-1)) = 0;
        iter+=1;
    }

//  print(T);

    return 0;
}</pre>
&nbsp;

<strong>RESULTS</strong>

In the figure below, a comparison of the run-time between the CPU and the GPU is attached. For the largest grid shown below, we obtained a speed up of 143X for the CUDA version compared to the CPU code.
<p style="text-align: center;"><a href="http://blog.accelereyes.com/blog/wp-content/uploads/2013/03/Screen-Shot-2013-03-26-at-5.46.18-PM.png"><img class="aligncenter  wp-image-3942" alt="Screen Shot 2013-03-26 at 5.46.18 PM" src="http://blog.accelereyes.com/blog/wp-content/uploads/2013/03/Screen-Shot-2013-03-26-at-5.46.18-PM-1024x781.png" width="614" height="469" /></a></p>
It will be apparent that the GPU version using CUDA is way faster than the CPU version for mesh sizes larger than 256 x 256. In reality, as mentioned above, the actual Jacobi iteration (or convolution) part is <em>always</em> faster than the CPU version (for both raw CUDA and ArrayFire). It will also be seen that the raw CUDA version is a little faster (~2X) than the ArrayFire version. This is mainly because the convolution operation in ArrayFire accounts for all neighbors of the point and the point itself (9 values) while the raw CUDA version skips these additional calculations (which sum to zero).

But the simplicity of the ArrayFire version is a big win for those new to GPU computing and those who wish to get the benefits of GPU computing without writing their own CUDA kernels.

To conclude, this example illustrated how a simple partial differential equation (PDE) can be solved on the GPU provided we can execute it as a "single instruction on multiple data" (SIMD) algorithm. The basic idea of creating a workspace (array) on the GPU and updating elements of this array in parallel can be applied even to more complex PDEs and algorithms, provided they can be written in a SIMD fashion. The only limitation in this approach is the amount of memory available on the GPU.

&nbsp;