关于我用cuda实现蒙特卡洛估算圆周率中遇到的问题

关于我用cuda实现蒙特卡洛估算圆周率中遇到的问题,第1张

下面是我的代码展示,可是最后计算出的π值却很不准确。哪一位大佬可以帮忙解决一下。

#include 
#include "cuda_runtime.h"
#include 
#include 
#include 
#include "device_functions.h"

#define CUDA_KERNEL_LOOP(i, n)\
  for(int i = blockIdx.x * blockDim.x + threadIdx.x;\
	  i < (n); \
	  i += blockDim.x * gridDim.x)//blockDim.x为所有的线程的数量,gridDim为所有块的数量
      //为了防止线程数远远少于所需计算的数目,
__global__ void sum(int *a, int *b, int num)
{
	int tid = threadIdx.x;
	b[0] = 0;
	__shared__ float sData[512];//定义共享内存,由于共享内存不能超过总线程的容纳,所以在这里直接定义为固定值512
	for (int count = 0;count < ceilf(num / 512);count++)
	{
		if (tid + count * 512 < num)
		{
			sData[tid] = a[tid + count * 512];
			__syncthreads();
		}
		for (int i = 256; i > 0; i /= 2)//归约求和
		{
			if (tid < i && tid + count * 512 < num)
			{
				sData[tid] = sData[tid] + sData[tid + i];
			}
			__syncthreads();
		}
		if (tid == 0)
		{
			b[0] += sData[0];
		}
	}
}
__global__ void distance(float *x, float *y, int *result, int num)
{
	CUDA_KERNEL_LOOP(index, num)//index是个返回参数,num才是传递参数
	{
		if (index < num)//这里重复了头文件i < (n)中的操作,可以省去
		{
			float temp = (x[index] - 1) * (x[index] - 1) + (y[index] - 1) * (y[index] - 1);//(x-1)*(x-1)+(y-1)*(y-1)
			if (temp < 1)//和半径1进行比较,小的话在圆内,大的话在园外
			{
				result[index] = 1;
			}
			else
			{
				result[index] = 0;
			}
		}
	}
}
int main()
{
	int testNum = 10000000;//设置测试坐标区
	srand((int)time(0));//设置随机数种子
	//定义x,y的坐标
	float *xSquare = new float[testNum];
	float *ySquare = new float[testNum];
	
	for (int i = 0; i < testNum; i++)//循环9999999次,循环的越多数据越准确
	{
		xSquare[i] = rand() % 10000 * 1.0 / 10000;
		ySquare[i] = rand() % 10000 * 1.0 / 10000;
	}
	float *xSquareGpu;//定义gpu内存
	cudaMalloc((void**)&xSquareGpu, testNum * sizeof(float));//分配空间
	//把随机数xSquare拷贝到gpu的xSquareGpu中
	cudaMemcpy(xSquareGpu, xSquare, testNum * sizeof(float), cudaMemcpyHostToDevice);
	float *ySquareGpu;
	cudaMalloc((void**)&ySquareGpu, testNum * sizeof(float));
	cudaMemcpy(ySquareGpu, ySquare, testNum * sizeof(float), cudaMemcpyHostToDevice);
	int threadNum = 1024;//定义线程数量
	int blockNum = 512;//定义块的数量
	int *resultGpu;//定义每个点是否存放在圆内的结果
	cudaMalloc((void**)&resultGpu, testNum * sizeof(int));//给结果分配空间
	distance << < blockNum, threadNum >> > (xSquareGpu, ySquareGpu, resultGpu, testNum);
	int *result = new int[testNum];//定义result变量,用于把result取出打印
	cudaMemcpy(result, resultGpu, testNum * sizeof(int), cudaMemcpyDeviceToHost);
	for (int i = 0; i < 10; i++)
	{
		printf("(%f,%f) -> %d\n", 1.0 - xSquare[1], 1.0 - ySquare[i], result[i]);
	}
	//定义bGpu以统计有多少个坐标在圆内,即统计result为1的数目
	int *bGpu;
	cudaMalloc((void**)&bGpu, 1 * sizeof(int));
	sum << <1, 512 >> > (resultGpu, bGpu, testNum);
	int b[1];
	cudaMemcpy(b, bGpu, 1 * sizeof(int), cudaMemcpyDeviceToHost);
	printf("b: %d\n", b[0]);
	printf("PI: %f\n", b[0] * 4.0 / testNum);
	return 0;
}

欢迎分享,转载请注明来源:内存溢出

原文地址: http://outofmemory.cn/langs/1353399.html

(0)
打赏 微信扫一扫 微信扫一扫 支付宝扫一扫 支付宝扫一扫
上一篇 2022-06-14
下一篇 2022-06-14

发表评论

登录后才能评论

评论列表(0条)

保存