GPU Computing for Quantitative Finance: CUDA vs OpenCL vs Vulkan Compute

1Platform Comparison:
2
3CUDA (NVIDIA only):
4+ Best performance on NVIDIA GPUs
5+ Mature ecosystem (cuBLAS, cuFFT, etc.)
6+ Excellent tooling (Nsight, nvprof)
7+ Most documentation and examples
8- NVIDIA hardware only
9- Proprietary
10
11OpenCL (Cross-platform):
12+ Works on NVIDIA, AMD, Intel
13+ Open standard
14+ CPU fallback support
15+ Good for portability
16- Lower performance than CUDA
17- Less mature ecosystem
18- Inconsistent vendor support
19
20Vulkan Compute (Modern):
21+ Cross-platform (GPU + mobile)
22+ Modern API design
23+ Low overhead
24+ Good for mixed graphics/compute
25- Verbose API
26- Less mature for compute
27- Smaller ecosystem
28
29Choice for Finance:
30- CUDA: Best for pure performance (NVIDIA hardware)
31- OpenCL: Best for portability
32- Vulkan: Best for mixed workloads
33

1#include <cuda_runtime.h>
2#include <curand_kernel.h>
3#include <iostream>
4#include <chrono>
5
6// Black-Scholes formula for European call option
7__device__ float black_scholes_call(
8    float S,      // Stock price
9    float K,      // Strike price
10    float r,      // Risk-free rate
11    float sigma,  // Volatility
12    float T       // Time to maturity
13) {
14    float d1 = (logf(S / K) + (r + 0.5f * sigma * sigma) * T) / (sigma * sqrtf(T));
15    float d2 = d1 - sigma * sqrtf(T);
16    
17    // Normal CDF approximation
18    auto normcdf = [](float x) {
19        return 0.5f * (1.0f + erff(x / sqrtf(2.0f)));
20    };
21    
22    float call_price = S * normcdf(d1) - K * expf(-r * T) * normcdf(d2);
23    return call_price;
24}
25
26// Monte Carlo kernel for European call option
27__global__ void monte_carlo_option_pricing(
28    float* prices,           // Output: option prices
29    float S0,                // Initial stock price
30    float K,                 // Strike price
31    float r,                 // Risk-free rate
32    float sigma,             // Volatility
33    float T,                 // Time to maturity
34    int num_paths,           // Number of Monte Carlo paths
35    unsigned long long seed
36) {
37    int idx = blockIdx.x * blockDim.x + threadIdx.x;
38    
39    if (idx >= num_paths) return;
40    
41    // Initialize RNG
42    curandState state;
43    curand_init(seed, idx, 0, &state);
44    
45    // Generate random path
46    float z = curand_normal(&state);
47    
48    // Geometric Brownian Motion
49    float ST = S0 * expf((r - 0.5f * sigma * sigma) * T + sigma * sqrtf(T) * z);
50    
51    // Payoff for call option
52    float payoff = fmaxf(ST - K, 0.0f);
53    
54    // Discounted payoff
55    prices[idx] = expf(-r * T) * payoff;
56}
57
58// Reduction kernel to sum option prices
59__global__ void reduce_sum(
60    const float* input,
61    float* output,
62    int n
63) {
64    extern __shared__ float sdata[];
65    
66    unsigned int tid = threadIdx.x;
67    unsigned int idx = blockIdx.x * blockDim.x + threadIdx.x;
68    
69    // Load data into shared memory
70    sdata[tid] = (idx < n) ? input[idx] : 0.0f;
71    __syncthreads();
72    
73    // Reduction in shared memory
74    for (unsigned int s = blockDim.x / 2; s > 0; s >>= 1) {
75        if (tid < s) {
76            sdata[tid] += sdata[tid + s];
77        }
78        __syncthreads();
79    }
80    
81    // Write result for this block
82    if (tid == 0) {
83        output[blockIdx.x] = sdata[tid];
84    }
85}
86
87class MonteCarloGPU {
88private:
89    float* d_prices;
90    float* d_reduction;
91    int num_paths;
92    
93public:
94    MonteCarloGPU(int paths) : num_paths(paths) {
95        // Allocate device memory
96        cudaMalloc(&d_prices, num_paths * sizeof(float));
97        cudaMalloc(&d_reduction, (num_paths / 256 + 1) * sizeof(float));
98    }
99    
100    ~MonteCarloGPU() {
101        cudaFree(d_prices);
102        cudaFree(d_reduction);
103    }
104    
105    float price_option(
106        float S0,
107        float K,
108        float r,
109        float sigma,
110        float T
111    ) {
112        // Launch configuration
113        int threadsPerBlock = 256;
114        int blocksPerGrid = (num_paths + threadsPerBlock - 1) / threadsPerBlock;
115        
116        // Generate random seed
117        unsigned long long seed = std::chrono::system_clock::now().time_since_epoch().count();
118        
119        // Launch Monte Carlo kernel
120        monte_carlo_option_pricing<<<blocksPerGrid, threadsPerBlock>>>(
121            d_prices, S0, K, r, sigma, T, num_paths, seed
122        );
123        
124        // Reduce to get average
125        int reductionBlocks = (num_paths / 256 + 1);
126        reduce_sum<<<reductionBlocks, 256, 256 * sizeof(float)>>>(
127            d_prices, d_reduction, num_paths
128        );
129        
130        // Copy final results back
131        std::vector<float> h_reduction(reductionBlocks);
132        cudaMemcpy(h_reduction.data(), d_reduction, 
133                   reductionBlocks * sizeof(float), cudaMemcpyDeviceToHost);
134        
135        // Final sum on CPU
136        float total = 0.0f;
137        for (float val : h_reduction) {
138            total += val;
139        }
140        
141        return total / num_paths;
142    }
143    
144    // Benchmark
145    void benchmark() {
146        float S0 = 100.0f;
147        float K = 100.0f;
148        float r = 0.05f;
149        float sigma = 0.2f;
150        float T = 1.0f;
151        
152        // Warmup
153        price_option(S0, K, r, sigma, T);
154        
155        // Benchmark
156        auto start = std::chrono::high_resolution_clock::now();
157        
158        const int iterations = 100;
159        float price = 0.0f;
160        for (int i = 0; i < iterations; ++i) {
161            price = price_option(S0, K, r, sigma, T);
162        }
163        
164        auto end = std::chrono::high_resolution_clock::now();
165        auto duration = std::chrono::duration_cast<std::chrono::microseconds>(end - start).count();
166        
167        std::cout << "=== GPU Monte Carlo Benchmark ===\n";
168        std::cout << "Paths: " << num_paths << "\n";
169        std::cout << "Iterations: " << iterations << "\n";
170        std::cout << "Option price: " << price << "\n";
171        std::cout << "Time per iteration: " << duration / iterations << " μs\n";
172        std::cout << "Throughput: " << (num_paths * iterations * 1000000.0) / duration << " paths/sec\n";
173    }
174};
175
176// CPU reference implementation
177float monte_carlo_cpu(
178    float S0, float K, float r, float sigma, float T, int num_paths
179) {
180    std::random_device rd;
181    std::mt19937 gen(rd());
182    std::normal_distribution<float> dis(0.0f, 1.0f);
183    
184    float sum = 0.0f;
185    
186    for (int i = 0; i < num_paths; ++i) {
187        float z = dis(gen);
188        float ST = S0 * expf((r - 0.5f * sigma * sigma) * T + sigma * sqrtf(T) * z);
189        float payoff = std::max(ST - K, 0.0f);
190        sum += expf(-r * T) * payoff;
191    }
192    
193    return sum / num_paths;
194}
195
196int main() {
197    // GPU version
198    MonteCarloGPU gpu(1000000);  // 1M paths
199    gpu.benchmark();
200    
201    // CPU version for comparison
202    auto start = std::chrono::high_resolution_clock::now();
203    float cpu_price = monte_carlo_cpu(100.0f, 100.0f, 0.05f, 0.2f, 1.0f, 1000000);
204    auto end = std::chrono::high_resolution_clock::now();
205    auto cpu_time = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
206    
207    std::cout << "\n=== CPU Comparison ===\n";
208    std::cout << "CPU price: " << cpu_price << "\n";
209    std::cout << "CPU time: " << cpu_time << " ms\n";
210    
211    return 0;
212}
213
214// Compile:
215// nvcc -O3 -arch=sm_80 monte_carlo_gpu.cu -o monte_carlo_gpu
216

