Day 34 - GPU Kernels Learning Journey

Description

### Files: `vec_add_cpu.c`, `vec_add_gpu.cu`, `benchmark.py` **Summary:** Today I compared the performance of vector addition implemented using two different approaches: CPU with MPI and Loop Unrolling - Parallelized across multiple CPU cores using MPI. GPU with CUDA - Executed on a GPU for massive parallelism. The script runs performance tests on different input sizes and generates a log-log scale plot comparing execution times.

Code Files

benchmark.py

import subprocess
import numpy as np
import matplotlib.pyplot as plt

# Define test sizes
sizes = [10**i for i in range(1, 10)]  # 10, 100, 1000, ..., 1 billion
cpu_times = []
gpu_times = []

# Compile CPU and GPU programs
subprocess.run("mpicc -o vec_add_cpu vec_add_cpu.c -O3", shell=True, check=True)
subprocess.run("nvcc -o vec_add_gpu vec_add_gpu.cu -O3", shell=True, check=True)

for size in sizes:
    print(f"Running for size: {size}")
    
    # Run MPI CPU version
    result_cpu = subprocess.run(
        f"mpirun --oversubscribe --allow-run-as-root -np 4 ./vec_add_cpu {size}",
        shell=True, capture_output=True, text=True
    )
    
    # Debugging output
    print("CPU Output:", result_cpu.stdout)
    print("Error (if any):", result_cpu.stderr)
    
    output = result_cpu.stdout.strip().split()
    if len(output) >= 2:
        cpu_time = float(output[-2])  # Extract time if available
    else:
        print("Error: Unexpected CPU output format ->", result_cpu.stdout)
        cpu_time = float('inf')  # Assign a large value to indicate failure
    cpu_times.append(cpu_time)
    
    # Run CUDA GPU version
    result_gpu = subprocess.run(
        f"./vec_add_gpu {size}", shell=True, capture_output=True, text=True
    )
    
    # Debugging output
    print("GPU Output:", result_gpu.stdout)
    print("Error (if any):", result_gpu.stderr)
    
    output = result_gpu.stdout.strip().split()
    if len(output) >= 2:
        gpu_time = float(output[-2])  # Extract time if available
    else:
        print("Error: Unexpected GPU output format ->", result_gpu.stdout)
        gpu_time = float('inf')  # Assign a large value to indicate failure
    gpu_times.append(gpu_time)

# Plot results
plt.figure(figsize=(10, 6))
plt.plot(sizes, cpu_times, marker='o', label='CPU (MPI + Unrolling)')
plt.plot(sizes, gpu_times, marker='s', label='GPU (CUDA)')
plt.xscale('log')
plt.yscale('log')
plt.xlabel("Input Size (log scale)")
plt.ylabel("Execution Time (ms, log scale)")
plt.legend()
plt.grid()
plt.title("Performance Comparison: CPU vs. GPU for Vector Addition")
plt.show()

vec_add_gpu.cu

#include <stdio.h>
#include <cuda.h>

#define UNROLL_FACTOR 4

__global__ void vector_add(double *A, double *B, double *C, int N) {
    int i = blockIdx.x * blockDim.x + threadIdx.x;
    int limit = N - (N % UNROLL_FACTOR);

    for (; i < limit; i += blockDim.x * gridDim.x) {
        C[i] = A[i] + B[i];
        C[i+1] = A[i+1] + B[i+1];
        C[i+2] = A[i+2] + B[i+2];
        C[i+3] = A[i+3] + B[i+3];
    }

    for (; i < N; i++) {
        C[i] = A[i] + B[i];
    }
}

int main(int argc, char **argv) {
    int N = atoi(argv[1]);
    double *h_A, *h_B, *h_C;
    double *d_A, *d_B, *d_C;

    size_t size = N * sizeof(double);
    h_A = (double*) malloc(size);
    h_B = (double*) malloc(size);
    h_C = (double*) malloc(size);

    for (int i = 0; i < N; i++) {
        h_A[i] = i;
        h_B[i] = i * 2;
    }

    cudaMalloc((void**)&d_A, size);
    cudaMalloc((void**)&d_B, size);
    cudaMalloc((void**)&d_C, size);

    cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);

    int threadsPerBlock = 256;
    int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;

    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    cudaEventRecord(start);

    vector_add<<<blocksPerGrid, threadsPerBlock>>>(d_A, d_B, d_C, N);
    
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);

    cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);

    float elapsedTime;
    cudaEventElapsedTime(&elapsedTime, start, stop);
    printf("GPU (CUDA): %f ms\n", elapsedTime);

    cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
    free(h_A); free(h_B); free(h_C);

    return 0;
}

vec_add_cpu.c

#include <stdio.h>
#include <stdlib.h>
#include <mpi.h>
#include <time.h>

#define UNROLL_FACTOR 4  // Loop unrolling factor

void vector_add(double *A, double *B, double *C, int size) {
    int i, limit = size - (size % UNROLL_FACTOR);
    for (i = 0; i < limit; i += UNROLL_FACTOR) {
        C[i]   = A[i]   + B[i];
        C[i+1] = A[i+1] + B[i+1];
        C[i+2] = A[i+2] + B[i+2];
        C[i+3] = A[i+3] + B[i+3];
    }
    for (; i < size; i++) {
        C[i] = A[i] + B[i];
    }
}

int main(int argc, char **argv) {
    int rank, size, N = atoi(argv[1]);  // Take vector size as argument
    MPI_Init(&argc, &argv);
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    MPI_Comm_size(MPI_COMM_WORLD, &size);

    int chunk_size = N / size;
    int start = rank * chunk_size;
    int end = (rank == size - 1) ? N : start + chunk_size;

    double *A = (double*) malloc(chunk_size * sizeof(double));
    double *B = (double*) malloc(chunk_size * sizeof(double));
    double *C = (double*) malloc(chunk_size * sizeof(double));

    for (int i = 0; i < chunk_size; i++) {
        A[i] = i + rank;
        B[i] = i - rank;
    }

    double t1 = MPI_Wtime();
    vector_add(A, B, C, chunk_size);
    double t2 = MPI_Wtime();

    if (rank == 0) {
        printf("CPU (MPI + Unrolling): %lf miliseconds\n", (t2 - t1)*1000.0);
    }

    free(A); free(B); free(C);
    MPI_Finalize();
    return 0;
}