Heterogeneous Parallel Programming: Dive into the World of CUDA

A previous article in this series titled ‘Introducing NVIDIAs CUDA’ covered the basics of the NVIDIA CUDA device architecture. This article covers parallel programming using CUDA C with sequential and parallel implementations of a vector addition program.

Parallel programming and general-purpose GPU computing are some of the hottest trends in computer science today due to the decreased prices of multi-core systems and the increase in compute efficiency. Various parallel programming languages like OpenCL and CUDA have been developed and evaluated over the years. This article will cover the basics of CUDA C, invoking kernels, threads and blocks with a vector addition program and it aims to give you an insight into beginning programming on your CUDA device.
A few important points before we begin. To use CUDA on your system, you will need to have the following installed:
1.     CUDA-capable GPU hardware
2.     A supported version of Linux with a GCC compiler and toolchain
3.     NVIDIA CUDA toolkit and drivers
It is presumed that if you already have a CUDA device within your system, you probably have the latest toolkit and drivers installed and configured correctly. In case you do not have the NVIDIA CUDA drivers configured or you have recently upgraded your hardware with a CUDA device, you could follow the simple steps given below to configure your device.
1.     Download the toolkit from the NVIDIA website, available at no cost from: http://www.nvidia.com/content/cuda/cuda-downloads.html. Select the correct product release depending on your operating system preferences. This download contains an all-in-one package, which includes the CUDA toolkit, SDK code samples and the required drivers. After downloading, follow the steps in the NVIDIA guide to install the drivers, CUDA samples and the toolkit, at:http://developer.download.nvidia.com/compute/cuda/5_0/rel/docs/CUDA_Getting_Started_Guide_For_Linux.pdf.This guide will help you set up the complete environment on the system.

In this article, I will cover the serial and parallel versions of a vector addition program. Once you have understood the basics of the parallel vector addition program, you could use the concepts pretty well in parallelising other algorithms as per your requirements.
First, let’s write a simple C program to perform vector addition of two arrays. Open your favourite editor and write a simple vector addition code that looks like what’s shown below:

#include<stdio.h>   static const int N=100;   //Add the vectors and store result in vector C void vector_add(int *a,int *b, int *c) {      int i=0;      for (i=1;i<=N;i++)     {          c[i]=a[i]+b[i];          }      }        int main()       {             int a[N], b[N], c[N];           int i=0;            //Initialize the vectors with values from 1 to 100 and its double in another array           {              a[i]=i;              b[i]=2*i;           }                //Call function vector_add to display the result              vector_add(a,b,c);                 //Print the resultant array.               for (i=1;i<=N;i++)               {                  printf("%d %d %d\n\n", a[i],b[i],c[i] );                  }                system("pause");               return 0;             }

This very simple serial vector addition program creates two arrays of integer values, and adds them using the vector_add function.
Compile the code using the following command:

gcc sequential_vector.c  –o sequential

Run it using the command given below:

./sequential

In the above program, the processor runs each task sequentially, one after the other. Looping in the above program works sequentially as it starts with the first index and computes consequentially till the last index, and then exits the program. It is a single-threaded execution. Now, CUDA gives us the functionality to perform the same operation in parallel. What it does essentially is offload the data parallel sections to the GPU device and send the result back after computation. In what follows, you will get an insight into launching kernels, writing a device and host code, and performing the same serial vector addition program given above, in parallel.
Here is a simple vector addition code in CUDA.

