GPU programming has evolved from graphics rendering to become a cornerstone of enterprise computing, enabling massive parallel processing capabilities for complex computational workloads. This comprehensive guide explores advanced GPU programming techniques using CUDA and OpenCL, focusing on enterprise-grade applications and production deployment strategies.

GPU Architecture Understanding

Modern GPUs are designed for massive parallelism, featuring thousands of cores organized in streaming multiprocessors (SMs) or compute units (CUs). Understanding this architecture is crucial for effective GPU programming.

CUDA Architecture Fundamentals

CUDA organizes execution into a hierarchy of threads, blocks, and grids:

// CUDA kernel launch configuration
__global__ void matrix_multiply_kernel(float *A, float *B, float *C, 
                                     int N, int M, int K) {
    int row = blockIdx.y * blockDim.y + threadIdx.y;
    int col = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (row < N && col < M) {
        float sum = 0.0f;
        for (int k = 0; k < K; k++) {
            sum += A[row * K + k] * B[k * M + col];
        }
        C[row * M + col] = sum;
    }
}

// Host code for kernel launch
void launch_matrix_multiply(float *h_A, float *h_B, float *h_C,
                           int N, int M, int K) {
    size_t size_A = N * K * sizeof(float);
    size_t size_B = K * M * sizeof(float);
    size_t size_C = N * M * sizeof(float);
    
    float *d_A, *d_B, *d_C;
    
    // Allocate device memory
    cudaMalloc(&d_A, size_A);
    cudaMalloc(&d_B, size_B);
    cudaMalloc(&d_C, size_C);
    
    // Copy data to device
    cudaMemcpy(d_A, h_A, size_A, cudaMemcpyHostToDevice);
    cudaMemcpy(d_B, h_B, size_B, cudaMemcpyHostToDevice);
    
    // Configure launch parameters
    dim3 blockSize(16, 16);
    dim3 gridSize((M + blockSize.x - 1) / blockSize.x,
                  (N + blockSize.y - 1) / blockSize.y);
    
    // Launch kernel
    matrix_multiply_kernel<<<gridSize, blockSize>>>(d_A, d_B, d_C, N, M, K);
    
    // Check for kernel launch errors
    cudaError_t err = cudaGetLastError();
    if (err != cudaSuccess) {
        fprintf(stderr, "CUDA kernel launch failed: %s\n", 
                cudaGetErrorString(err));
        return;
    }
    
    // Copy result back to host
    cudaMemcpy(h_C, d_C, size_C, cudaMemcpyDeviceToHost);
    
    // Cleanup
    cudaFree(d_A);
    cudaFree(d_B);
    cudaFree(d_C);
}

OpenCL Architecture and Cross-Platform Development

OpenCL provides a standardized approach to heterogeneous computing across different vendors and device types:

#include <CL/cl.h>

typedef struct {
    cl_context context;
    cl_command_queue queue;
    cl_program program;
    cl_kernel kernel;
    cl_device_id device;
    cl_platform_id platform;
} opencl_context_t;

// OpenCL initialization and setup
int initialize_opencl(opencl_context_t *ctx, const char *kernel_source) {
    cl_int err;
    cl_uint num_platforms, num_devices;
    
    // Get platform
    err = clGetPlatformIDs(1, &ctx->platform, &num_platforms);
    if (err != CL_SUCCESS) return -1;
    
    // Get device
    err = clGetDeviceIDs(ctx->platform, CL_DEVICE_TYPE_GPU, 
                        1, &ctx->device, &num_devices);
    if (err != CL_SUCCESS) {
        // Fallback to CPU if no GPU available
        err = clGetDeviceIDs(ctx->platform, CL_DEVICE_TYPE_CPU,
                            1, &ctx->device, &num_devices);
        if (err != CL_SUCCESS) return -2;
    }
    
    // Create context
    ctx->context = clCreateContext(NULL, 1, &ctx->device, 
                                  NULL, NULL, &err);
    if (err != CL_SUCCESS) return -3;
    
    // Create command queue
    ctx->queue = clCreateCommandQueue(ctx->context, ctx->device, 
                                     CL_QUEUE_PROFILING_ENABLE, &err);
    if (err != CL_SUCCESS) return -4;
    
    // Create program from source
    ctx->program = clCreateProgramWithSource(ctx->context, 1,
                                            &kernel_source, NULL, &err);
    if (err != CL_SUCCESS) return -5;
    
    // Build program
    err = clBuildProgram(ctx->program, 1, &ctx->device, 
                        "-cl-fast-relaxed-math", NULL, NULL);
    if (err != CL_SUCCESS) {
        // Get build log for debugging
        size_t log_size;
        clGetProgramBuildInfo(ctx->program, ctx->device, 
                             CL_PROGRAM_BUILD_LOG, 0, NULL, &log_size);
        char *log = malloc(log_size);
        clGetProgramBuildInfo(ctx->program, ctx->device,
                             CL_PROGRAM_BUILD_LOG, log_size, log, NULL);
        fprintf(stderr, "OpenCL build error:\n%s\n", log);
        free(log);
        return -6;
    }
    
    return 0;
}

