CUDA 5.0 - optimising vector addition

This time I plan to concentrate just on optimising my vector addition a little bit. In my previous example I had a statically defined number of blocks and threads used by my kernel 'AddVectorsKernel'. This time I'd like to make it a bit more adjustable. I also got a slightly better syntax highlighter because the old one was pretty much unreadable.

To begin with, this is my vanilla code without the CPU vector addition, just plain CUDA. I will be introducing and explaining some changes to it as I go along.

#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <time.h>

#pragma comment(lib, "cudart") 

typedef struct 
{
    float *content;
    const unsigned int size;
} pjVector_t;

__global__ void AddVectorsKernel(float *firstVector, float *secondVector, float *resultVector)
{
    unsigned int index = threadIdx.x + blockIdx.x * blockDim.x;
    resultVector[index] = firstVector[index] + secondVector[index];
}

int main(void)
{
    const unsigned int vectorLength = 1000000;
    const unsigned int blocks = 1000;
    const unsigned int threads = 1000;
    const unsigned int vectorSize = sizeof(float) * vectorLength;

    pjVector_t firstVector = { (float *)calloc(vectorLength, sizeof(float)), vectorLength };
    pjVector_t secondVector = { (float *)calloc(vectorLength, sizeof(float)), vectorLength };
    pjVector_t resultVector = { (float *)calloc(vectorLength, sizeof(float)), vectorLength };

    float *d_firstVector;
    float *d_secondVector;
    float *d_resultVector;

    cudaMalloc((void **)&d_firstVector, vectorSize);
    cudaMalloc((void **)&d_secondVector, vectorSize);
    cudaMalloc((void **)&d_resultVector, vectorSize);

    for (unsigned int i = 0; i < vectorLength; i++)
    {
        firstVector.content[i] = 1.0f;
        secondVector.content[i] = 2.0f;
    }

    cudaMemcpy(d_firstVector, firstVector.content, vectorSize, cudaMemcpyHostToDevice);
    cudaMemcpy(d_secondVector, secondVector.content, vectorSize, cudaMemcpyHostToDevice);
    AddVectorsKernel<<<blocks, threads>>>(d_firstVector, d_secondVector, d_resultVector);
    cudaMemcpy(resultVector.content, d_resultVector, vectorSize, cudaMemcpyDeviceToHost);

    free(firstVector.content);
    free(secondVector.content);
    free(resultVector.content);

    cudaFree(d_firstVector);
    cudaFree(d_secondVector);
    cudaFree(d_resultVector);
    cudaDeviceReset();

    return 0;
}

First of all, I'm going to modify the number of threads per block to be as high as possible on my machine to decrease the number of used blocks. To do this, I'm calling cudaGetDeviceProperties function and reading the maxThreadsPerBlock field from the populated structure cudaDeviceProp. I also got rid of const keyword in for blocks and threads since they are populated dynamically.

#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <math.h>

#define VECTOR_LENGTH 1000000

#pragma comment(lib, "cudart") 

typedef struct 
{
    float *content;
    const unsigned int size;
} pjVector_t;

__global__ void AddVectorsKernel(float *firstVector, float *secondVector, float *resultVector)
{
    unsigned int index = threadIdx.x + blockIdx.x * blockDim.x;

    if(index < VECTOR_LENGTH)
    {
        resultVector[index] = firstVector[index] + secondVector[index];
    }
}

