MPI小试牛刀

有朋友问我一个Fortran程序的问题,Fortran语言以前曾学过,只是时间上太过久远,所有的语法全还给老师了。看着类似天书的代码,不禁想在单机系统上已经有如Matlab这样的如此强大的数值计算软件为什么还要费力用Fortran来实现呢?!在多机或集群环境上可能会不一样,只是我对并行计算所知甚少。

集群环境上MPI使用相当多,也只是了解一些,从未实际使用过,正好趁现在有些兴趣便做了做实验。

在网上看到有个计算Pi的程序 <http://chpc.wustl.edu/mpi-fortran.html>,用Fortran写的,其实就是计算定积分:

image

文中采用的是蒙特卡洛方法,我重写了一个C语言的实现,程序MPI_mc_pi.c,并在VMware环境下做了做实验,使用了2E9个采样点,分别做了单机不同并行进程数及双结点的对比,随并行度的增加,所需时间基本上是线性减少,由于测试机是双核四线程的Intel i5 CPU,当并行度超过4时便没有什么意义了,纯粹是为测试着玩。

后面又用差分的办法重新做了同样的计算,将[0, 1)分成2E9个网格,每个进程平分空间[0, 1),如进程1为[0, 1/nprocs),进程0处理最后剩下的部分[1-(nprocs-1)/nprocs, 1),见程序MPI_dc_pi.c,计算结果基本遵循线性减少规律,但是时间上却要比蒙特卡洛方法高出很多,特别是单进程的情况下,效率差了近50%。

于是又改进了差分程序,将空间平均分布,改成各进程散布在[0,1)空间上,即每个进程在网格采样上心nprocs为步进,参见程序MPI_dc_pi2.c,如此改进后,结果和均分采样空间基本差不多,没有显著变化。

后面又用双机做了Beowulf Cluster测试(VMware + X61),每个节点限制为2个进程,因为另一台机器是双核的CPU(X61),发现双节点情况下性能明显好于单节点。

将结果汇总后做成图表如下:

MPI-mc_vs_dc

如图所示,两种差分方式基本一致,这点完全可以理解,但是蒙特卡洛方法为什么能快这么多,确实让人费解,难道是drand48()的效率比double型除法还快?想起来总觉不太可能,那就让数据说话吧:

分别对drand48()和double型除法做了测试,结果正如预料,drand48()的效率更差,分别执行1E10次的时间为:

drand48():         4m21.065s double型除法: 3m27.399s

排除了这个原因之后,剩下的可能只能在drand48()里面了,drand48()虽是高精度的48位的伪随机数产生器,其结果应该还是不够平均和随机。评估drand48()的事咱不干,所以就做了个简单实验,只要能验证计算时间和drand48()生成的随机数不平均度显著相关即可。验证程序还是采用蒙特卡洛方法,只是将随机数的产生范围限定在[0.5,1]之内,不是之前的[0, 1],结果单进程计算时间为1m59.160s,双进程为1m35.185s,这足以说明随机数的随机度对计算时间有较大影响(Pi的值对与不对没关系,此时只关注计算过程所花的时间)。

最后在Mathematica上也做了个测试,出人意料的是Mathematica只用几了秒钟就算出了1E10网格的差分,即便是定积分计算,也是完完全全的秒杀,太牛x了,不服不行!

image

本文相关代码:

mpiuser@Ubunut-X64:~/MPI$ cat MPI_mc_pi.c #include <stdio.h> #include <stdlib.h> #include <time.h> #include <mpi.h> int main(int argc, char** argv) { int myrank, nprocs, rc; unsigned long t, times, avg_times; double l_sum = 0, g_sum = 0, x = 1.0; times = 2000000000; MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myrank); avg_times = times / nprocs; if (0 == myrank) { avg_times = times - avg_times * (nprocs - 1); } /* printf("%d: points: %lu\n", myrank, avg_times); */ srand48((myrank + 1) * (time(NULL) & 0xFFF)); for (t = 0; t < avg_times; t++) { x = drand48(); l_sum = l_sum + 4.0/(x * x + 1.0); } rc = MPI_Reduce(&l_sum, &g_sum, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD); if (0 != rc) { printf("MPI_Reduce failed with error code: %d\n", rc); } else if (0 == myrank) { printf("PI result: %f/%lu =  %f\n", g_sum, times, g_sum / times); } MPI_Finalize(); return 0; } mpiuser@Ubunut-X64:~/MPI$ cat MPI_dc_pi.c #include <stdio.h> #include <stdlib.h> #include <time.h> #include <mpi.h> int main(int argc, char** argv) { int myrank, nprocs, rc; unsigned long t, grids, avg_grids, start, end; double l_sum = 0, g_sum = 0, x = 1.0; grids = 2000000000; MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myrank); avg_grids = grids / nprocs; if (0 == myrank) { start = avg_grids * (nprocs - 1); end = grids; } else { start = avg_grids * (myrank - 1); end = start + avg_grids; } /* printf("%d: [%lu - %lu)\n", myrank, start, end); */ for (t = start; t < end; t++) { x = 1.0 * t / grids; l_sum = l_sum + 4.0/(x * x + 1.0); } rc = MPI_Reduce(&l_sum, &g_sum, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD); if (0 != rc) { printf("MPI_Reduce failed with error code: %d\n", rc); } else if (0 == myrank) { printf("PI result: %f with %lu grids.\n", g_sum/grids, grids); } MPI_Finalize(); return 0; } mpiuser@Ubunut-X64:~/MPI$ cat MPI_dc_pi2.c #include <stdio.h> #include <stdlib.h> #include <time.h> #include <mpi.h> int main(int argc, char** argv) { int myrank, nprocs, rc; unsigned long t, grids; double l_sum = 0, g_sum = 0, x = 1.0; grids = 2000000000; MPI_Init(&argc, &argv); MPI_Comm_size(MPI_COMM_WORLD, &nprocs); MPI_Comm_rank(MPI_COMM_WORLD, &myrank); for (t = myrank; t < grids; t += nprocs) { x = 1.0 * t / grids; l_sum = l_sum + 4.0/(x * x + 1.0); } rc = MPI_Reduce(&l_sum, &g_sum, 1, MPI_DOUBLE_PRECISION, MPI_SUM, 0, MPI_COMM_WORLD); if (0 != rc) { printf("MPI_Reduce failed with error code: %d\n", rc); } else if (0 == myrank) { printf("PI result: %f with %lu grids.\n", g_sum/grids, grids); } MPI_Finalize(); return 0; } 编译命令: mpiuser@Ubunut-X64:~/MPI$ cat build mpicc MPI_mc_pi.c -std=c99 -W -Wall -pedantic -D_SVID_SOURCE -g -o MPI_mc_pi mpicc MPI_dc_pi.c -std=c99 -W -Wall -pedantic -g -o MPI_dc_pi mpicc MPI_dc_pi2.c -std=c99 -W -Wall -pedantic -g -o MPI_dc_pi2
Continue reading » · Rating: · Written on: 02-26-13 · 1 Comment »

ERESOURCE:和你玩躲猫猫

我是一直都知晓你的存在的,虽然你从未现身,任何时候都像操控着一切的隐形之手,但在我心里,你是世界的主宰,规则的制定者,高高在上的上帝。

我只是一个小螺丝钉,每天只是重复地做着我被要求的工作。关于我的出生,我并不清楚,因为我的创造者-“瓶瓶”老爷向来沉默寡言,从不多说话,并且一直都是忙忙碌碌。我也只能在他有空的时候,和他有一句没一句得聊聊天,这才从他口中得知关于我身世一星半点的信息。