// OpenCL kernel for parallel reduction
const char *reduction_kernel_source = R"(
__kernel void parallel_reduction(__global float *input,
                                __global float *output,
                                __local float *local_mem,
                                int n) {
    int global_id = get_global_id(0);
    int local_id = get_local_id(0);
    int local_size = get_local_size(0);
    int group_id = get_group_id(0);
    
    // Load data into local memory
    if (global_id < n) {
        local_mem[local_id] = input[global_id];
    } else {
        local_mem[local_id] = 0.0f;
    }
    
    barrier(CLK_LOCAL_MEM_FENCE);
    
    // Perform reduction in local memory
    for (int stride = local_size / 2; stride > 0; stride >>= 1) {
        if (local_id < stride) {
            local_mem[local_id] += local_mem[local_id + stride];
        }
        barrier(CLK_LOCAL_MEM_FENCE);
    }
    
    // Write result to global memory
    if (local_id == 0) {
        output[group_id] = local_mem[0];
    }
}
)";

Memory Management and Optimization

Effective memory management is critical for GPU performance, requiring understanding of different memory types and access patterns.

CUDA Memory Hierarchy Optimization

// Optimized matrix multiplication using shared memory
__global__ void optimized_matrix_multiply(float *A, float *B, float *C,
                                        int N, int M, int K) {
    __shared__ float tile_A[16][16];
    __shared__ float tile_B[16][16];
    
    int bx = blockIdx.x, by = blockIdx.y;
    int tx = threadIdx.x, ty = threadIdx.y;
    
    int row = by * 16 + ty;
    int col = bx * 16 + tx;
    
    float sum = 0.0f;
    
    // Process tiles
    for (int tile = 0; tile < (K + 15) / 16; tile++) {
        // Load tile into shared memory
        if (row < N && tile * 16 + tx < K) {
            tile_A[ty][tx] = A[row * K + tile * 16 + tx];
        } else {
            tile_A[ty][tx] = 0.0f;
        }
        
        if (col < M && tile * 16 + ty < K) {
            tile_B[ty][tx] = B[(tile * 16 + ty) * M + col];
        } else {
            tile_B[ty][tx] = 0.0f;
        }
        
        __syncthreads();
        
        // Compute partial result
        for (int k = 0; k < 16; k++) {
            sum += tile_A[ty][k] * tile_B[k][tx];
        }
        
        __syncthreads();
    }
    
    // Write result
    if (row < N && col < M) {
        C[row * M + col] = sum;
    }
}

// Memory coalescing optimization
__global__ void coalesced_transpose(float *input, float *output,
                                   int width, int height) {
    __shared__ float tile[32][32];
    
    int x = blockIdx.x * 32 + threadIdx.x;
    int y = blockIdx.y * 32 + threadIdx.y;
    
    // Coalesced read from global memory
    if (x < width && y < height) {
        tile[threadIdx.y][threadIdx.x] = input[y * width + x];
    }
    
    __syncthreads();
    
    // Calculate transposed coordinates
    x = blockIdx.y * 32 + threadIdx.x;
    y = blockIdx.x * 32 + threadIdx.y;
    
    // Coalesced write to global memory
    if (x < height && y < width) {
        output[y * height + x] = tile[threadIdx.x][threadIdx.y];
    }
}