1// Antithetic variates: use both z and -z
2__global__ void monte_carlo_antithetic(
3    float* prices,
4    float S0, float K, float r, float sigma, float T,
5    int num_paths,
6    unsigned long long seed
7) {
8    int idx = blockIdx.x * blockDim.x + threadIdx.x;
9    int pair_idx = idx / 2;  // Process pairs
10    
11    if (pair_idx >= num_paths / 2) return;
12    
13    curandState state;
14    curand_init(seed, pair_idx, 0, &state);
15    
16    float z = curand_normal(&state);
17    
18    // Path 1: use z
19    float ST1 = S0 * expf((r - 0.5f * sigma * sigma) * T + sigma * sqrtf(T) * z);
20    float payoff1 = fmaxf(ST1 - K, 0.0f);
21    
22    // Path 2: use -z (antithetic)
23    float ST2 = S0 * expf((r - 0.5f * sigma * sigma) * T + sigma * sqrtf(T) * (-z));
24    float payoff2 = fmaxf(ST2 - K, 0.0f);
25    
26    // Average of both paths
27    float avg_payoff = (payoff1 + payoff2) / 2.0f;
28    
29    prices[pair_idx] = expf(-r * T) * avg_payoff;
30}
31
32// Control variates: use known analytical solution
33__global__ void monte_carlo_control_variates(
34    float* prices,
35    float S0, float K, float r, float sigma, float T,
36    int num_paths,
37    unsigned long long seed
38) {
39    int idx = blockIdx.x * blockDim.x + threadIdx.x;
40    
41    if (idx >= num_paths) return;
42    
43    curandState state;
44    curand_init(seed, idx, 0, &state);
45    
46    float z = curand_normal(&state);
47    
48    // Simulate path
49    float ST = S0 * expf((r - 0.5f * sigma * sigma) * T + sigma * sqrtf(T) * z);
50    float payoff = fmaxf(ST - K, 0.0f);
51    
52    // Control variate: use geometric average (has analytical solution)
53    float geometric_avg = S0;  // Simplified
54    float control_payoff = fmaxf(geometric_avg - K, 0.0f);
55    
56    // Known expected value of control
57    float control_expected = 10.0f;  // Placeholder - calculate analytically
58    
59    // Adjust using control variate
60    float beta = 0.9f;  // Optimal beta found empirically
61    float adjusted = payoff - beta * (control_payoff - control_expected);
62    
63    prices[idx] = expf(-r * T) * adjusted;
64}
65

