蒙特卡罗算法实验报告

六、
为了改进并行算法，得到更高的加速比，有两种途径可以尝试：减少线程状态转化次数和使用可并行的随机数产生算法。简介如下：
6.1
此方法具体为：在并行程序中使用互斥锁，当某一线程进入临界区后，一次性产生m个随机点，然后再退出临界区，开始对m个点进行计算；与此同时，若另一线程也要进入临界区，则被挂起，等待该线程退出。如此循环，直至两个线程均计算完所要求的点的个数，则计算输出 值，程序结束。
7、通过（2）式计算 的值。
3.2
在这个实验中，采用Linux操作系统pthread接口来实现程序的并行化。这些接口函数和数据类型都在头文件<pthread.h>中声明。因为pthread并没有包含在C的标准库中，编译的时候需要加上-lpthread选项，使程序链接到libpthread，才能编译成功。
return (double)Y/(double)(d-1);
}
通过改变种子的值，算法可生成不同的伪随机数列并且可以满足多个处理器同时调用。但调用所需时间略大于调用系统库函数rand（）。(调用myrand()函数的串行算法，见附件Smyrand.c)
示例程序见附件Pmyrand.c
20
22.033
20.757
20.120
20.647
19.918
20.798
20.160
50
49.725
49.785
51.535
49.420
50.992
52.379
47.015
表2
图4
而对同样的运算量多次运行串行算法得到如下表3所示结果。（图5）
规模
1
2
3
4
5
6
7
1
0.814
0.814
0.814
0.813
示例见附件Serial.c
三、
3.1
算法步骤：
1、确定需要产生的点的个数n，参与运行的处理器数m；
2、对每一个处理器，生成两个随机数x，y，范围[0，1]；
3、判断两个随机数x，y是否满足 ；
4、若满足，则变量COUNTi++；
5、重复步骤2-4，直至每个处理器均生成n/m个随机点；
6、收集COUNTi的值，并累加至变量COUNT中，此即为随机点落在圆弧内的数量；
0.810
0.815
0.813
5
4.024Biblioteka Baidu
4.053
4.062
4.057
4.044
4.090
4.053
10
8.089
8.152
8.153
8.134
8.160
8.095
8.108
20
16.121
16.289
16.290
16.318
16.288
16.245
16.292
50
40.726
40.721
40.701
示例程序参见附件Pmutex.c。
6.2
生成随机数最常用的方法为线性同余法，其C语言源代码如下：
//myrand()用的种子
unsigned static Y =568731;
unsigned d=1<<31;
//生成伪随机数算法
double inline myrand()
{
Y=(15625*Y+22221)%d;
0.8488
2.5
1.618
2.673
0.6053
5
4.024
5.716
0.7039
7.5
6.069
7.376
0.8228
10
8.089
10.001
0.8088
12.5
10.105
12.227
0.8264
15
12.115
14.842
0.8162
17.5
14.119
19.522
0.7232
20
16.121
对同一运算量多次运行并行算法得到如下表2所示结果。（图4）
规模
1
2
3
4
5
6
7
1
0.959
1.205
0.963
1.043
1.002
1.053
1.011
5
5.716
4.877
5.094
4.761
5.212
4.875
5.296
10
10.001
9.892
9.990
10.151
9.941
10.168
10.169
例子程序参见附件Parallel.c。
3.3
本并行算法只是简单的把独立的任务进行分派，经多次试验测试，结果正确。
四、
硬件平台：惠普刀片集群
编译器：gcc&g++
操作系统：Linux
测试数据集合：由随机数函数产生的数据集合
4.1
N（千万）
串行算法运行时间（秒）
并行算法运行时间（秒）
加速比
1
0.814
0.959
40.696
40.706
40.695
40.694
表3
图5
如图4图5所示，对同一计算量，串行算法每次运行时间相差较小，而并行算法则相差明显。因此，通过分析源代码可得出以下结论：
程序所用的rand（）函数在同一时间只允许一个处理器调用，当两个处理器都需调用rand（）函数时，后调用的将被挂起，等待另一个处理器运行完毕。两线程在就绪和执行态之间不断变化，浪费了大量CPU时间，因此对同一运算量，并行程序运行时间反而比串行程序慢，而且线程状态转换次数范围为[0，n]，平均为 次，因此，相比于串行程序的无状态转换，并行算法的运行时间才会有如此大的波动。
算法：
1、确定产生点n的个数和缓冲区m（m<=n）的值，声明互斥锁
2、某一线程进入临界区，上锁
3、该线程一次性生成m个数，其他线程若想进入则挂起等待
4、该线程退出临界区，解锁，开始对刚才生成的随机点进行计算
5、重复2-4步，直至每个线程均完成对所要求点的操作
6、统计COUNTi的值
7、计算 的值
在此算法中，每一线程因为争用rand（）函数而产生的状态转化次数范围为[0， ]，平均次数为 ，调整m的值，使生成m个随机点的时间与对m个随机点进行计算的时间相等时，则算法执行速度可达到最大值，即加速比最大。
intrand(void);
此函数产生随机数列，每次调用时均返回0到RAND_MAX之间的一个整数。
voidsrand(unsignedintseed);
此函数为rand（）函数所生成的伪随机数序列设置起始点，使之产生不同的伪随机数。
算法：
产生2n个随机数据，范围[0，1]，对每个数据点计算其坐标是否满足 ，统计满足此关系的点的数量count，则
（1）
由此可得
（2）
因此，只要计算出落在圆弧内的点的数量在点总数中所占的比例，就能求出 的值。
由图1可知，所有点均落在正方形范围内，因此点的x坐标满足 。又，当点落在圆弧范围内，则点的二维坐标关系满足 。检验每一个点是否满足此关系即可判定改点是否落在圆弧内。
二、
本项目中使用了标准C语言库中的产生随机数函数。该函数原型为：
多核软件设计——实验指导
蒙特卡洛算法
开发者：
开发时间：
版本号：
蒙特卡洛算法可理解为通过大量实验，模拟实际行为，来收集统计数据。本例中，算法随机产生一系列点，模拟这些点落在如下图所示的正方形区域内的情况。其几何解释如下
图1
如图1所示，正方形边长为1，左下顶点与原点重合，两边分别与x，y轴重合。曲线为1/4圆弧，圆心位于原点，与正方形左下定点重合，半径为1。正方形面积S1=1，圆弧内面积S2= 。算法模拟大量点随机落在此正方形区域内，落在圆弧内的点的数量（n2）与点的总数（n1）的比例与面积成正比关系。即
22.033
0.7316
30
24.183
32.592
0.7419
40
32.259
41.542
0.7765
50
40.726
49.725
0.8109
表1
[注]：N：算法生成随机点的个数
算法运行时间为某一次运行时间，非多次运行之平均时间
4.2
并行、串行算法运算量时间比、加速比如下图所示
图2
图3
五、
如表1、图3所示，加速比在（0.6,0.9）区间，与理论上的值2相去甚远。