pi_cuda.cu

#include <iostream>
#include <cuda.h>
#include <algorithm>
#include <curand_kernel.h>

const int N = 10000;

__global__ void
func(double *out){
  int id = blockIdx.x * blockDim.x + threadIdx.x;
  curandState s;
  curand_init(0,id,0,&s);
  int sum = 0;
  for(int i=0;i<N;i++){
    double x = curand_uniform_double(&s);
    double y = curand_uniform_double(&s);
    if(x*x + y*y < 1.0)sum++;
  }
  out[id] = static_cast<double>(sum)/static_cast<double>(N)*4;
}

int
main(void){
  const int n_thread = 256;
  const int n_block = 256;
  const int n = n_thread * n_block;
  static double r[n];
  std::fill_n(r,n,0);
  double *d_r;
  cudaSetDeviceFlags(cudaDeviceMapHost);
  size_t size = n*sizeof(double);
  cudaMalloc((void**)&d_r,size);
  func<<<n_block,n_thread>>>(d_r);
  cudaThreadSynchronize();
  cudaMemcpy(r,d_r,size,cudaMemcpyDeviceToHost);
  double a = 0.0;
  double v = 0.0;
  for(int i=0;i<n;i++){
    a += r[i];
    v += r[i]*r[i];
  }
  const double s = static_cast<double>(n);
  a /= s;
  v /= s;
  v = sqrt((v- a*a)/(s-1.0));
  std::cout << "Total Samples " << N*n << std::endl;
  std::cout << a << std::endl;
  std::cout << v << std::endl;
}