Graphics Processing Units (GPUs) have evolved from specialized graphics hardware to general-purpose parallel processors capable of executing thousands of threads simultaneously. GPU programming enables orders-of-magnitude speedups for data-parallel workloads. This comprehensive guide covers the fundamentals of GPU programming in C using NVIDIA CUDA and the cross-platform OpenCL framework.

Understanding GPU Architecture

GPUs are designed for massive parallelism with thousands of cores optimized for throughput rather than latency.

┌─────────────────────────────────────────────────────────────────┐
│                          CPU Host                               │
│  ┌─────────┐  ┌─────────┐  ┌─────────┐  ┌─────────┐            │
│  │ Core 0  │  │ Core 1  │  │ Core 2  │  │ Core 3  │            │
│  └─────────┘  └─────────┘  └─────────┘  └─────────┘            │
│                         │                                       │
│                   PCIe Bus                                       │
│                         │                                       │
├─────────────────────────┼───────────────────────────────────────┤
│                         ▼                                       │
│                        GPU Device                               │
│  ┌─────────────────────────────────────────────────────────┐   │
│  │                   Global Memory                          │   │
│  └─────────────────────────────────────────────────────────┘   │
│         │              │              │              │          │
│         ▼              ▼              ▼              ▼          │
│  ┌──────────┐   ┌──────────┐   ┌──────────┐   ┌──────────┐    │
│  │  SM 0    │   │  SM 1    │   │  SM 2    │   │  SM 3    │    │
│  │ ┌──────┐ │   │ ┌──────┐ │   │ ┌──────┐ │   │ ┌──────┐ │    │
│  │ │Cores │ │   │ │Cores │ │   │ │Cores │ │   │ │Cores │ │    │
│  │ │128   │ │   │ │128   │ │   │ │128   │ │   │ │128   │ │    │
│  │ └──────┘ │   │ └──────┘ │   │ └──────┘ │   │ └──────┘ │    │
│  └──────────┘   └──────────┘   └──────────┘   └──────────┘    │
└─────────────────────────────────────────────────────────────────┘

Part 1: CUDA Programming

CUDA (Compute Unified Device Architecture) is NVIDIA's parallel computing platform.

1. CUDA Setup

# Install CUDA Toolkit
# Download from developer.nvidia.com/cuda-downloads
# Verify installation
nvcc --version
nvidia-smi
# Check CUDA samples
ls /usr/local/cuda/samples

2. First CUDA Program

// hello_cuda.cu
#include <stdio.h>
#include <cuda_runtime.h>
// Kernel function (runs on GPU)
__global__ void hello_kernel() {
int thread_id = threadIdx.x + blockIdx.x * blockDim.x;
printf("Hello from GPU thread %d\n", thread_id);
}
int main() {
// Launch kernel with 16 threads
hello_kernel<<<1, 16>>>();
// Wait for GPU to finish
cudaDeviceSynchronize();
printf("Hello from CPU\n");
return 0;
}

Compile and run:

nvcc -o hello_cuda hello_cuda.cu
./hello_cuda

3. Vector Addition (The "Hello World" of GPU Computing)

// vector_add.cu
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <time.h>
#define N 1000000
#define BLOCK_SIZE 256
// CUDA kernel for vector addition
__global__ void vector_add(float *a, float *b, float *c, int n) {
int index = blockIdx.x * blockDim.x + threadIdx.x;
if (index < n) {
c[index] = a[index] + b[index];
}
}
int main() {
float *h_a, *h_b, *h_c;     // Host arrays
float *d_a, *d_b, *d_c;     // Device arrays
size_t size = N * sizeof(float);
// Allocate host memory
h_a = (float*)malloc(size);
h_b = (float*)malloc(size);
h_c = (float*)malloc(size);
// Initialize host arrays
for (int i = 0; i < N; i++) {
h_a[i] = rand() / (float)RAND_MAX;
h_b[i] = rand() / (float)RAND_MAX;
}
// Allocate device memory
cudaMalloc(&d_a, size);
cudaMalloc(&d_b, size);
cudaMalloc(&d_c, size);
// Copy data from host to device
cudaMemcpy(d_a, h_a, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_b, h_b, size, cudaMemcpyHostToDevice);
// Launch kernel
int blocks = (N + BLOCK_SIZE - 1) / BLOCK_SIZE;
vector_add<<<blocks, BLOCK_SIZE>>>(d_a, d_b, d_c, N);
// Wait for kernel to finish
cudaDeviceSynchronize();
// Copy result back to host
cudaMemcpy(h_c, d_c, size, cudaMemcpyDeviceToHost);
// Verify result
int errors = 0;
for (int i = 0; i < N; i++) {
if (fabs(h_c[i] - (h_a[i] + h_b[i])) > 1e-5) {
errors++;
}
}
printf("Vector addition completed with %d errors\n", errors);
// Cleanup
free(h_a); free(h_b); free(h_c);
cudaFree(d_a); cudaFree(d_b); cudaFree(d_c);
return 0;
}