1struct OptionGreeks {
2    float price;
3    float delta;
4    float gamma;
5    float vega;
6    float theta;
7    float rho;
8};
9
10__global__ void calculate_greeks_gpu(
11    OptionGreeks* greeks,
12    float S0, float K, float r, float sigma, float T,
13    int num_options
14) {
15    int idx = blockIdx.x * blockDim.x + threadIdx.x;
16    
17    if (idx >= num_options) return;
18    
19    // Adjust parameters for this option (example: different strikes)
20    float K_adjusted = K + idx * 1.0f;
21    
22    // Calculate d1, d2
23    float d1 = (logf(S0 / K_adjusted) + (r + 0.5f * sigma * sigma) * T) / (sigma * sqrtf(T));
24    float d2 = d1 - sigma * sqrtf(T);
25    
26    // Normal CDF and PDF
27    auto N = [](float x) {
28        return 0.5f * (1.0f + erff(x / sqrtf(2.0f)));
29    };
30    
31    auto n = [](float x) {
32        return expf(-0.5f * x * x) / sqrtf(2.0f * 3.14159265f);
33    };
34    
35    // Price
36    float price = S0 * N(d1) - K_adjusted * expf(-r * T) * N(d2);
37    
38    // Delta
39    float delta = N(d1);
40    
41    // Gamma
42    float gamma = n(d1) / (S0 * sigma * sqrtf(T));
43    
44    // Vega
45    float vega = S0 * n(d1) * sqrtf(T) / 100.0f;  // Divide by 100 for 1% vol change
46    
47    // Theta
48    float theta = -(S0 * n(d1) * sigma) / (2.0f * sqrtf(T)) - 
49                  r * K_adjusted * expf(-r * T) * N(d2);
50    theta /= 365.0f;  // Per day
51    
52    // Rho
53    float rho = K_adjusted * T * expf(-r * T) * N(d2) / 100.0f;  // Per 1% rate change
54    
55    // Store results
56    greeks[idx].price = price;
57    greeks[idx].delta = delta;
58    greeks[idx].gamma = gamma;
59    greeks[idx].vega = vega;
60    greeks[idx].theta = theta;
61    greeks[idx].rho = rho;
62}
63
64class GreeksCalculator {
65private:
66    OptionGreeks* d_greeks;
67    OptionGreeks* h_greeks;
68    int num_options;
69    
70public:
71    GreeksCalculator(int n) : num_options(n) {
72        cudaMalloc(&d_greeks, n * sizeof(OptionGreeks));
73        h_greeks = new OptionGreeks[n];
74    }
75    
76    ~GreeksCalculator() {
77        cudaFree(d_greeks);
78        delete[] h_greeks;
79    }
80    
81    void calculate(float S0, float K, float r, float sigma, float T) {
82        int threadsPerBlock = 256;
83        int blocksPerGrid = (num_options + threadsPerBlock - 1) / threadsPerBlock;
84        
85        calculate_greeks_gpu<<<blocksPerGrid, threadsPerBlock>>>(
86            d_greeks, S0, K, r, sigma, T, num_options
87        );
88        
89        // Copy results back
90        cudaMemcpy(h_greeks, d_greeks, num_options * sizeof(OptionGreeks),
91                   cudaMemcpyDeviceToHost);
92    }
93    
94    void print_results(int num_to_print = 5) {
95        std::cout << "=== Greeks Results (first " << num_to_print << ") ===\n";
96        for (int i = 0; i < std::min(num_to_print, num_options); ++i) {
97            std::cout << "Option " << i << ":\n";
98            std::cout << "  Price: " << h_greeks[i].price << "\n";
99            std::cout << "  Delta: " << h_greeks[i].delta << "\n";
100            std::cout << "  Gamma: " << h_greeks[i].gamma << "\n";
101            std::cout << "  Vega:  " << h_greeks[i].vega << "\n";
102            std::cout << "  Theta: " << h_greeks[i].theta << "\n";
103            std::cout << "  Rho:   " << h_greeks[i].rho << "\n\n";
104        }
105    }
106};
107