// Unified memory management for simplified programming
void unified_memory_example(int n) {
    float *data;
    size_t size = n * sizeof(float);
    
    // Allocate unified memory
    cudaMallocManaged(&data, size);
    
    // Initialize on CPU
    for (int i = 0; i < n; i++) {
        data[i] = i * 0.5f;
    }
    
    // Launch kernel - data automatically migrated to GPU
    int threadsPerBlock = 256;
    int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
    
    vector_scale<<<blocksPerGrid, threadsPerBlock>>>(data, 2.0f, n);
    
    // Synchronize before CPU access
    cudaDeviceSynchronize();
    
    // Access results on CPU - data automatically migrated back
    float sum = 0.0f;
    for (int i = 0; i < n; i++) {
        sum += data[i];
    }
    
    cudaFree(data);
}

Advanced Memory Patterns

// Pinned memory for faster transfers
void* allocate_pinned_memory(size_t size) {
    void *ptr;
    cudaError_t err = cudaMallocHost(&ptr, size);
    if (err != cudaSuccess) {
        return NULL;
    }
    return ptr;
}

// Asynchronous memory transfers with streams
void async_memory_transfer_example(float *h_data, float *d_data,
                                  size_t size, int num_streams) {
    cudaStream_t *streams = malloc(num_streams * sizeof(cudaStream_t));
    
    // Create streams
    for (int i = 0; i < num_streams; i++) {
        cudaStreamCreate(&streams[i]);
    }
    
    size_t chunk_size = size / num_streams;
    
    // Launch asynchronous transfers
    for (int i = 0; i < num_streams; i++) {
        size_t offset = i * chunk_size;
        size_t current_chunk = (i == num_streams - 1) ? 
                              size - offset : chunk_size;
        
        cudaMemcpyAsync(&d_data[offset], &h_data[offset],
                       current_chunk * sizeof(float),
                       cudaMemcpyHostToDevice, streams[i]);
    }
    
    // Synchronize all streams
    for (int i = 0; i < num_streams; i++) {
        cudaStreamSynchronize(streams[i]);
        cudaStreamDestroy(streams[i]);
    }
    
    free(streams);
}

// Memory pool for frequent allocations
typedef struct {
    void **free_blocks;
    size_t *block_sizes;
    int num_blocks;
    int capacity;
    size_t total_memory;
    pthread_mutex_t mutex;
} gpu_memory_pool_t;

void* gpu_pool_alloc(gpu_memory_pool_t *pool, size_t size) {
    pthread_mutex_lock(&pool->mutex);
    
    // Find suitable block
    for (int i = 0; i < pool->num_blocks; i++) {
        if (pool->block_sizes[i] >= size) {
            void *ptr = pool->free_blocks[i];
            
            // Remove from free list
            pool->num_blocks--;
            pool->free_blocks[i] = pool->free_blocks[pool->num_blocks];
            pool->block_sizes[i] = pool->block_sizes[pool->num_blocks];
            
            pthread_mutex_unlock(&pool->mutex);
            return ptr;
        }
    }
    
    pthread_mutex_unlock(&pool->mutex);
    
    // No suitable block found, allocate new one
    void *ptr;
    cudaError_t err = cudaMalloc(&ptr, size);
    return (err == cudaSuccess) ? ptr : NULL;
}

Parallel Algorithm Patterns

Effective GPU programming requires understanding fundamental parallel patterns and their implementations.

Reduction Patterns

// Optimized parallel reduction with warp-level primitives
__global__ void warp_reduction_kernel(float *input, float *output, int n) {
    __shared__ float sdata[256];
    
    int tid = threadIdx.x;
    int bid = blockIdx.x;
    int i = bid * blockDim.x + tid;
    
    // Load data
    sdata[tid] = (i < n) ? input[i] : 0.0f;
    __syncthreads();
    
    // Reduction in shared memory
    for (int s = blockDim.x / 2; s > 32; s >>= 1) {
        if (tid < s) {
            sdata[tid] += sdata[tid + s];
        }
        __syncthreads();
    }
    
    // Warp-level reduction
    if (tid < 32) {
        volatile float *vmem = sdata;
        vmem[tid] += vmem[tid + 32];
        vmem[tid] += vmem[tid + 16];
        vmem[tid] += vmem[tid + 8];
        vmem[tid] += vmem[tid + 4];
        vmem[tid] += vmem[tid + 2];
        vmem[tid] += vmem[tid + 1];
    }
    
    if (tid == 0) {
        output[bid] = sdata[0];
    }
}