int main(void)
{
    const unsigned int vectorSize = sizeof(float) * VECTOR_LENGTH;
    int threads = 0;
    unsigned int blocks = 0;
    cudaDeviceProp deviceProperties;

    cudaGetDeviceProperties(&deviceProperties, 0);

    threads = deviceProperties.maxThreadsPerBlock;
    blocks = (unsigned int)ceil(VECTOR_LENGTH / (double)threads);

    pjVector_t firstVector = { (float *)calloc(VECTOR_LENGTH, sizeof(float)), VECTOR_LENGTH };
    pjVector_t secondVector = { (float *)calloc(VECTOR_LENGTH, sizeof(float)), VECTOR_LENGTH };
    pjVector_t resultVector = { (float *)calloc(VECTOR_LENGTH, sizeof(float)), VECTOR_LENGTH };

    float *d_firstVector;
    float *d_secondVector;
    float *d_resultVector;

    cudaMalloc((void **)&d_firstVector, vectorSize);
    cudaMalloc((void **)&d_secondVector, vectorSize);
    cudaMalloc((void **)&d_resultVector, vectorSize);

    for (unsigned int i = 0; i < VECTOR_LENGTH; i++)
    {
        firstVector.content[i] = 1.0f;
        secondVector.content[i] = 2.0f;
    }

    cudaMemcpy(d_firstVector, firstVector.content, vectorSize, cudaMemcpyHostToDevice);
    cudaMemcpy(d_secondVector, secondVector.content, vectorSize, cudaMemcpyHostToDevice);
    AddVectorsKernel<<<blocks, threads>>>(d_firstVector, d_secondVector, d_resultVector);
    cudaMemcpy(resultVector.content, d_resultVector, vectorSize, cudaMemcpyDeviceToHost);

    free(firstVector.content);
    free(secondVector.content);
    free(resultVector.content);

    cudaFree(d_firstVector);
    cudaFree(d_secondVector);
    cudaFree(d_resultVector);
    cudaDeviceReset();

    return 0;
}

After running this code, performance dropped to 833[µs] (~5%). 'Did I do something wrong?', I was asking myself. And actually I did - that 'if' statement in the kernel code was that 'something wrong'. Or maybe not entirely wrong, but not very well placed. I thought initially: how come now, since I'm utilising all possible threads in every block (well, maybe except the last one - it will have 448 unused threads, 977 * 1024 - 1000000 = 448), my kernel uses more time to do its job than previously, where the number of unused threads was a lot (~53 times!) higher? The answer was that 'if' - I quickly realised that every thread in every block had to evaluate it 1 million (+448) times. Why do that since I can precisely adjust the length of my vector (+ that 448 elements) so that 'if' is not needed anymore and I will never go outside my vector in the kernel?

I'm obviously just playing with the toolkit right now so it's not a problem to add something here and tweak something there to achieve my goal, but even in the real life scenario (which I highly doubt this code will ever be used for!) you could do the same thing - add some extra elements to you vector and fill them with zeros not to check if you're within your vector in the kernel code.

Let's do some (integer!) math:

• Unused threads in my previous approach: (1024 - 1000) * 1000 = 24000 (!)
  Redundant threads: 0 (all working threads are important)
  index checked: 0 times
• Unused threads in my current approach: ((1000000 / 1024) * 1024) - 1000000 + 1024 = -576 + 1024 = 448 (they are doing nothing)
  Redundant threads: 0 (all working threads are important)
  index checked: 1000448 times (!)
• Unused threads in my new approach: 0 (all are used)
  Redundant threads: 1000448 - 1000000 = 448 (they are just adding zeros)
  index checked: 0 times

Now let's see some code (I'm not filling these last 448 vector items with zeros, but you get the idea):

#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <math.h>

#define VECTOR_LENGTH 1000448

#pragma comment(lib, "cudart") 

typedef struct 
{
    float *content;
    const unsigned int size;
} pjVector_t;

__global__ void AddVectorsKernel(float *firstVector, float *secondVector, float *resultVector)
{
    unsigned int index = threadIdx.x + blockIdx.x * blockDim.x;
    resultVector[index] = firstVector[index] + secondVector[index];
}