4. Matrix Multiplication

// matrix_mul.cu
#include <stdio.h>
#include <stdlib.h>
#include <cuda_runtime.h>
#include <math.h>
#define TILE_SIZE 16
// Naive matrix multiplication kernel
__global__ void matmul_naive(float *A, float *B, float *C, int N) {
int row = blockIdx.y * blockDim.y + threadIdx.y;
int col = blockIdx.x * blockDim.x + threadIdx.x;
if (row < N && col < N) {
float sum = 0.0f;
for (int k = 0; k < N; k++) {
sum += A[row * N + k] * B[k * N + col];
}
C[row * N + col] = sum;
}
}
// Tiled matrix multiplication kernel (optimized)
__global__ void matmul_tiled(float *A, float *B, float *C, int N) {
__shared__ float As[TILE_SIZE][TILE_SIZE];
__shared__ float Bs[TILE_SIZE][TILE_SIZE];
int bx = blockIdx.x;
int by = blockIdx.y;
int tx = threadIdx.x;
int ty = threadIdx.y;
int row = by * TILE_SIZE + ty;
int col = bx * TILE_SIZE + tx;
float sum = 0.0f;
for (int k = 0; k < (N + TILE_SIZE - 1) / TILE_SIZE; k++) {
// Load tile from A
if (row < N && k * TILE_SIZE + tx < N) {
As[ty][tx] = A[row * N + k * TILE_SIZE + tx];
} else {
As[ty][tx] = 0.0f;
}
// Load tile from B
if (col < N && k * TILE_SIZE + ty < N) {
Bs[ty][tx] = B[(k * TILE_SIZE + ty) * N + col];
} else {
Bs[ty][tx] = 0.0f;
}
__syncthreads();
// Compute partial sum
for (int i = 0; i < TILE_SIZE; i++) {
sum += As[ty][i] * Bs[i][tx];
}
__syncthreads();
}
if (row < N && col < N) {
C[row * N + col] = sum;
}
}
// Host function to initialize matrix
void init_matrix(float *mat, int N) {
for (int i = 0; i < N * N; i++) {
mat[i] = rand() / (float)RAND_MAX;
}
}
// Host function to verify result
float verify_matrix(float *C, float *A, float *B, int N) {
float max_error = 0.0f;
for (int i = 0; i < N; i++) {
for (int j = 0; j < N; j++) {
float expected = 0.0f;
for (int k = 0; k < N; k++) {
expected += A[i * N + k] * B[k * N + j];
}
float error = fabs(C[i * N + j] - expected);
if (error > max_error) max_error = error;
}
}
return max_error;
}
int main() {
int N = 1024;
size_t size = N * N * sizeof(float);
// Allocate host memory
float *h_A = (float*)malloc(size);
float *h_B = (float*)malloc(size);
float *h_C = (float*)malloc(size);
// Initialize matrices
srand(time(NULL));
init_matrix(h_A, N);
init_matrix(h_B, N);
// Allocate device memory
float *d_A, *d_B, *d_C;
cudaMalloc(&d_A, size);
cudaMalloc(&d_B, size);
cudaMalloc(&d_C, size);
// Copy to device
cudaMemcpy(d_A, h_A, size, cudaMemcpyHostToDevice);
cudaMemcpy(d_B, h_B, size, cudaMemcpyHostToDevice);
// Configure kernel launch
dim3 blockSize(TILE_SIZE, TILE_SIZE);
dim3 gridSize((N + TILE_SIZE - 1) / TILE_SIZE,
(N + TILE_SIZE - 1) / TILE_SIZE);
// Measure time for tiled kernel
cudaEvent_t start, stop;
cudaEventCreate(&start);
cudaEventCreate(&stop);
cudaEventRecord(start);
matmul_tiled<<<gridSize, blockSize>>>(d_A, d_B, d_C, N);
cudaEventRecord(stop);
cudaEventSynchronize(stop);
float milliseconds = 0;
cudaEventElapsedTime(&milliseconds, start, stop);
// Copy result back
cudaMemcpy(h_C, d_C, size, cudaMemcpyDeviceToHost);
// Verify result
float max_error = verify_matrix(h_C, h_A, h_B, N);
printf("Matrix multiplication (N=%d)\n", N);
printf("  Time: %.2f ms\n", milliseconds);
printf("  Max error: %e\n", max_error);
printf("  Performance: %.2f GFLOPS\n",
2.0 * N * N * N / (milliseconds / 1000) / 1e9);
// Cleanup
free(h_A); free(h_B); free(h_C);
cudaFree(d_A); cudaFree(d_B); cudaFree(d_C);
cudaEventDestroy(start);
cudaEventDestroy(stop);
return 0;
}