#include <cuda.h> #include<stdio.h> #define N 100 #define numThread 1 // in this example we keep one thread in one block #define numBlock 100 // in this example we use 100 blocks   __global__ void vector_add( int *a, int *b, int *c ) {      // keep track of the index    int tid = blockIdx.x;      while (tid < N) {         c[tid] = a[tid] + b[tid];         tid = tid+ numBlock; // shift by the total number of blocks, i.e. 100 in our case      } }   int main( void ) {     int *a, *b, *c;     int *dev_a, *dev_b, *dev_c;       // allocate the memory on the CPU    a = (int*)malloc( N * sizeof(int) );    b = (int*)malloc( N * sizeof(int) );    c = (int*)malloc( N * sizeof(int) );     //Initialize the vectors with values from 1 to 100 and its double in another array  for (int i=0; i<N; i++) {         a[i] = i;         b[i] = 2 * i;    }      // allocate the memory on the GPU    cudaMalloc( (void**)&dev_a, N * sizeof(int) ) ;   cudaMalloc( (void**)&dev_b, N * sizeof(int) ); cudaMalloc( (void**)&dev_c, N * sizeof(int) );        // copy the arrays 'a' and 'b' to the GPU     cudaMemcpy( dev_a, a, N * sizeof(int),cudaMemcpyHostToDevice );    cudaMemcpy( dev_b, b, N * sizeof(int),cudaMemcpyHostToDevice );      vector_add<<<numBlock,numThread>>>( dev_a, dev_b, dev_c );      // copy the array 'c' back from the GPU to the CPU    cudaMemcpy( c, dev_c, N * sizeof(int),cudaMemcpyDeviceToHost ) ;   //Prints the results for (int i=0; i<N; i++) { printf("%d %d %d \n\n", a[i],b[i],c[i] ); }        // free the memory we allocated on the CPU     free( a );     free( b );     free( c );      // free the memory we allocated on the GPU    cudaFree( dev_a ) ;    cudaFree( dev_b ) ;     cudaFree( dev_c);       return 0; }

Let’s begin with analysing each part of the code, and then compile the code to get our results.

From the previous article (Part 1 of this series), you know that CUDA programs execute in two places—the host (your CPU) and the device (GPU). You might be a bit surprised by the fact that writing the device code is much simpler than writing the CPU host code. Hence, let’s begin with analysing the device code first.

Since we have 100 array values in this code, to simplify things, let’s have 100 blocks (kernels) running simultaneously, where each kernel runs a single thread. Hence, let’s set numThread to 1 and numBlock to 100, and use these variables later while calling the device from the host.
The device code is:

__global__ void vector_add(int *a, int *b, int *c)      // keep track of the index     int tid = blockIdx.x;       while (tid < N) {         c[tid] = a[tid] + b[tid];         tid = tid+ numBlock; // shift by the total number of blocks, i.e. 100 in our case     } }

As shown above, add a __global__ qualifier to the function vector_add. Notice that there are very few changes in the function vector_add of the serial and parallel sections. The __global__ qualifier indicates that this is a device function that would be called from the host.

blockIdx is a built-in CUDA runtime variable, which is a three-component vector to identify threads in a one, two and three dimension index. Imagine a block as a 3-D matrix and to access the different components in this vector, use blockIdx.x, blockIdx.y and blockIdx.z. In this code, we will be using 100 blocks with a single thread on every grid, which will be seen while we analyse the host code, and hence we use only blockIdx.x which returns the current block number.

The condition while tid<N checks that the bounds for array computation have not been reached and computes the array sum taking the block value as an index, i.e., tid.
Add numBlock to the tid value, to shift the index by the number of blocks, as each block would be computing just one array index, and we have 100 blocks for 100 array indexes. This explanation pretty much sums up the device code for the program.
Now move on to the host code, which prepares the GPU for execution and invokes the kernel. It works by allocating memory to the GPU and CPU, transfers the input vectors to the GPU, launches the kernel and transfers the result back to the host (CPU).