1class OptimizedMemoryTransfer {
2public:
3    // Pattern 1: Pinned Memory (faster transfers)
4    void use_pinned_memory() {
5        const int N = 100000000;  // 100M floats
6        
7        float* h_data;
8        float* d_data;
9        
10        // Allocate pinned (page-locked) host memory
11        cudaMallocHost(&h_data, N * sizeof(float));
12        cudaMalloc(&d_data, N * sizeof(float));
13        
14        // Initialize data
15        for (int i = 0; i < N; ++i) {
16            h_data[i] = i;
17        }
18        
19        // Transfer is 2-3x faster with pinned memory
20        auto start = std::chrono::high_resolution_clock::now();
21        cudaMemcpy(d_data, h_data, N * sizeof(float), cudaMemcpyHostToDevice);
22        auto end = std::chrono::high_resolution_clock::now();
23        
24        auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
25        float bandwidth = (N * sizeof(float)) / (duration / 1000.0) / 1e9;
26        
27        std::cout << "Pinned memory transfer: " << duration << " ms\n";
28        std::cout << "Bandwidth: " << bandwidth << " GB/s\n";
29        
30        cudaFree(d_data);
31        cudaFreeHost(h_data);
32    }
33    
34    // Pattern 2: Asynchronous Transfers (overlap with computation)
35    void async_transfer_compute() {
36        const int N = 10000000;
37        const int num_streams = 4;
38        
39        float* h_data[num_streams];
40        float* d_data[num_streams];
41        cudaStream_t streams[num_streams];
42        
43        int chunk_size = N / num_streams;
44        
45        // Setup
46        for (int i = 0; i < num_streams; ++i) {
47            cudaMallocHost(&h_data[i], chunk_size * sizeof(float));
48            cudaMalloc(&d_data[i], chunk_size * sizeof(float));
49            cudaStreamCreate(&streams[i]);
50            
51            // Initialize
52            for (int j = 0; j < chunk_size; ++j) {
53                h_data[i][j] = i * chunk_size + j;
54            }
55        }
56        
57        auto start = std::chrono::high_resolution_clock::now();
58        
59        // Pipeline: transfer and compute overlap
60        for (int i = 0; i < num_streams; ++i) {
61            // Async copy
62            cudaMemcpyAsync(d_data[i], h_data[i], chunk_size * sizeof(float),
63                           cudaMemcpyHostToDevice, streams[i]);
64            
65            // Launch kernel on this stream
66            int threadsPerBlock = 256;
67            int blocksPerGrid = (chunk_size + threadsPerBlock - 1) / threadsPerBlock;
68            
69            // Example kernel: square each element
70            // square_kernel<<<blocksPerGrid, threadsPerBlock, 0, streams[i]>>>(d_data[i], chunk_size);
71            
72            // Async copy results back
73            cudaMemcpyAsync(h_data[i], d_data[i], chunk_size * sizeof(float),
74                           cudaMemcpyDeviceToHost, streams[i]);
75        }
76        
77        // Wait for all streams
78        for (int i = 0; i < num_streams; ++i) {
79            cudaStreamSynchronize(streams[i]);
80        }
81        
82        auto end = std::chrono::high_resolution_clock::now();
83        auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
84        
85        std::cout << "Async pipelined transfer+compute: " << duration << " ms\n";
86        
87        // Cleanup
88        for (int i = 0; i < num_streams; ++i) {
89            cudaFree(d_data[i]);
90            cudaFreeHost(h_data[i]);
91            cudaStreamDestroy(streams[i]);
92        }
93    }
94    
95    // Pattern 3: Unified Memory (simplest, auto-managed)
96    void unified_memory() {
97        const int N = 10000000;
98        float* data;
99        
100        // Allocate unified memory (accessible from both CPU and GPU)
101        cudaMallocManaged(&data, N * sizeof(float));
102        
103        // Initialize on CPU
104        for (int i = 0; i < N; ++i) {
105            data[i] = i;
106        }
107        
108        // Use on GPU (data automatically transferred)
109        int threadsPerBlock = 256;
110        int blocksPerGrid = (N + threadsPerBlock - 1) / threadsPerBlock;
111        
112        // square_kernel<<<blocksPerGrid, threadsPerBlock>>>(data, N);
113        cudaDeviceSynchronize();
114        
115        // Access results on CPU (automatically transferred back)
116        float sum = 0.0f;
117        for (int i = 0; i < std::min(100, N); ++i) {
118            sum += data[i];
119        }
120        
121        std::cout << "Unified memory sum: " << sum << "\n";
122        
123        cudaFree(data);
124    }
125    
126    // Pattern 4: Zero-Copy Memory (avoid transfers for small data)
127    void zero_copy_memory() {
128        const int N = 1000;  // Small dataset
129        float* h_data;
130        float* d_data;
131        
132        // Allocate mapped memory
133        cudaHostAlloc(&h_data, N * sizeof(float), cudaHostAllocMapped);
134        cudaHostGetDevicePointer(&d_data, h_data, 0);
135        
136        // Initialize on CPU
137        for (int i = 0; i < N; ++i) {
138            h_data[i] = i;
139        }
140        
141        // GPU accesses host memory directly (no copy)
142        // Slower per access, but no transfer overhead
143        // Good for small/irregular access patterns
144        
145        // square_kernel<<<1, 256>>>(d_data, N);
146        cudaDeviceSynchronize();
147        
148        cudaFreeHost(h_data);
149    }
150};
151