5. CUDA Memory Types

// memory_types.cu
#include <stdio.h>
#include <cuda_runtime.h>
// Global memory (global scope)
__device__ int global_counter = 0;
// Shared memory example
__global__ void shared_memory_example(int *input, int *output, int n) {
// Shared memory allocated per block
__shared__ int shared_data[256];
int idx = threadIdx.x + blockIdx.x * blockDim.x;
int tid = threadIdx.x;
// Load data into shared memory
if (idx < n) {
shared_data[tid] = input[idx];
}
__syncthreads();
// Process using shared memory
if (tid < n) {
int sum = 0;
for (int i = 0; i < blockDim.x; i++) {
sum += shared_data[i];
}
output[idx] = sum;
}
}
// Constant memory (cached, read-only)
__constant__ float constant_array[256];
__global__ void constant_memory_example(float *output) {
int idx = threadIdx.x + blockIdx.x * blockDim.x;
output[idx] = constant_array[idx % 256];
}
// Texture memory (cached, read-only)
texture<float, 1, cudaReadModeElementType> texture_data;
__global__ void texture_memory_example(float *output, int n) {
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx < n) {
output[idx] = tex1Dfetch(texture_data, idx);
}
}
int main() {
const int N = 1024;
size_t size = N * sizeof(float);
// Allocate host memory
float *h_data = (float*)malloc(size);
for (int i = 0; i < N; i++) {
h_data[i] = i * 1.0f;
}
// Constant memory example
cudaMemcpyToSymbol(constant_array, h_data, size);
// Texture memory example
float *d_texture;
cudaMalloc(&d_texture, size);
cudaMemcpy(d_texture, h_data, size, cudaMemcpyHostToDevice);
cudaBindTexture(0, texture_data, d_texture, size);
// Launch kernels...
// Cleanup
cudaUnbindTexture(texture_data);
cudaFree(d_texture);
free(h_data);
return 0;
}

6. CUDA Streams for Concurrency

// streams.cu
#include <stdio.h>
#include <cuda_runtime.h>
#define N 10000000
#define NUM_STREAMS 4
__global__ void process_data(float *data, int n) {
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx < n) {
data[idx] = sqrtf(data[idx] * data[idx]);
}
}
int main() {
// Create CUDA streams
cudaStream_t streams[NUM_STREAMS];
for (int i = 0; i < NUM_STREAMS; i++) {
cudaStreamCreate(&streams[i]);
}
size_t size = N * sizeof(float);
size_t stream_size = size / NUM_STREAMS;
// Allocate pinned host memory
float *h_data;
cudaHostAlloc(&h_data, size, cudaHostAllocDefault);
// Initialize host data
for (int i = 0; i < N; i++) {
h_data[i] = i * 1.0f;
}
// Allocate device memory
float *d_data;
cudaMalloc(&d_data, size);
// Launch kernels in streams (overlapping computation and transfer)
for (int i = 0; i < NUM_STREAMS; i++) {
int offset = i * stream_size / sizeof(float);
// Async memory copy
cudaMemcpyAsync(d_data + offset, h_data + offset,
stream_size, cudaMemcpyHostToDevice, streams[i]);
// Kernel launch
int blocks = (stream_size / sizeof(float) + 255) / 256;
process_data<<<blocks, 256, 0, streams[i]>>>(d_data + offset,
stream_size / sizeof(float));
// Async copy back
cudaMemcpyAsync(h_data + offset, d_data + offset,
stream_size, cudaMemcpyDeviceToHost, streams[i]);
}
// Synchronize all streams
for (int i = 0; i < NUM_STREAMS; i++) {
cudaStreamSynchronize(streams[i]);
}
// Cleanup
for (int i = 0; i < NUM_STREAMS; i++) {
cudaStreamDestroy(streams[i]);
}
cudaFreeHost(h_data);
cudaFree(d_data);
return 0;
}