听“瓶瓶”老爷说,他是最不愿意制造我们这些螺丝钉的。他说,每制造出一个,他的负担就会增加了一个,虽然我们都是他的长工,每时每刻不停得为他工作,但我们的衣食住行还是要“瓶瓶”老爷为我们打点才行,这也是忙碌的他在空闲时经常向我们发的牢骚。不过,尽管他爱发些不着边的牢骚,我们还是很喜欢他,毕竟他创造了并一直照顾着我们,我们真心地愿意为他打长工,哪怕是一辈子,只要这个世界存在并运转着。

但就在前些日子,世界不在像以前那么可靠了,原因是那只隐形的上帝之手在规则上所做的一个很小的改动,改动虽小,对我来说却是灾难,从此我的世界整个都变成瓦蓝瓦蓝的了(Blue Screen Of Death)。生活如此简单的我,与世无争,咱这井水犯不着任何一条河水,真想不到哎,躺着也会中枪。

我只是重复着做一丁点儿的事,说来也很简单:以Paging I/O的方式去从文件里读数据出来。至于为什么要以Paging I/O方式,我曾问过“瓶瓶”老爷,他说的话很长很长,很深奥很深奥,我只记住了一点点,希望我的复述你能听得明白:

我要读的文件是于用户层中打开的,而我却生活在内核里,我的另一个哥们会将用户层的HANDLE转成FileObject并交给我,之后用户层会将这个HANDLE释放掉(CloseHandle),这时候我虽然手上有Fileobject这把钥匙,但通向Cached I/O及Direct I/O的路全给封上了,只有Paging I/O这条黑路让我去走,我的职责如此,即便是刀山火海,也只能认了。但自从世界变成Win7之后,我总会遇到麻烦,连Paging I/O这条黑道有时也不让我走了,特别是看起来更炫的Win8世界,任何时候都是瓦蓝瓦蓝的。

当时我很伤心,世界一片惨蓝,我以为再也见不着太阳了。可“瓶瓶”老爷又让我重获新生。所以我更感激他了,我发誓愿意为他打二辈子长工,最多只能二辈子,不能再多了,因为我们所处的这个世界一直是很二很二的。

在世界又恢复至正常后,“瓶瓶”老爷曾对我说,或者是老爷他自言自语地说:Win7世界的NTFS动了点手脚,在获取文件ERESOURCE锁后,它会改变此ERESOURCE的Owner Thread’s Pointer(参见ExSetResourceOwnerPointer)!虽然我知道在Filter中是不应该操作此ERESOURCE锁的,但针对这种情况,也是不得已而为之,做个简单的Paging I/O也不是什么危险或惹事的操作。

“瓶瓶”老爷现在让我获取两次共享锁并在释放锁时多做个判断,你想看看我最得意时的状态吗?

