curand_gpu_random_monte_carlo_pi_estimation.py

1import numpy as np
2from cuda import cuda, curand
3
4def check_cuda_error(result):
5    if isinstance(result, cuda.CUresult):
6        if result != cuda.CUresult.CUDA_SUCCESS:
7            raise RuntimeError(f"CUDA Error: {result}")
8    elif isinstance(result, curand.curandStatus):
9        if result != curand.curandStatus.CURAND_STATUS_SUCCESS:
10            raise RuntimeError(f"cuRAND Error: {result}")
11
12def main():
13    n = 1000000  # Number of points
14    
15    # Initialize CUDA Driver API
16    check_cuda_error(cuda.cuInit(0))
17    device_ptr = cuda.cuDeviceGet(0)[1]
18    context = cuda.cuCtxCreate(0, device_ptr)[1]
19
20    # Allocate memory on GPU
21    size = n * np.dtype(np.float32).itemsize
22    d_output = cuda.cuMemAlloc(size)[1]
23
24    # Create cuRAND generator
25    generator = curand.curandCreateGenerator(curand.curandRngType.CURAND_RNG_PSEUDO_DEFAULT)[1]
26    
27    # Set seed
28    check_cuda_error(curand.curandSetPseudoRandomGeneratorSeed(generator, 1234))
29
30    # Generate n float32 random numbers between 0.0 and 1.0
31    check_cuda_error(curand.curandGenerateUniform(generator, d_output, n))
32
33    # Copy results back to host
34    h_output = np.zeros(n, dtype=np.float32)
35    check_cuda_error(cuda.cuMemcpyDtoH(h_output, d_output, size))
36
37    # Simple Pi estimation logic (using first half as X, second half as Y)
38    half = n // 2
39    x = h_output[:half]
40    y = h_output[half:]
41    inside_circle = np.sum(x**2 + y**2 <= 1.0)
42    pi_estimate = 4.0 * inside_circle / half
43
44    print(f"Estimated Pi: {pi_estimate}")
45
46    # Cleanup
47    curand.curandDestroyGenerator(generator)
48    cuda.cuMemFree(d_output)
49    cuda.cuCtxDestroy(context)
50
51if __name__ == "__main__":
52    main()