1class MultiGPUBacktest {
2private:
3    int num_gpus;
4    float** d_prices;  // Array of device pointers
5    
6public:
7    MultiGPUBacktest() {
8        cudaGetDeviceCount(&num_gpus);
9        std::cout << "Found " << num_gpus << " GPUs\n";
10        
11        d_prices = new float*[num_gpus];
12    }
13    
14    ~MultiGPUBacktest() {
15        for (int i = 0; i < num_gpus; ++i) {
16            cudaSetDevice(i);
17            cudaFree(d_prices[i]);
18        }
19        delete[] d_prices;
20    }
21    
22    void run_multi_gpu_simulation(int total_simulations) {
23        int sims_per_gpu = total_simulations / num_gpus;
24        
25        // Allocate memory on each GPU
26        for (int i = 0; i < num_gpus; ++i) {
27            cudaSetDevice(i);
28            cudaMalloc(&d_prices[i], sims_per_gpu * sizeof(float));
29        }
30        
31        // Launch kernels on each GPU
32        cudaStream_t streams[num_gpus];
33        
34        auto start = std::chrono::high_resolution_clock::now();
35        
36        for (int i = 0; i < num_gpus; ++i) {
37            cudaSetDevice(i);
38            cudaStreamCreate(&streams[i]);
39            
40            int threadsPerBlock = 256;
41            int blocksPerGrid = (sims_per_gpu + threadsPerBlock - 1) / threadsPerBlock;
42            
43            unsigned long long seed = i * 1000000;
44            
45            monte_carlo_option_pricing<<<blocksPerGrid, threadsPerBlock, 0, streams[i]>>>(
46                d_prices[i], 100.0f, 100.0f, 0.05f, 0.2f, 1.0f, sims_per_gpu, seed
47            );
48        }
49        
50        // Wait for all GPUs
51        for (int i = 0; i < num_gpus; ++i) {
52            cudaSetDevice(i);
53            cudaStreamSynchronize(streams[i]);
54        }
55        
56        auto end = std::chrono::high_resolution_clock::now();
57        auto duration = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
58        
59        std::cout << "Multi-GPU simulation: " << total_simulations << " paths\n";
60        std::cout << "Time: " << duration << " ms\n";
61        std::cout << "Throughput: " << (total_simulations / (duration / 1000.0)) / 1e6 
62                  << " million paths/sec\n";
63        
64        // Cleanup
65        for (int i = 0; i < num_gpus; ++i) {
66            cudaSetDevice(i);
67            cudaStreamDestroy(streams[i]);
68        }
69    }
70};
71

