CUDA果然好厉害！8过也各种无法理解ww

也许有人还记得咱曾经做过Mente Carlo法球圆周率的程序，当时那个是用的C++写的；现在为了练习CUDA，咱也用CUDA写了一个——8过现在这个程序运行倒是运行了，但是有一个奇怪的问题是：完全一模一样的程序，直接在Terminal里面用./打开得到的误差远高于直接双击图标打开运行得到的误差orz原因不明

如上图所示，同样是512×10⁶个随机点，左边是直接在Terminal里面用./打开得到的结果，得到的圆周率对比实际的圆周率的误差高达7825×10^-6，而右边是直接双击打开得到的结果，误差只有24×10^-6，而且这不是偶然，无论测试多少次./得到的结果都是高达4位数，而直接打开则一般都是两位数，或者运气好的时候也有一位数的时候orz可以肯定./打开的时候绝对有什么奇怪的效果，但是目前为止查了各种资料未果orz虽然大致可以推测应该是CUDA的随机数函数curand的问题ww CUDA到底是怎么回事www

啊当然，对于Monte Carlo法来说，GPU的运算要远快于CPU，同样是512×10⁶个随机点，咱女仆本上的9400m只要不到两秒钟，而Core 2 Duo（2.26GHz）则要将近一分钟www

嘛，这里把CUDA的代码贴上来，欢迎有兴趣的童鞋研究www

#include <stdio.h>
#include <stdlib.h>
#include <math.h>
#include <time.h>
#include <curand.h>
#include <curand_kernel.h>
#define MAX_THREAD  256
#define MAX_COUNT	2000000
/*	最大线程MAX_THREAD	256
	每个线程生成随机点个数MAX_COUNT	2000000	*/

/* GPU端圆周率函数piGenerator宣言 */
__global__ void piGenerator (float *res, int *count, int *time);

/* 主程序 */
int main (int argc, const char * argv[]){
	/* CPU端圆周率h_pi、GPU端圆周率d_pi、最终圆周率（平均值）pi、偏差err */
	float *h_pi, *d_pi, pi, err;
	
	/* 最大线程数maxThread，随机点个数h_count/d_count，随机数初期值设定用时间h_time/d_time，for循环次数i */
	int maxThread = MAX_THREAD, *h_count, *d_count, *h_time, *d_time, i;
	
	/* CPU端h_pi、h_count、h_time内存确保 */
	h_pi = (float *)malloc(sizeof(float) * maxThread);
	h_count = (int *)malloc(sizeof(int) * 1);
	h_time = (int *)malloc(sizeof(int) * 1);
	
	/* GPU端d_pi、d_count、d_time内存确保 */
	cudaMalloc((void **)&d_pi, sizeof(float) * maxThread);
	cudaMalloc((void **)&d_count, sizeof(int) * 1);
	cudaMalloc((void **)&d_time, sizeof(int) * 1);
	
	/* CPU端随机点个数值h_count及时间值h_time带入 */
	h_count[0] = MAX_COUNT;
	h_time[0] = (int)time(NULL);
	
	/* 将h_count及h_time的值带入GPU端d_count及d_time */
	cudaMemcpy(d_count, h_count, sizeof(int) * 1, cudaMemcpyHostToDevice);
	cudaMemcpy(d_time, h_time, sizeof(int) * 1, cudaMemcpyHostToDevice);
	
	/* 运行时间测试 */
	clock_t start, end;
	
	/* 运行时间测试开始 */
	start = clock();
	
	/* GPU端函数piGenerator开始运行，1个block，maxThread = 512个thread */
	piGenerator<<<1, maxThread>>> (d_pi, d_count, d_time);
	
	/* 将GPU生成的圆周率d_pi带入CPU端c_pi */
	cudaMemcpy(h_pi, d_pi, sizeof(float) * maxThread, cudaMemcpyDeviceToHost);
	
	/* 512个线程生成的圆周率求平均值 */
	for (i = 0; i < maxThread; i++) {
		pi += h_pi[i];
	}
	pi = pi / maxThread;
	
	/* 运行时间测试停止 */
	end = clock();
	
	/* 求圆周率误差 */
	err = pi - (float)M_PI;
	if (err < 0) {
		err = -err;
	}
	
	/* 结果输出 */
	printf("Points: %d, Generated π: %f, Error: %.0fe-6\n", h_count[0] * maxThread, pi, err * 1000000);
	printf("Timer: %f sec\n", (float)(end - start)/CLOCKS_PER_SEC);
	
	/* CPU端内存释放 */
	free(h_pi);
	free(h_count);
	free(h_time);
	
	/* GPU端内存释放 */
	cudaFree(d_pi);
	cudaFree(d_count);
	cudaFree(d_time);
	
	/* 结束 */
    return 0;
}

__global__ void piGenerator (float *res, int *count, int *time) {
	/* 1/4圆内点个数a初期值=0，线程ID index_x=threadIdx.x，随机点个数n带入，时间值t带入，for循环次数i */
	int a=0, index_x = threadIdx.x, n = *count, t = *time, i;
	
	/* 随机点x、y宣言 */
	float x, y;
	
	/* 圆周率内存保存地址确认 */
	res += index_x;
	
	/* curand随机函数初期化 */
	curandState s;
	curand_init(t, index_x, 0, &s);
	
	/* for循环次数=随机点生成个数 */
	for (i = 1; i <= n; i++) {
		/* 随机点x坐标、y坐标（0<x<1，0<y<1）的随机生成 */
		x = curand_uniform(&s);
		y = curand_uniform(&s);
		
		/* 如果随机点在1/4圆圈以内，1/4圆内点个数a+1 */
		if (pow(x, 2) + pow(y, 2) <= 1) {
			a++;
		}
		
		/* 圆周率计算（1/4圆内点个数 / 随机点总数 * 4） */
		*res = 4 * (float)a / (float)n;
		
		/* 线程同步 */
		__syncthreads();
	}
}

嘛CUDA的事就说到这儿w话说最近宇都宫的大型电器店Yodobashi Camera的店内Layout稍微改了一下，然后似乎成为了BOSE的授权零售商，于是音响区增加了一个BOSE专区！好想要Companion 20啊＞＜好想要Companion 5啊＞＜

话说前几天购物回家的时候在等红灯的时候看到旁边的停车场里面的车——人生＼(^o^)／ｵﾜﾀ!!

大天使千反田恵瑠＞ω＜

为了迎接ゆるゆり♪♪而专程购入的ラムレーズン！ラムレーズンうめぇ！！

以上XD

p.s. 虽然很多人都称偶为大大——8过咱真心不觉得咱是大大，至少咱现在还远远达不到咱心里“大大”的标准——虽然咱正在努力成为“大大”www 8过至少现在来说咱真的不是大大！！！！！wwwww

Author: 星野恵瑠

Mac user, Niji-Ota, Chinese, Now working in Japan at MAGES. Inc., Future's aim is that one day my name can be listed in Wikipedia View all posts by 星野恵瑠