int main( void ) {     int *a, *b, *c;     int *dev_a, *dev_b, *dev_c;       // allocate the memory on the CPU     a = (int*)malloc( N * sizeof(int) );     b = (int*)malloc( N * sizeof(int) );     c = (int*)malloc( N * sizeof(int) );       // fill the arrays 'a' and 'b' on the CPU    for (int i=0; i<N; i++) {         a[i] = i;        b[i] = 2 * i;     }       // allocate the memory on the GPU    cudaMalloc( (void**)&dev_a, N * sizeof(int) ) ;    cudaMalloc( (void**)&dev_b, N * sizeof(int) ); cudaMalloc( (void**)&dev_c, N * sizeof(int) );      // copy the arrays 'a' and 'b' to the GPU     cudaMemcpy( dev_a, a, N * sizeof(int),cudaMemcpyHostToDevice );   cudaMemcpy( dev_b, b, N * sizeof(int),cudaMemcpyHostToDevice );

As in the above code, similar to the allocation in C, variables a, b and c are allocated memory on the CPU. cudaMalloc() is a standard sub-routine of the CUDA API to allocate memory on the device. It works similar to the C malloc() function we used earlier.
Since you cannot modify the memory allocated to the device from the host directly in CUDA, cudaMemcpy() is used to transfer data to the device. This method takes a pointer to the local memory, a pointer to the GPU memory being copied to, the number of bytes that will be copied, and a flag which determines the direction of the memory transfer, respectively.

vector_add<<<numBlock,numThread>>>( dev_a, dev_b, dev_c );

This line is a call to the device kernel from the host to execute the function on the device. It is similar to the serial function call with some additional code. Blocks are organised in three dimensional grids and threads are organised in three dimensional blocks. numBlock and numThread are passed as arguments to let the device know about the structure to be adopted for computation.

cudaMemcpy( c, dev_c, N * sizeof(int),cudaMemcpyDeviceToHost ) ;

This copies the resultant data back from the GPU to the host. As you can see, it is similar to the cudaMemcpy() we used above, with the last variable being changed to cudaMemcpyDeviceToHost to indicate data transfer is between device to host.
The rest of the code is pretty much self-explanatory, except that you use cudaFree() to free the memory allocated to the GPU.
Now that you understand the code well, compile the code to verify your results.
Save the code and parallel_vector.cu, and type the following command in the terminal:

nvcc parallel_vector.cu –o parallel

Now, to execute the code, run the following command:

./parallel

If you have followed this guide correctly, it will print the vectors ‘a’ and ‘b’ along with their additional resultant ‘c’. This would verify that the first GPU code which you wrote has worked correctly. Still confused? Well, have a look at Figure 2, which will give you a clear understanding of how things work in parallel, in this case.

Performance
You may wonder how the GPU code can perform about 100 times faster than the CPU code, since we have created 100 blocks that are executing in parallel. This is not the case, since there is an overhead involved in copying data between the CPU and the GPU and the resultant data back to the CPU. Hence, CUDA is generally used for computing algorithms that are significantly data intensive, as it would then make sense to spend some time for data transfer. GPUs are, therefore, generally known as data intensive computational devices.

As a next step, you could try programming matrix addition on the GPU in parallel to get a good grip on kernels, threads and parallel execution. Your best companion for this would be the links and the books mentioned in the ‘References’ section at the end of this article. I also recommend you visit the NVIDIA website and documentation, as it will give you a good idea about the power of CUDA if you are not already impressed by what this simple GPU device on your laptop can do. Next up in this series, I might cover an advanced CUDA program with multiple threads and blocks on a grid, and analyse the running time of the code. I will follow it up with a discussion on OpenACC and other simpler parallel programming models that have come up recently.
Till then,  start thinking of algorithms in parallel. The world of parallel computing is here to stay!

References
[1]    ‘CUDA C Programming Guide’ by NVIDIA; http://docs.nvidia.com/cuda/cuda-c-programming-guide/index.html
[2]     ‘CUDA Application Design and Development’ by Rob Farber
[3]     ‘Programming Massively Parallel Processors’ by David B. Kirk and Wen-mei W. Hwu

Be Cautious while using Bit Fields for Programming

In this article, the author tells embedded C programmers why writing code involving bit fields needs to be done carefully, failing which the results may not be as expected. However, bit fields are handy as they involve only low level programming and result in efficient data storage.

Embedded C programmers who have worked with structures should also be familiar with the use of bit fields structures. Use of bit fields is one of the key optimisation methods in embedded C programming, because these allow one to pack together several related entities, where each set of bits and single bits can be addressed. Of course, the usage of bit fields is ‘easy’ and comes handy, especially during low level programming. Though considered as one of the unique features of C programming, bit fields do have some limitations. Let us look at these by exploring the example problems in this article.

Data types and bit fields
Let us look into the signed qualifiers affecting the output of the bit field structure. Please note that the code snippets provided here are tested with the GCC compiler [gcc version 4.7.3] running under a Linux environment.
Let us consider a simple small C code snippet as shown below, with a structure named bit field, with three integer fields—hours, mins and secs, of bit field sizes 5, 6 and 6, respectively:

typedef struct bit_field
{
int hours : 5;
int mins : 6;
int secs : 6;
}time_t;

Now let us declare a variable alarm of type time_t and set values as 22, 12 and 20, respectively:

//Declaration of the variable of type time_t
time_t alarm;
/Assigning the values to the different members of the bit-field structures
alarm.hours = 22;
alarm.mins  = 12;
alarm.secs  = 20;

When we print these values using a simple printf statement, what could be the output? At first, most of us will envision the answers to be 22, 12 and 20, for hours, mins and secs respectively. Whereas when we actually compile and run the run the code, the value printed for the hours would be different – 10, 12 and 20 (as shown in Figure 1).

Figure 1: Actual output

Where did we go wrong?
1. We all know that the default signed qualifier for the ‘int’ is ‘signed int’.
2.     We reserved 5 bits for storing the hours field assuming we were using the 24-hour format. From among 5 bits, 1 bit was used for storing the sign of the number, which means only 4 bits were then available for storing the actual value. In these 4 bits, we can store the numbers ranging from -16 to +15 according to the formula (-2^k) to ([+2^k] -1)) including ‘0’, where  ‘k’  indicates the number of bits.