1class GPUDecisionGuide {
2public:
3    // Anti-pattern 1: Small datasets
4    void small_dataset_bad() {
5        // DON'T: Transfer overhead > computation time
6        const int N = 100;
7        float* d_data;
8        cudaMalloc(&d_data, N * sizeof(float));
9        
10        // Transfer time: ~50 μs
11        // Computation time: ~1 μs
12        // Total: 51 μs
13        
14        // Better: Just do it on CPU (~0.1 μs)
15    }
16    
17    // Anti-pattern 2: Branchy code
18    void branchy_code_bad() {
19        // DON'T: Different threads take different paths
20        // Causes thread divergence and serialization
21        
22        /*
23        __global__ void bad_kernel(float* data, int n) {
24            int idx = threadIdx.x + blockIdx.x * blockDim.x;
25            if (idx >= n) return;
26            
27            // BAD: Lots of branching
28            if (data[idx] > 100) {
29                // Complex computation path A
30            } else if (data[idx] > 50) {
31                // Complex computation path B
32            } else if (data[idx] > 25) {
33                // Complex computation path C
34            } else {
35                // Complex computation path D
36            }
37        }
38        */
39    }
40    
41    // Anti-pattern 3: Sequential dependencies
42    void sequential_bad() {
43        // DON'T: Each step depends on previous
44        // Can't parallelize effectively
45        
46        /*
47        for (int i = 1; i < n; ++i) {
48            data[i] = data[i-1] + func(data[i]);
49        }
50        */
51        
52        // This is inherently sequential - CPU better
53    }
54    
55    // Good use case: Embarrassingly parallel
56    void embarrassingly_parallel_good() {
57        // DO: Independent computations
58        // Each thread processes independent data
59        
60        /*
61        __global__ void good_kernel(float* input, float* output, int n) {
62            int idx = threadIdx.x + blockIdx.x * blockDim.x;
63            if (idx >= n) return;
64            
65            // Each thread independent
66            output[idx] = expensive_computation(input[idx]);
67        }
68        */
69    }
70    
71    // Decision flowchart
72    bool should_use_gpu(
73        int dataset_size,
74        bool has_branching,
75        bool has_dependencies,
76        float computation_time_cpu_ms
77    ) {
78        // Rule 1: Dataset too small?
79        if (dataset_size < 10000) {
80            std::cout << "Dataset too small - use CPU\n";
81            return false;
82        }
83        
84        // Rule 2: Too much branching?
85        if (has_branching) {
86            std::cout << "Too much branching - use CPU\n";
87            return false;
88        }
89        
90        // Rule 3: Sequential dependencies?
91        if (has_dependencies) {
92            std::cout << "Sequential dependencies - use CPU\n";
93            return false;
94        }
95        
96        // Rule 4: Computation time vs transfer time
97        float transfer_time_ms = (dataset_size * sizeof(float) * 2) / (10.0 * 1e9) * 1000;  // 10 GB/s PCIe
98        
99        if (transfer_time_ms > computation_time_cpu_ms * 0.5) {
100            std::cout << "Transfer overhead too high - use CPU\n";
101            return false;
102        }
103        
104        std::cout << "Good candidate for GPU\n";
105        return true;
106    }
107};
108