// Scan (prefix sum) implementation
__global__ void prefix_sum_kernel(float *input, float *output, int n) {
    __shared__ float temp[512];
    
    int tid = threadIdx.x;
    int offset = 1;
    
    // Load data
    if (2 * tid < n) temp[2 * tid] = input[2 * tid];
    else temp[2 * tid] = 0;
    
    if (2 * tid + 1 < n) temp[2 * tid + 1] = input[2 * tid + 1];
    else temp[2 * tid + 1] = 0;
    
    // Up-sweep (build sum tree)
    for (int d = n >> 1; d > 0; d >>= 1) {
        __syncthreads();
        if (tid < d) {
            int ai = offset * (2 * tid + 1) - 1;
            int bi = offset * (2 * tid + 2) - 1;
            temp[bi] += temp[ai];
        }
        offset *= 2;
    }
    
    // Clear last element
    if (tid == 0) temp[n - 1] = 0;
    
    // Down-sweep
    for (int d = 1; d < n; d *= 2) {
        offset >>= 1;
        __syncthreads();
        if (tid < d) {
            int ai = offset * (2 * tid + 1) - 1;
            int bi = offset * (2 * tid + 2) - 1;
            float t = temp[ai];
            temp[ai] = temp[bi];
            temp[bi] += t;
        }
    }
    
    __syncthreads();
    
    // Write results
    if (2 * tid < n) output[2 * tid] = temp[2 * tid];
    if (2 * tid + 1 < n) output[2 * tid + 1] = temp[2 * tid + 1];
}

Sorting and Searching Algorithms

// Bitonic sort for power-of-2 sized arrays
__global__ void bitonic_sort_step(float *data, int j, int k, int n) {
    int i = threadIdx.x + blockDim.x * blockIdx.x;
    int ixj = i ^ j;
    
    if (ixj > i && i < n && ixj < n) {
        bool ascending = ((i & k) == 0);
        if ((data[i] > data[ixj]) == ascending) {
            // Swap elements
            float temp = data[i];
            data[i] = data[ixj];
            data[ixj] = temp;
        }
    }
}

void bitonic_sort_gpu(float *data, int n) {
    float *d_data;
    cudaMalloc(&d_data, n * sizeof(float));
    cudaMemcpy(d_data, data, n * sizeof(float), cudaMemcpyHostToDevice);
    
    int threadsPerBlock = 256;
    int blocksPerGrid = (n + threadsPerBlock - 1) / threadsPerBlock;
    
    for (int k = 2; k <= n; k *= 2) {
        for (int j = k / 2; j > 0; j /= 2) {
            bitonic_sort_step<<<blocksPerGrid, threadsPerBlock>>>(
                d_data, j, k, n);
            cudaDeviceSynchronize();
        }
    }
    
    cudaMemcpy(data, d_data, n * sizeof(float), cudaMemcpyDeviceToHost);
    cudaFree(d_data);
}

// Parallel binary search
__global__ void binary_search_kernel(float *sorted_array, float *targets,
                                    int *results, int array_size, 
                                    int num_targets) {
    int tid = blockIdx.x * blockDim.x + threadIdx.x;
    
    if (tid < num_targets) {
        float target = targets[tid];
        int left = 0, right = array_size - 1;
        int result = -1;
        
        while (left <= right) {
            int mid = (left + right) / 2;
            float mid_val = sorted_array[mid];
            
            if (mid_val == target) {
                result = mid;
                break;
            } else if (mid_val < target) {
                left = mid + 1;
            } else {
                right = mid - 1;
            }
        }
        
        results[tid] = result;
    }
}

Multi-GPU Programming and Scaling

Enterprise applications often require multi-GPU coordination for maximum performance:

// Multi-GPU context management
typedef struct {
    int num_devices;
    int *device_ids;
    cudaStream_t *streams;
    void **device_ptrs;
    size_t *device_memory_sizes;
} multi_gpu_context_t;