Part 2: OpenCL Programming

OpenCL (Open Computing Language) is a cross-platform framework for heterogeneous computing.

1. OpenCL Setup

# Install OpenCL headers and libraries
sudo apt-get install opencl-headers ocl-icd-opencl-dev
# Check OpenCL platforms
clinfo

2. Basic OpenCL Program

// vector_add_opencl.c
#include <stdio.h>
#include <stdlib.h>
#include <CL/cl.h>
#define N 1000000
#define KERNEL_SOURCE \
"__kernel void vector_add(__global float *a, __global float *b, __global float *c, int n) {\n" \
"    int i = get_global_id(0);\n" \
"    if (i < n) {\n" \
"        c[i] = a[i] + b[i];\n" \
"    }\n" \
"}\n"
int main() {
cl_platform_id platform;
cl_device_id device;
cl_context context;
cl_command_queue queue;
cl_program program;
cl_kernel kernel;
cl_mem d_a, d_b, d_c;
float *h_a, *h_b, *h_c;
size_t size = N * sizeof(float);
// Get platform
clGetPlatformIDs(1, &platform, NULL);
// Get GPU device
clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);
// Create context
context = clCreateContext(NULL, 1, &device, NULL, NULL, NULL);
// Create command queue
queue = clCreateCommandQueue(context, device, 0, NULL);
// Allocate host memory
h_a = (float*)malloc(size);
h_b = (float*)malloc(size);
h_c = (float*)malloc(size);
// Initialize arrays
for (int i = 0; i < N; i++) {
h_a[i] = i * 1.0f;
h_b[i] = i * 2.0f;
}
// Allocate device memory
d_a = clCreateBuffer(context, CL_MEM_READ_ONLY, size, NULL, NULL);
d_b = clCreateBuffer(context, CL_MEM_READ_ONLY, size, NULL, NULL);
d_c = clCreateBuffer(context, CL_MEM_WRITE_ONLY, size, NULL, NULL);
// Copy data to device
clEnqueueWriteBuffer(queue, d_a, CL_TRUE, 0, size, h_a, 0, NULL, NULL);
clEnqueueWriteBuffer(queue, d_b, CL_TRUE, 0, size, h_b, 0, NULL, NULL);
// Create program from source
program = clCreateProgramWithSource(context, 1, (const char**)&KERNEL_SOURCE, NULL, NULL);
// Build program
clBuildProgram(program, 1, &device, NULL, NULL, NULL);
// Create kernel
kernel = clCreateKernel(program, "vector_add", NULL);
// Set kernel arguments
clSetKernelArg(kernel, 0, sizeof(cl_mem), &d_a);
clSetKernelArg(kernel, 1, sizeof(cl_mem), &d_b);
clSetKernelArg(kernel, 2, sizeof(cl_mem), &d_c);
clSetKernelArg(kernel, 3, sizeof(int), &N);
// Execute kernel
size_t global_size = N;
clEnqueueNDRangeKernel(queue, kernel, 1, NULL, &global_size, NULL, 0, NULL, NULL);
// Read result
clEnqueueReadBuffer(queue, d_c, CL_TRUE, 0, size, h_c, 0, NULL, NULL);
// Verify result
int errors = 0;
for (int i = 0; i < N; i++) {
if (h_c[i] != h_a[i] + h_b[i]) errors++;
}
printf("Vector addition completed with %d errors\n", errors);
// Cleanup
free(h_a); free(h_b); free(h_c);
clReleaseMemObject(d_a);
clReleaseMemObject(d_b);
clReleaseMemObject(d_c);
clReleaseKernel(kernel);
clReleaseProgram(program);
clReleaseCommandQueue(queue);
clReleaseContext(context);
return 0;
}

3. OpenCL with Kernel File

// kernel.cl - OpenCL kernel file
__kernel void matmul(
__global const float *A,
__global const float *B,
__global float *C,
int N) {
int row = get_global_id(1);
int col = get_global_id(0);
if (row < N && col < N) {
float sum = 0.0f;
for (int k = 0; k < N; k++) {
sum += A[row * N + k] * B[k * N + col];
}
C[row * N + col] = sum;
}
}