3.    We will see how 22 is stored in binary form in 5 bits through pictorial representation (Figure 2).
4.    From the table(as shown in Figure 2), it is very clear that sign bit (b4) is SET, which indicates the value is negative. So, when printed using the printf statement, we will get -10 (the decimal value of 10110), because of which we got an unexpected output.

Figure 2 : Pictorial representation of binary value of 24 in ‘5’ bits

Now that we have understood the problem, how do we fix it? It is very simple; just qualify ‘int’ to ‘unsigned int’ just before the hours in the bit field structure, as shown below. The corrected output is shown in Figure 3.

Figure 3 : Correct output after correct usage of datatype

#include <stdio.h>
typedef struct bit_field
{
unsigned int hours  : 5;
unsigned int mins   : 6;
unsigned int secs    : 6;
}time_t;
int main()
{
//Declaration of the variable of type time_t
time_t alarm;
//Assigning the values to the different members of the bit-field structures
alarm.hours = 22;
alarm.mins  = 12;
alarm.secs  = 20;
printf(“Hours : %d\nMins : %d\nSecs : %d\n”, alarm.hours, alarm.mins, alarm.secs);
}

Bit wise operators definitely provide advantages, but they need to be used a ‘bit’ carefully. In the embedded programming environment, they might lead to major issues in case they are not handled properly.

Endianess of the architecture and bit fields
In this problem, we will see how Endianess affects the bit fields. Bit fields in C always start at Bit 0, which is the least significant bit (LSB) on Little Endian. But most compilers on Big Endian systems inconveniently consider the most significant bit (MSB)—Bit 0.
Note: Big Endian machines pack bit fields from the most significant byte to the least significant.

Little Endian machines pack bit fields from the least significant byte to the most.
To start with, let us consider the code( Labelled as byte_order.c) given below:

1 #include <stdio.h>
2 typedef union {
3         unsigned int value;
4         struct {
5                 unsigned char one   : 8;
6                 unsigned char two   : 8;
7                 unsigned char three : 8;
8                 unsigned char four  : 8;
9         } bit_field;
10 } data_t;
11
12 int main() {
13
14         data_t var = {0x1A1B1C1D};
15         unsigned char *ptr = (unsigned char *)(&var);
16
17         printf(“The entire hex value is 0x%X\n”, var.value);
18         printf(“The first byte is 0x%X @ %p\n”, *(ptr + 0), ptr + 0);
19         printf(“The second byte is 0x%X @ %p\n”, *(ptr + 1), ptr + 1);
20         printf(“The third byte is 0x%X @ %p\n”, *(ptr + 2), ptr + 2);
21         printf(“The fourth byte is 0x%X @ %p\n”, *(ptr + 3), ptr + 3);
22
23         return 0;
24 }

Figure 4 : Output of the code byte_order.c

When I run this code in my system, I get the output shown in Figure 4.

Figure 5 : Byte-ordering in Little – Endianess Machine

Figure 6 : Byte-ordering in Big – Endianess Machine

From Figure 4, it is very clear that the underlying architecture is following the little Endian. When the same code is run under a different architecture, which follows Big Endian, the result will be different. So, portability issues need to be considered while using bit fields.
Let’s look at one more example to understand how bits are packed in Big Endian and Little Endian.
To start with, let us consider the sample code(Labelled as bit_order.c) given below:

1 #include <stdio.h>
2 typedef union {
3         unsigned short value;
4         struct {
5                 unsigned short v1   : 1;
6                 unsigned short v2   : 2;
7                 unsigned short v3   : 3;
8                 unsigned short v4   : 4;
9                 unsigned short v5   : 5;
10         } bit;
11 } data_t;
12
13 int main() {
14
15         data_t var ;
16         unsigned char *ptr = (unsigned char*)(&var);
17         var.bit.v1   = 1;
18         var.bit.v2   = 2;
19         var.bit.v3   = 3;
20         var.bit.v4   = 4;
21         var.bit.v5   = 5;
22
23         printf(“The Entire hex value is 0x%X\n”, var.value);
24         printf(“The first byte is 0x%X @ %p\n”, *(ptr + 0), ptr + 0);
25         printf(“The second byte is 0x%X @ %p\n”, *(ptr + 1), ptr + 1);
26
27         return 0;
28 }

Figure 7 : output of the code bit_order.c

When I run this code in my system, I get the output as shown in Figure 7.

From this figure, one can see that the bits are packed from the least significant on a little Endian machine. Figure 8 helps us understand how the bits ordering takes place.

Figure 8 : Bit – ordering in small Endianess architecture

If you run the same code in big Endian architecture, you will get the output given in Figure 9.

Figure 9 : Expected output of the code bit_order.c when run in big-endian architecture.

For more clarity, see Figure 10.
From the last two examples, it is very clear that bit fields pose serious portability issues. When the same programs are compiled on different systems, they may not work properly. This is because some C compilers use the left-to-right order, while other C compilers use the right-to-left order. They also have architecture-specific bit orders and packing issues.

Figure 10 : Bit- ordering in Big Endianess Architecture

As a concluding note, let us list the advantages and limitations of bit fields structures.
Advantages
1.     Efficiency – Storage of data structures by packing.
2.     Readability – Members can be easily addressed by the names assigned to them.
3.     Low level programming – The biggest advantage of bit fields is that one does not have to keep track of how flags and masks actually map to the memory. Once the structure is defined, one is completely abstracted from the memory representation as in the case of bit-wise operations, during which one has to keep track of all the shifts and masks.

Limitations
1.     As we saw earlier, bit fields result in non-portable code. Also, the bit field length has a high dependency on word size.
2.    Reading (using scanf) and using pointers on bit fields is not possible due to non-addressability.
3.    Bit fields are used to pack more variables into a smaller data space, but cause the compiler to generate additional code to manipulate these variables. This results in an increase in both space as well as time complexities.
4. The sizeof() operator cannot be applied to the bit fields, since sizeof() yields the result in bytes and not in bits.