int initialize_multi_gpu(multi_gpu_context_t *ctx) {
    cudaGetDeviceCount(&ctx->num_devices);
    
    if (ctx->num_devices < 2) {
        return -1; // Not enough GPUs
    }
    
    ctx->device_ids = malloc(ctx->num_devices * sizeof(int));
    ctx->streams = malloc(ctx->num_devices * sizeof(cudaStream_t));
    ctx->device_ptrs = malloc(ctx->num_devices * sizeof(void*));
    ctx->device_memory_sizes = malloc(ctx->num_devices * sizeof(size_t));
    
    // Initialize each device
    for (int i = 0; i < ctx->num_devices; i++) {
        ctx->device_ids[i] = i;
        cudaSetDevice(i);
        cudaStreamCreate(&ctx->streams[i]);
        
        // Enable peer access for direct GPU-to-GPU transfers
        for (int j = 0; j < ctx->num_devices; j++) {
            if (i != j) {
                int can_access;
                cudaDeviceCanAccessPeer(&can_access, i, j);
                if (can_access) {
                    cudaDeviceEnablePeerAccess(j, 0);
                }
            }
        }
    }
    
    return 0;
}

// Data distribution across multiple GPUs
void distribute_data_multi_gpu(multi_gpu_context_t *ctx, 
                              float *host_data, size_t total_size) {
    size_t chunk_size = total_size / ctx->num_devices;
    
    for (int i = 0; i < ctx->num_devices; i++) {
        cudaSetDevice(i);
        
        size_t current_chunk = (i == ctx->num_devices - 1) ?
                              total_size - i * chunk_size : chunk_size;
        
        ctx->device_memory_sizes[i] = current_chunk * sizeof(float);
        cudaMalloc(&ctx->device_ptrs[i], ctx->device_memory_sizes[i]);
        
        // Asynchronous transfer to each GPU
        cudaMemcpyAsync(ctx->device_ptrs[i], 
                       &host_data[i * chunk_size],
                       ctx->device_memory_sizes[i],
                       cudaMemcpyHostToDevice,
                       ctx->streams[i]);
    }
    
    // Synchronize all transfers
    for (int i = 0; i < ctx->num_devices; i++) {
        cudaSetDevice(i);
        cudaStreamSynchronize(ctx->streams[i]);
    }
}

// Multi-GPU reduction operation
float multi_gpu_reduction(multi_gpu_context_t *ctx) {
    float *partial_results = malloc(ctx->num_devices * sizeof(float));
    
    // Launch reduction on each GPU
    for (int i = 0; i < ctx->num_devices; i++) {
        cudaSetDevice(i);
        
        float *d_result;
        cudaMalloc(&d_result, sizeof(float));
        
        // Launch single-GPU reduction kernel
        int num_elements = ctx->device_memory_sizes[i] / sizeof(float);
        launch_reduction_kernel(ctx->device_ptrs[i], d_result, num_elements);
        
        // Copy result back
        cudaMemcpy(&partial_results[i], d_result, sizeof(float),
                  cudaMemcpyDeviceToHost);
        
        cudaFree(d_result);
    }
    
    // Final reduction on CPU
    float final_result = 0.0f;
    for (int i = 0; i < ctx->num_devices; i++) {
        final_result += partial_results[i];
    }
    
    free(partial_results);
    return final_result;
}

Performance Profiling and Optimization

Systematic performance analysis is essential for enterprise GPU applications:

// CUDA profiling integration
void profile_kernel_execution(void (*kernel_func)(void), 
                             const char *kernel_name) {
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    // Start timing
    cudaEventRecord(start);
    
    // Execute kernel
    kernel_func();
    
    // Stop timing
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float milliseconds = 0;
    cudaEventElapsedTime(&milliseconds, start, stop);
    
    printf("Kernel %s execution time: %.3f ms\n", 
           kernel_name, milliseconds);
    
    // Get additional metrics
    size_t free_mem, total_mem;
    cudaMemGetInfo(&free_mem, &total_mem);
    printf("GPU memory usage: %.1f%% (%.1f MB / %.1f MB)\n",
           100.0 * (total_mem - free_mem) / total_mem,
           (total_mem - free_mem) / 1048576.0,
           total_mem / 1048576.0);
    
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}