// opencl_matmul.c
#include <stdio.h>
#include <stdlib.h>
#include <CL/cl.h>
#include <time.h>
#define N 1024
// Load kernel from file
char* load_kernel_source(const char *filename) {
FILE *fp = fopen(filename, "r");
if (!fp) return NULL;
fseek(fp, 0, SEEK_END);
long size = ftell(fp);
rewind(fp);
char *source = (char*)malloc(size + 1);
fread(source, 1, size, fp);
source[size] = '\0';
fclose(fp);
return source;
}
int main() {
cl_platform_id platform;
cl_device_id device;
cl_context context;
cl_command_queue queue;
cl_program program;
cl_kernel kernel;
cl_mem d_A, d_B, d_C;
float *h_A, *h_B, *h_C;
size_t size = N * N * sizeof(float);
// Get platform and device
clGetPlatformIDs(1, &platform, NULL);
clGetDeviceIDs(platform, CL_DEVICE_TYPE_GPU, 1, &device, NULL);
// Create context and queue
context = clCreateContext(NULL, 1, &device, NULL, NULL, NULL);
queue = clCreateCommandQueue(context, device, 0, NULL);
// Allocate host memory
h_A = (float*)malloc(size);
h_B = (float*)malloc(size);
h_C = (float*)malloc(size);
// Initialize matrices
srand(time(NULL));
for (int i = 0; i < N * N; i++) {
h_A[i] = rand() / (float)RAND_MAX;
h_B[i] = rand() / (float)RAND_MAX;
}
// Allocate device memory
d_A = clCreateBuffer(context, CL_MEM_READ_ONLY, size, NULL, NULL);
d_B = clCreateBuffer(context, CL_MEM_READ_ONLY, size, NULL, NULL);
d_C = clCreateBuffer(context, CL_MEM_WRITE_ONLY, size, NULL, NULL);
// Copy to device
clEnqueueWriteBuffer(queue, d_A, CL_TRUE, 0, size, h_A, 0, NULL, NULL);
clEnqueueWriteBuffer(queue, d_B, CL_TRUE, 0, size, h_B, 0, NULL, NULL);
// Load and build kernel
char *source = load_kernel_source("kernel.cl");
program = clCreateProgramWithSource(context, 1, (const char**)&source, NULL, NULL);
free(source);
clBuildProgram(program, 1, &device, NULL, NULL, NULL);
kernel = clCreateKernel(program, "matmul", NULL);
// Set kernel arguments
clSetKernelArg(kernel, 0, sizeof(cl_mem), &d_A);
clSetKernelArg(kernel, 1, sizeof(cl_mem), &d_B);
clSetKernelArg(kernel, 2, sizeof(cl_mem), &d_C);
clSetKernelArg(kernel, 3, sizeof(int), &N);
// Measure time
cl_event event;
size_t global_size[2] = {N, N};
clEnqueueNDRangeKernel(queue, kernel, 2, NULL, global_size, NULL,
0, NULL, &event);
clWaitForEvents(1, &event);
cl_ulong start, end;
clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_START,
sizeof(cl_ulong), &start, NULL);
clGetEventProfilingInfo(event, CL_PROFILING_COMMAND_END,
sizeof(cl_ulong), &end, NULL);
double time_ms = (end - start) / 1000000.0;
// Read result
clEnqueueReadBuffer(queue, d_C, CL_TRUE, 0, size, h_C, 0, NULL, NULL);
printf("Matrix multiplication (N=%d) completed in %.2f ms\n", N, time_ms);
printf("Performance: %.2f GFLOPS\n",
2.0 * N * N * N / (time_ms / 1000) / 1e9);
// Cleanup
free(h_A); free(h_B); free(h_C);
clReleaseMemObject(d_A);
clReleaseMemObject(d_B);
clReleaseMemObject(d_C);
clReleaseKernel(kernel);
clReleaseProgram(program);
clReleaseCommandQueue(queue);
clReleaseContext(context);
return 0;
}

Advanced GPU Optimization Techniques

1. Warp-Level Primitives