int main(void)
{
    const unsigned int vectorSize = sizeof(float) * VECTOR_LENGTH;
    int threads = 0;
    unsigned int blocks = 0;
    cudaDeviceProp deviceProperties;

    cudaGetDeviceProperties(&deviceProperties, 0);

    threads = deviceProperties.maxThreadsPerBlock;
    blocks = (unsigned int)ceil(VECTOR_LENGTH / (double)threads);

    pjVector_t firstVector = { (float *)calloc(VECTOR_LENGTH, sizeof(float)), VECTOR_LENGTH };
    pjVector_t secondVector = { (float *)calloc(VECTOR_LENGTH, sizeof(float)), VECTOR_LENGTH };
    pjVector_t resultVector = { (float *)calloc(VECTOR_LENGTH, sizeof(float)), VECTOR_LENGTH };

    float *d_firstVector;
    float *d_secondVector;
    float *d_resultVector;

    cudaMalloc((void **)&d_firstVector, vectorSize);
    cudaMalloc((void **)&d_secondVector, vectorSize);
    cudaMalloc((void **)&d_resultVector, vectorSize);

    for (unsigned int i = 0; i < VECTOR_LENGTH; i++)
    {
        firstVector.content[i] = 1.0f;
        secondVector.content[i] = 2.0f;
    }

    cudaMemcpy(d_firstVector, firstVector.content, vectorSize, cudaMemcpyHostToDevice);
    cudaMemcpy(d_secondVector, secondVector.content, vectorSize, cudaMemcpyHostToDevice);
    AddVectorsKernel<<<blocks, threads>>>(d_firstVector, d_secondVector, d_resultVector);
    cudaMemcpy(resultVector.content, d_resultVector, vectorSize, cudaMemcpyDeviceToHost);

    free(firstVector.content);
    free(secondVector.content);
    free(resultVector.content);

    cudaFree(d_firstVector);
    cudaFree(d_secondVector);
    cudaFree(d_resultVector);
    cudaDeviceReset();

    return 0;
}

And there it is: 756.479[µs] - time needed to process more data was actually ~5% shorter!



That's it for now, thanks for reading, I hope it was useful or at least interesting.

CUDA 5.0 first calculations

OK - now, since I have my 'Hello World' program done, I think it would make sense to write something that actually uses CUDA processing power. Inspired by the CUDA introductory video, I want to write a simple vector addition kernel and run a quick test on how it compares to the CPU processing. I haven't mentioned it earlier, but I'm writing and running this code on my GeForce GTX 660M laptop GPU. Let's get started!

My first piece of code runs on the CPU and its only purpose is to add two vectors (of size one million, 1,000,000) together.

#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <time.h>

#pragma comment(lib, "cudart") 

typedef struct 
{
 float * const content;
 const unsigned int size;
} pjVector_t;

void AddVectors(const pjVector_t * const firstVector, const pjVector_t * const secondVector, pjVector_t * const resultVector);

int main(void)
{
 unsigned int vectorSize = 1000000;
 double cpuTime;
 clock_t begin, end;
 pjVector_t firstVector = { (float *)calloc(vectorSize, sizeof(float)), vectorSize };
 pjVector_t secondVector = { (float *)calloc(vectorSize, sizeof(float)), vectorSize };
 pjVector_t resultVector = { (float *)calloc(vectorSize, sizeof(float)), vectorSize };

 for (unsigned int i = 0; i < vectorSize; i++)
 {
  firstVector.content[i] = 1.0f;
  secondVector.content[i] = 2.0f;
 }

 begin = clock();

 AddVectors(&firstVector, &secondVector, &resultVector);

 end = clock();
 cpuTime = (double)(end - begin) / CLOCKS_PER_SEC;

 printf("Result vector calculated in: %f[sec]\n", cpuTime);
 getchar();


 free(firstVector.content);
 free(secondVector.content);
 free(resultVector.content);
 
 return 0;
}

void AddVectors(const pjVector_t * const firstVector, const pjVector_t * const secondVector, pjVector_t * const resultVector)
{
 for (unsigned int i = 0; i < firstVector -> size; i++)
 {
  resultVector -> content[i] = firstVector -> content[i] + secondVector -> content[i];
 }
}

Even though it's single-threaded it's pretty fast and takes from 5 to 6 milliseconds to calculate the result vector.

My second piece of code combines GPU and CPU code to compare the performance. A couple of things before I show the code though:

• To simplify passing parameters to my kernel, I ditched the structures and chose plain float arrays.
• Although I'm sure it can be done a lot more efficient, I'm running my CUDA code using 1000 blocks and 1000 threads each (1,000 * 1,000 = 1,000,000) just to show how the problem is divided into sections.
• Performance is measured by an external tool, NVIDIA Visual Profiler which comes with the CUDA Toolkit.

#include <cuda.h>
#include <cuda_runtime.h>
#include <device_launch_parameters.h>
#include <stdio.h>
#include <time.h>

#pragma comment(lib, "cudart") 

typedef struct 
{
 float *content;
 const unsigned int size;
} pjVector_t;

__global__ void AddVectorsKernel(float *firstVector, float *secondVector, float *resultVector)
{
 unsigned int index = threadIdx.x + blockIdx.x * blockDim.x;
 resultVector[index] = firstVector[index] + secondVector[index];
}

void AddVectors(const pjVector_t * const firstVector, const pjVector_t * const secondVector, pjVector_t * const resultVector);

int main(void)
{
 const unsigned int vectorLength = 1000000;
 const unsigned int blocks = 1000;
 const unsigned int threads = 1000;
 const unsigned int vectorSize = sizeof(float) * vectorLength;
 double cpuTime;
 clock_t begin, end;

 pjVector_t firstVector = { (float *)calloc(vectorLength, sizeof(float)), vectorLength };
 pjVector_t secondVector = { (float *)calloc(vectorLength, sizeof(float)), vectorLength };
 pjVector_t resultVector = { (float *)calloc(vectorLength, sizeof(float)), vectorLength };

 float *d_firstVector;
 float *d_secondVector;
 float *d_resultVector;

 cudaMalloc((void **)&d_firstVector, vectorSize);
 cudaMalloc((void **)&d_secondVector, vectorSize);
 cudaMalloc((void **)&d_resultVector, vectorSize);

 for (unsigned int i = 0; i < vectorLength; i++)
 {
  firstVector.content[i] = 1.0f;
  secondVector.content[i] = 2.0f;
 }

 // CPU calculatons
 begin = clock();

 AddVectors(&firstVector, &secondVector, &resultVector);

 end = clock();
 cpuTime = (double)(end - begin) / CLOCKS_PER_SEC;
 // - CPU calculatons

 // GPU calculatons
 cudaMemcpy(d_firstVector, firstVector.content, vectorSize, cudaMemcpyHostToDevice);
 cudaMemcpy(d_secondVector, secondVector.content, vectorSize, cudaMemcpyHostToDevice);

 AddVectorsKernel<<<blocks, threads>>>(d_firstVector, d_secondVector, d_resultVector);

 cudaMemcpy(resultVector.content, d_resultVector, vectorSize, cudaMemcpyDeviceToHost);
 // - GPU calculatons

 free(firstVector.content);
 free(secondVector.content);
 free(resultVector.content);

 cudaFree(d_firstVector);
 cudaFree(d_secondVector);
 cudaFree(d_resultVector);
 cudaDeviceReset();

 printf("CPU result vector calculated in: %f[ms]\n", cpuTime * 1000.0);

 getchar();
 
 return 0;
}

void AddVectors(const pjVector_t * const firstVector, const pjVector_t * const secondVector, pjVector_t * const resultVector)
{
 for (unsigned int i = 0; i < firstVector -> size; i++)
 {
  resultVector -> content[i] = firstVector -> content[i] + secondVector -> content[i];
 }
}

This time the magic really happened: vector addition takes only 792.166[µs], which is roughly 15% of the previous value. That's simply amazing for the first try. Let me show some screenshots from the profiler:



But what really made me smile was the warning messages from the profiler:



It actually complains about the time needed to copy the data being longer than the time of the calculations! And see how much time the cudaMalloc needed? As it turns out, adding two vectors of size 1,000,000 is the smallest problem here! I am really pleased with my first real CUDA test and it certainly inspires me even more.