// Memory bandwidth measurement
void measure_memory_bandwidth(size_t data_size) {
    float *h_data = malloc(data_size);
    float *d_data;
    
    cudaMalloc(&d_data, data_size);
    
    // Initialize data
    for (size_t i = 0; i < data_size / sizeof(float); i++) {
        h_data[i] = i * 0.5f;
    }
    
    cudaEvent_t start, stop;
    cudaEventCreate(&start);
    cudaEventCreate(&stop);
    
    // Measure host to device transfer
    cudaEventRecord(start);
    cudaMemcpy(d_data, h_data, data_size, cudaMemcpyHostToDevice);
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float h2d_time;
    cudaEventElapsedTime(&h2d_time, start, stop);
    
    // Measure device to host transfer
    cudaEventRecord(start);
    cudaMemcpy(h_data, d_data, data_size, cudaMemcpyDeviceToHost);
    cudaEventRecord(stop);
    cudaEventSynchronize(stop);
    
    float d2h_time;
    cudaEventElapsedTime(&d2h_time, start, stop);
    
    double data_gb = data_size / 1073741824.0;
    printf("H2D bandwidth: %.2f GB/s\n", data_gb / (h2d_time / 1000.0));
    printf("D2H bandwidth: %.2f GB/s\n", data_gb / (d2h_time / 1000.0));
    
    free(h_data);
    cudaFree(d_data);
    cudaEventDestroy(start);
    cudaEventDestroy(stop);
}

// Occupancy optimization
__global__ void occupancy_test_kernel(float *data, int n) {
    int idx = blockIdx.x * blockDim.x + threadIdx.x;
    if (idx < n) {
        // Simulate computation
        float value = data[idx];
        for (int i = 0; i < 100; i++) {
            value = sinf(value) + cosf(value);
        }
        data[idx] = value;
    }
}

void optimize_occupancy(float *d_data, int n) {
    int min_grid_size, block_size;
    
    // Calculate optimal block size for maximum occupancy
    cudaOccupancyMaxPotentialBlockSize(&min_grid_size, &block_size,
                                      occupancy_test_kernel, 0, 0);
    
    int grid_size = (n + block_size - 1) / block_size;
    
    printf("Optimal block size: %d\n", block_size);
    printf("Grid size: %d\n", grid_size);
    
    // Calculate theoretical occupancy
    int max_active_blocks;
    cudaOccupancyMaxActiveBlocksPerMultiprocessor(&max_active_blocks,
                                                 occupancy_test_kernel,
                                                 block_size, 0);
    
    int device;
    cudaGetDevice(&device);
    
    cudaDeviceProp prop;
    cudaGetDeviceProperties(&prop, device);
    
    float occupancy = (max_active_blocks * block_size / 
                      (float)prop.maxThreadsPerMultiProcessor) * 100;
    
    printf("Theoretical occupancy: %.1f%%\n", occupancy);
    
    // Launch with optimized parameters
    occupancy_test_kernel<<<grid_size, block_size>>>(d_data, n);
}

Enterprise Integration and Deployment

Production GPU applications require robust integration patterns and deployment strategies:

// GPU resource management system
typedef struct {
    int total_gpus;
    int *gpu_utilization;  // Percentage utilization per GPU
    int *gpu_memory_usage; // Memory usage per GPU
    pthread_mutex_t *gpu_locks;
    int *gpu_queue_lengths;
    time_t *last_activity;
} gpu_resource_manager_t;

gpu_resource_manager_t* initialize_gpu_manager() {
    gpu_resource_manager_t *manager = malloc(sizeof(gpu_resource_manager_t));
    
    cudaGetDeviceCount(&manager->total_gpus);
    
    manager->gpu_utilization = calloc(manager->total_gpus, sizeof(int));
    manager->gpu_memory_usage = calloc(manager->total_gpus, sizeof(int));
    manager->gpu_locks = malloc(manager->total_gpus * sizeof(pthread_mutex_t));
    manager->gpu_queue_lengths = calloc(manager->total_gpus, sizeof(int));
    manager->last_activity = calloc(manager->total_gpus, sizeof(time_t));
    
    for (int i = 0; i < manager->total_gpus; i++) {
        pthread_mutex_init(&manager->gpu_locks[i], NULL);
    }
    
    return manager;
}

int acquire_gpu_resource(gpu_resource_manager_t *manager, 
                        int required_memory_mb) {
    int best_gpu = -1;
    int lowest_utilization = 100;
    
    for (int i = 0; i < manager->total_gpus; i++) {
        cudaSetDevice(i);
        
        size_t free_mem, total_mem;
        cudaMemGetInfo(&free_mem, &total_mem);
        
        int available_mb = free_mem / 1048576;
        
        if (available_mb >= required_memory_mb &&
            manager->gpu_utilization[i] < lowest_utilization) {
            best_gpu = i;
            lowest_utilization = manager->gpu_utilization[i];
        }
    }
    
    if (best_gpu >= 0) {
        pthread_mutex_lock(&manager->gpu_locks[best_gpu]);
        manager->gpu_utilization[best_gpu]++;
        manager->last_activity[best_gpu] = time(NULL);
        pthread_mutex_unlock(&manager->gpu_locks[best_gpu]);
    }
    
    return best_gpu;
}