1Monte Carlo Pricing (1M paths):
2- CPU (AVX-512):    10,200 ms
3- GPU (RTX 4090):   0.82 ms
4- Speedup:          12,439x
5
6Option Greeks (100k options):
7- CPU:              8,400 ms
8- GPU:              1.02 ms
9- Speedup:          8,235x
10
11VaR Calculation (50k positions, 10k scenarios):
12- CPU:              142,000 ms
13- GPU:              37 ms
14- Speedup:          3,838x
15
16Backtest (10 years daily data, 1000 simulations):
17- CPU:              2,840,000 ms (47 min)
18- GPU:              1,350 ms (1.35 sec)
19- Speedup:          2,104x
20

1Hardware:
2- CPU Server: Intel Xeon 8380 (40 cores) = $10,000
3- GPU Server: NVIDIA A100 (80GB) = $15,000
4
5Performance per Dollar:
6- CPU: 100 GFLOPS / $
7- GPU: 1,950 GFLOPS / $ (19.5x better)
8
9Power Efficiency:
10- CPU: 2.1 GFLOPS/Watt
11- GPU: 12.8 GFLOPS/Watt (6.1x better)
12
13Total Cost of Ownership (3 years):
14- CPU: $35,000 (hardware + power)
15- GPU: $22,000 (hardware + power)
16- Savings: 37% with GPU
17