0: kd> dt _ERESOURCE 0xfffffa8002b62050; !locks -v 0xfffffa8002b62050; nt!_ERESOURCE +0x000 SystemResourcesList : _LIST_ENTRY [ 0xfffffa8002b620f8 - 0xfffffa80029030d0 ] +0x010 OwnerTable       : 0xfffffa80`03ba48c0 _OWNER_ENTRY +0x018 ActiveCount      : 1 +0x01a Flag             : 0 +0x020 SharedWaiters    : (null) +0x028 ExclusiveWaiters : (null) +0x030 OwnerEntry       : _OWNER_ENTRY +0x040 ActiveEntries    : 2 +0x044 ContentionCount  : 0 +0x048 NumberOfSharedWaiters : 0 +0x04c NumberOfExclusiveWaiters : 0 +0x050 Reserved2        : (null) +0x058 Address          : (null) +0x058 CreatorBackTraceIndex : 0 +0x060 SpinLock         : 0

Resource @ 0xfffffa8002b62050    Shared 2 owning threads      Threads: fffffa800196ab53-01<*> *** Actual Thread fffffa800196ab50

THREAD fffffa800196ab50  Cid 1364.1368  Teb: 000007fffffde000 Win32Thread: fffff900c229fc20 RUNNING on processor 0 IRP List: fffffa80042c2c40: (0006,0118) Flags: 00060043  Mdl: fffffa80045c7340 Not impersonating DeviceMap                 fffff8a00256f460 Owning Process            fffffa80019b9b30       Image:         notepad.exe Attached Process          N/A            Image:         N/A Wait Start TickCount      8907           Ticks: 0 Context Switch Count      219                 LargeStack UserTime                  00:00:00.000 KernelTime                00:00:00.592 Win32 Start Address 0x00000000ff383570 Stack Init fffff880171f1db0 Current fffff880171f0940 Base fffff880171f2000 Limit fffff880171e9000 Call 0 Priority 10 BasePriority 8 UnusualBoost 0 ForegroundBoost 2 IoPriority 2 PagePriority 5 Child-SP          RetAddr           Call Site fffff880171f17b0 fffff880159d69de Foobar!FbDirectReadWriteSingle+0x2a1 [Foobar\read.c @ 332] fffff880171f1840 fffff880159d73a8 Foobar!FbDirectReadWrite+0x85e [Foobar\read.c @ 628] fffff880171f1930 fffff80001cbb1b5 Foobar!FbDispatchRead+0x18 [Foobar\read.c @ 1172] fffff880171f1960 fffff80001cbac89 nt!IoPageRead+0x255 fffff880171f19f0 fffff80001ca165a nt!MiIssueHardFault+0x255 fffff880171f1ac0 fffff80001c920ee nt!MmAccessFault+0x146a fffff880171f1c20 00000000ff384061 nt!KiPageFault+0x16e (TrapFrame @ fffff880171f1c20) 00000000000fde80 00000000`00000000 0xff384061

              fffffa800196ab50-01<*>

THREAD fffffa800196ab50  Cid 1364.1368  Teb: 000007fffffde000 Win32Thread: fffff900c229fc20 ……

哈哈,不论你怎么躲,也躲不开我了,因为我知道了你所有的藏身之所。现在想改变规则了?你想怎么改呢?

我是个简单的小螺丝钉,我不惧怕任何权威,我吃苦耐劳,我要求很少,只要不让我的世界一片惨蓝,做再多的事我也愿意。

现在我的生命力更顽强了,我不仅会点自检,还会些自愈。面对这个美好的世界,只想说一句话,那就是:活着真好。

最后,我的id是“螺丝钉小P”。你的呢?我知道你有多个马甲的,嘿嘿!

摘自“螺丝钉小P”的日记(DEC/12/2012)

 

Continue reading » · Rating: · Written on: 02-01-13 · No Comments »

南山滑雪

密云南山滑雪场SKI-P1060953谢鹏,同事兼滑雪教练 :)SKI-P1060961SKI-P1060991SKI-P1060992SKI-P1060994铠信全家福

Continue reading » · Rating: · Written on: 02-01-13 · No Comments »

无题

一夜的风声呼啸
吹散了黑沉的阴云
带走了连日的雾霾

天,晴了
太阳出来了
就连月亮也笑盈盈得挂在
    西南方的蓝色天空上

美梦中醒来
仿佛又走进了另一个 美梦


在这样的好天气里
我在上班的路上

在这样的好天气里
我在拥挤的地铁里

在这样的好天气里
我坐在电脑前编写程序

上班的路上
清冽的寒风绕过脖颈扑向脊梁
让我想起夏天里的梦

拥挤的地铁里
人和人 贴得很近
但每一丝笑容里面 全都是陌生

在办公室里
我将程序编成了麻花
头和尾 缠绕在一起

望着窗外的阳光
开始想像着远方的夏天

只是不知远方的夏天
此时此刻
是否在期盼着我的到来


在这样的好天气里
我应该 做点什么

在这样的好天气里
我 应该 做点什么

Continue reading » · Rating: · Written on: 02-01-13 · No Comments »