MC calculation of Pi (OpenACC+CURAND)

pi_curand.cc

#include <iostream>
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <curand.h>
//----------------------------------------------------------------------
const int N = 10000;
const int SAMPLE = 65536;
//----------------------------------------------------------------------
double
myrand(void){
  return static_cast<double>(rand())/static_cast<double>(RAND_MAX);
}
//----------------------------------------------------------------------
double
one_trial(void){
  int sum = 0;
  for(int i=0;i<N;i++){
    const double x = myrand();
    const double y = myrand();
    if (x*x+y*y < 1.0) sum ++;
  }
  return static_cast<double>(sum)/static_cast<double>(N)*4.0;
}
//----------------------------------------------------------------------
double r[2*N];
//----------------------------------------------------------------------
double
one_trial_acc(){
  int sum = 0;
#pragma acc kernels pcopyin(r[1:2*N])
#pragma acc loop reduction(+:sum)
  for(int i=0;i<N;i++){
    const double x = r[i*2];
    const double y = r[i*2+1];
    if (x*x+y*y < 1.0) sum ++;
  }
  return static_cast<double>(sum)/static_cast<double>(N)*4.0;
}
//----------------------------------------------------------------------
int
main(void){
  double a = 0.0;
  double v = 0.0;
  double *d_r;
  curandGenerator_t g;
  cudaSetDeviceFlags(cudaDeviceMapHost);
  curandCreateGenerator(&g, CURAND_RNG_PSEUDO_DEFAULT);
  curandSetPseudoRandomGeneratorSeed(g,0);
  size_t size = sizeof(double)*(N*2);
  cudaMalloc((void**)&d_r, size);
  for(int i=0;i<SAMPLE;i++){
    curandGenerateUniformDouble(g, d_r, N*2);
    cudaMemcpy(r, d_r, size, cudaMemcpyDeviceToHost);
    /*
    for(int i=0;i<2*N;i++){
      r[i] = myrand();
    }
    */
    const double p = one_trial_acc();
    a += p;
    v += p*p;
  }
  const double s = static_cast<double>(SAMPLE);
  a /= s;
  v /= s;
  v = sqrt((v -a*a)/(s-1.0));
  std::cout << "Total Samples " << N*SAMPLE << std::endl;
  std::cout << a << std::endl;
  std::cout << v << std::endl;
  curandDestroyGenerator(g);
  cudaFree(d_r);
}