// Error handling and recovery
typedef enum {
    GPU_ERROR_NONE = 0,
    GPU_ERROR_OUT_OF_MEMORY,
    GPU_ERROR_LAUNCH_FAILED,
    GPU_ERROR_INVALID_VALUE,
    GPU_ERROR_DEVICE_LOST,
    GPU_ERROR_TIMEOUT
} gpu_error_type_t;

typedef struct {
    gpu_error_type_t error_type;
    int device_id;
    char error_message[512];
    time_t timestamp;
    int retry_count;
} gpu_error_context_t;

int handle_gpu_error(gpu_error_context_t *error_ctx, cudaError_t cuda_error) {
    error_ctx->timestamp = time(NULL);
    error_ctx->retry_count++;
    
    switch (cuda_error) {
        case cudaErrorMemoryAllocation:
            error_ctx->error_type = GPU_ERROR_OUT_OF_MEMORY;
            strcpy(error_ctx->error_message, "GPU out of memory");
            
            // Attempt garbage collection
            cudaDeviceSynchronize();
            
            // Try to free unused allocations
            if (error_ctx->retry_count < 3) {
                return 1; // Retry allowed
            }
            break;
            
        case cudaErrorLaunchFailure:
            error_ctx->error_type = GPU_ERROR_LAUNCH_FAILED;
            strcpy(error_ctx->error_message, "Kernel launch failed");
            
            // Reset device and retry
            cudaDeviceReset();
            if (error_ctx->retry_count < 2) {
                return 1; // Retry allowed
            }
            break;
            
        case cudaErrorDevicesUnavailable:
            error_ctx->error_type = GPU_ERROR_DEVICE_LOST;
            strcpy(error_ctx->error_message, "GPU device unavailable");
            
            // Switch to CPU fallback
            return -1; // No retry, use fallback
            
        default:
            error_ctx->error_type = GPU_ERROR_INVALID_VALUE;
            snprintf(error_ctx->error_message, sizeof(error_ctx->error_message),
                    "CUDA error: %s", cudaGetErrorString(cuda_error));
            return -1; // No retry
    }
    
    return 0; // No retry
}

// Production monitoring and logging
void log_gpu_metrics(gpu_resource_manager_t *manager) {
    for (int i = 0; i < manager->total_gpus; i++) {
        cudaSetDevice(i);
        
        size_t free_mem, total_mem;
        cudaMemGetInfo(&free_mem, &total_mem);
        
        cudaDeviceProp prop;
        cudaGetDeviceProperties(&prop, i);
        
        printf("GPU %d (%s):\n", i, prop.name);
        printf("  Memory: %.1f MB / %.1f MB (%.1f%% used)\n",
               (total_mem - free_mem) / 1048576.0,
               total_mem / 1048576.0,
               100.0 * (total_mem - free_mem) / total_mem);
        printf("  Utilization: %d%%\n", manager->gpu_utilization[i]);
        printf("  Queue length: %d\n", manager->gpu_queue_lengths[i]);
        
        // Temperature monitoring (if available)
        int temp;
        if (cudaDeviceGetAttribute(&temp, cudaDevAttrGpuOverlapCanMap, i) 
            == cudaSuccess) {
            printf("  Temperature: %d°C\n", temp);
        }
    }
}

Conclusion

GPU programming for enterprise workloads requires mastery of parallel computing principles, memory optimization techniques, and robust production deployment strategies. The techniques presented in this guide provide a comprehensive foundation for developing high-performance, scalable GPU applications that can handle enterprise-scale computational demands.

Key considerations for enterprise GPU development include proper resource management, comprehensive error handling, multi-GPU coordination, and systematic performance optimization. By leveraging both CUDA and OpenCL technologies appropriately, developers can create portable, efficient solutions that maximize the computational potential of modern GPU hardware.

The examples and patterns demonstrated here form the basis for building sophisticated parallel applications that can scale from single-GPU workstations to large-scale data centers with hundreds of GPUs, enabling organizations to tackle previously intractable computational challenges with confidence and efficiency.