用cuRand生成蒙特卡罗法所需的随机数.md

本文链接：https://blog.csdn.net/wxid2798226/article/details/86664973

背景

基于cuda的GPU异构计算为机器学习插上了翅膀。cuda实现的方法中，我最喜欢蒙特卡罗法（montecarlo方法），简单，“暴力”，但是非常简洁有效。

蒙特卡罗法（montecarlo），我给起了一个中文名字，叫***“雨滴法”***，比如，求pi的值，我们可以假设有一个1x1米见方的正方形面积，里面套了一个直径1米的圆，根据公式圆面积=pi*半径的平方，因为半径等于1/2米，那么pi就等于圆的面积*4。
那么怎么算圆的面积？我们已知正方形面积是1，假设有无数滴的雨滴，均匀落到正方形上，那么距离中心点在半径之内的肯定在圆上，反之则在圆外，假设有1亿点雨滴落下，其中8千万点落在圆内，那么圆的面积就等于0.8，然后通过圆的面积可以算出pi值

蒙特卡罗法，最关键的是生成随机数。比如，求pi的值，算出来是3.1415，问题是cpu生成随机数，然后拷贝到gpu，效率很低
考虑直接由gpu生成随机数并计算蒙特卡洛方法命中次数。
cuda中的cuRand可以生成类似普通C中的随机数，用于cuda蒙特卡罗法计算。

参考

curand最简单例子
如下我的实现，得益于如上例子，实现简洁但是不快
后来考虑的方法是积累法，过程中将每1.5亿次计算结果打印出来，后来看到所用的伪随机数明显有问题，要么一阵子偏大，要么一阵子偏小，用的是线性同余法，考虑更换
知乎讨论用于蒙特卡洛的随机数

/*
 * curand.cu
 *
 * Created on: 2019年1月14日
 * Author: vincent
 */

#include <stdio.h>
#include <stdlib.h>
#include <curand_kernel.h>
#include <time.h>
#include <sys/timeb.h>

#define BLK 12
#define THD 256
#define TOTAL BLK*THD
#define KERNEL 65536
#define PHASES 10
void print_t(void);

__global__ void generateRandom(long seed, unsigned short int* results) {
 unsigned short int hits = 0;
 curandState state;
 unsigned int id = blockIdx.x * blockDim.x + threadIdx.x;
 curand_init(seed, id, 0, &state);

 double x = curand_uniform_double(&state);
 for (int i = 0; i < KERNEL; i++) {
  double y = curand_uniform_double(&state);
  //double y = curand_uniform_double(&state);
  if (x * x + y * y < 1.0) {
   hits++;
  }
  x = y;
 }
 *(results + id) = hits;
 //printf("ID : %d, results:%u \n", id, hits);
}

int main()

{
 print_t();
 unsigned long sum = 0;

 unsigned int nBytes = BLK * THD * sizeof(unsigned short int);
 unsigned short int* results;
 cudaMallocManaged((void**) &results, nBytes);

 srand((unsigned int) time(NULL));
 cudaSetDevice(0);

 for (unsigned long n = 1; n <= PHASES; n++) {
  generateRandom<<<BLK, THD>>>(rand(), results);
  // 同步device 保证结果能正确访问
  cudaDeviceSynchronize();

  for (int i = 0; i < TOTAL; i++) {
   sum += *(results + i);
   //printf("下标：%d, 内容: %u \n",i, *(results+i));
  }
  double pi = ((double) (sum)) / (TOTAL * KERNEL * n) * 4;
  printf("finally n [%lu] sum[%lu] pi[%lf]\n", n, sum, pi);
  print_t();
 }
 cudaFree(results);
 cudaDeviceReset();
 return 0;
}

void print_t(void) {
 struct timeb tp1;
 struct tm *s_tm1;
 //记录交易时间
 ftime(&tp1);
 s_tm1 = localtime(&(tp1.time));
 printf("Start : %02d:%02d:%02d\n\n", s_tm1->tm_hour, s_tm1->tm_min,
   s_tm1->tm_sec);
}

现在GPU的速度是cpu的3倍,每秒2亿次计算，cpu大概只有7千万次；而且编程简单，不用搞pthread多线程；不占用界面操作，后台开着即可

几点经验：