// warp_shuffle.cu
#include <stdio.h>
#include <cuda_runtime.h>
// Warp shuffle for reduction
__device__ float warp_reduce_sum(float val) {
for (int offset = warpSize / 2; offset > 0; offset >>= 1) {
val += __shfl_down_sync(0xFFFFFFFF, val, offset);
}
return val;
}
__global__ void reduction_shuffle(float *input, float *output, int n) {
float sum = 0.0f;
// Load data
for (int i = threadIdx.x; i < n; i += blockDim.x) {
sum += input[i];
}
// Warp-level reduction
sum = warp_reduce_sum(sum);
// Final reduction across warps
__shared__ float shared[32];
int warp_id = threadIdx.x / warpSize;
int lane_id = threadIdx.x % warpSize;
if (lane_id == 0) {
shared[warp_id] = sum;
}
__syncthreads();
if (threadIdx.x < 32) {
sum = shared[threadIdx.x];
sum = warp_reduce_sum(sum);
}
if (threadIdx.x == 0) {
output[blockIdx.x] = sum;
}
}

2. Dynamic Parallelism

// dynamic_parallelism.cu
#include <stdio.h>
#include <cuda_runtime.h>
// Child kernel for recursive decomposition
__global__ void child_kernel(float *data, int n) {
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx < n) {
data[idx] = sqrtf(data[idx]);
}
}
// Parent kernel that launches child kernels
__global__ void parent_kernel(float *data, int n) {
if (blockIdx.x == 0 && threadIdx.x == 0) {
// Dynamically launch child kernel
int blocks = (n + 255) / 256;
child_kernel<<<blocks, 256>>>(data, n);
}
}
int main() {
// Enable dynamic parallelism
cudaDeviceSetLimit(cudaLimitDevRuntimeSyncDepth, 10);
const int N = 1000000;
size_t size = N * sizeof(float);
float *d_data;
cudaMalloc(&d_data, size);
// Launch parent kernel
parent_kernel<<<1, 1>>>(d_data, N);
cudaDeviceSynchronize();
cudaFree(d_data);
return 0;
}

3. Unified Memory

// unified_memory.cu
#include <stdio.h>
#include <cuda_runtime.h>
__global__ void unified_memory_kernel(float *data, int n) {
int idx = threadIdx.x + blockIdx.x * blockDim.x;
if (idx < n) {
data[idx] = sqrtf(data[idx]);
}
}
int main() {
const int N = 10000000;
size_t size = N * sizeof(float);
// Allocate unified memory (accessible by both CPU and GPU)
float *data;
cudaMallocManaged(&data, size);
// Initialize on CPU
for (int i = 0; i < N; i++) {
data[i] = i * 1.0f;
}
// Launch kernel - data automatically migrates to GPU
int blocks = (N + 255) / 256;
unified_memory_kernel<<<blocks, 256>>>(data, N);
// Synchronize to ensure kernel completes
cudaDeviceSynchronize();
// Data now available on CPU again
printf("Result[0]: %f\n", data[0]);
// Cleanup
cudaFree(data);
return 0;
}

Performance Profiling

# NVIDIA profiling tools
nvprof ./vector_add
nvprof --metrics gld_efficiency,gst_efficiency ./matrix_mul
nvprof --events shared_load,shared_store ./matrix_mul
# Visual profiler
nvvp
# Compute profiler
nvprof --print-gpu-trace ./application

Common Optimization Patterns

Best Practices

1. Minimize Host-Device Transfers: Transfer data in bulk, not per-element
2. Use Asynchronous Operations: Overlap computation with transfers
3. Optimize Memory Access: Ensure coalesced access patterns
4. Maximize Occupancy: Balance register usage and block size
5. Use Shared Memory: Cache frequently accessed data
6. Avoid Divergence: Keep warp threads on same execution path
7. Profile Early: Identify bottlenecks with profiling tools
8. Check Error Codes: Always check CUDA/OpenCL return values
9. Use Appropriate Data Types: float for throughput, double for precision
10. Consider Power Limits: Balance performance with power consumption

Conclusion

GPU programming in C opens up possibilities for massive parallel computation. Key takeaways:

• CUDA provides NVIDIA-specific, high-performance access
• OpenCL offers cross-platform compatibility
• Memory hierarchy (global, shared, local) is critical for performance
• Optimization requires understanding of hardware architecture
• Profiling is essential for identifying bottlenecks

Whether you're accelerating scientific simulations, training neural networks, or processing large datasets, GPU programming is an essential skill for modern high-performance computing. Start with simple kernels, profile extensively, and gradually apply optimization techniques to achieve peak performance.