即用计算机产生统计独立的、有代表性的、对求平均贡献最大的组态，满足Boltzmann分布

即用计算机产生统计独立的、有代表

性的、对求平均贡献最大的组态，满足Boltzmann分

布

第28卷第2期

2004年2月

高能物理与核物理

HIGH ENERGY PHYSICS AND NUCLEAR PHYSICS

Vo1．28，No．2

Feb．，2004

Ising模型的并行计算*

刘军 沈扬

(中山大学物理系

罗向前D

广州 510275)

摘要 以Ising模型为例，介绍有关格点系统的Monte Carlo数值模拟并行算法的设计和编程

，

并给

出在本组建造的PC集群式高性能并行计算系统上的测量结果．本文的结果对格点量子色动力学

的大规模数值模拟研究有一定的参考价值．

关键词 蒙特卡罗模拟 Ising模型 格点规范理论 并行算法 高性能并行计算系统

1 引言

最近获批准的国家自然科学基金重点项目“格

点规范理论的大规模数值模拟研究⋯，成员由中山

大学、北京大学、浙江大学和中科院高能物理研究所

等单位组成．目的是利用国内日益发展的高性能计

算系统，从第一原理出发，探索计算非微扰物理量的

新方法，对胶球和混杂态等新强子态性质、强子散射

过程、规范场真空拓扑性质、夸克禁闭机制、夸克胶

子等离子体相变、量子瞬子、量子混沌等热门课题开

展格点规范理论的数值模拟研究．本项目的主要实

验工具是大规模计算机模拟．单靠一个或多个电脑

独立工作，无法处理更大体积的格点，需要用到高性

能计算机和大量计算机时，而并行计算能较好地解

决这些问题，

本文以二维Ising模型为例，介绍数值模拟及其

并行化程序的分析、设计和实现，同时简单介绍格点

量子色动力学的基本思想和在本组的高性能并行计

算系统一些测试结果，

2 Ising模型和Monte Carlo模拟

Ising模型 是在统计物理中最简单的格点系统

之一．二维Ising模型有精确解，可以用此来作为基

准检验其他近似方法的可靠性．并着重讨论如何用

Monte Carlo模拟方法来研究并如何提高其计算效

率．要得到更加精确的结果需要很大数目的格点，

这就需要并行计算机．

二维的Ising模型可描述为：设有V=L x L个

自旋，处于晶格格点位置，每个自旋只能取向上或向

下两个态，如图1所示，并只考虑近邻自旋之间的相

互作用，其哈密顿量为

V

H=一_，ΣS 一／zBΣS ， (1)

!

⋯

0

⋯ ⋯

1

⋯ 一 ⋯⋯

3

⋯ ．

图1 二维4×4的格点和自旋组态

其中si代表第i个格点位置的自旋，取值为+1或

一

1，分别对应于自旋向上或向下．_，为耦合常数．

这里令J>0，H描述铁磁体(因 为B=0时，其基态

是全部自旋取向相同的态)．若J<0，则描述反铁

磁体．式中第二项代表在外磁场B中的塞曼能，

2003—05—26收稿，2003—07—27收修改稿

*国家自然科学基金(10235040)，广东省教育厅，中山大学高等学术中心基金资助

1)通讯联系人，E-mail：stslxq@zsu．edu．On

l22— 128

第2期 刘军等：Ising模型的并行计算

为与自旋相应的磁矩．真正的物理系统对应于热力

学极限 一∞的情形，但实际计算只能对有限的

进行．为了考虑日中的第一项在边界的效应，一般

采用周期边界条件．

统计物理通常要计算一些物理量的统计平均

值，如系统内能和比热：

U =(日)，

， 、

a U (2)

L ，

其中某一物理量A的统计平均值为(A)定义为

(3)

以上公式中，分母是配分函数， 是温度， =

1／(k )，k 是Boltzmann常数，而求和对所有可能

的自旋组态进行．

对于 个格点的系统，这样的求和的项就有

2 ，即随着格点的大小指数地增长，要产生所有这些

组态可能远超出所有最先进的计算设备的能力．如

果用计算机求解，最有效是采用Monte Carlo重点抽

样法来计算．即用计算机产生统计独立的、有代表

性的、对求平均贡献最大的组态，满足Boltzmann分

布

s ． (4)

从平衡态开始，抽出其中J7、r个组态并测量有关物理

量．这样，A的统计平均值可近似为

(A)= ΣA({s} )． (5)

Metropolis抽样方法b 是实现产生组态的重点抽样

法之一，除此之外，还有热浴法(Heatbath)和杂化

(Hybrid)Monte Carlo法等．Metropolis算法的基本思

想是，通过一个适当的跃迁概率 ({s} 一

{S} + )：

w(I S I s 1)-min[1， 】'(6)

其中6日=H({s} + )一日({s} )．可以证明，按上

式满足细致平衡条件．以上过程的伪码见附录A．

3 格点规范理论

粒子物理和规范场论是探索最深层次物质世界

的前沿科学，量子色动力学(QCD)是描述物质结构

最基本的单元、即夸克与胶子间强作用的规范理论．

QCD的特点是低能区相互作用非常强，微扰论完全

失效，只能用非微扰方法处理．由诺贝尔奖得主

Wilson所建立的格点规范理论，是处理强作用非微

扰的最可靠工具．它还与计算科学紧密结合起来，

推动物理学和其他学科的发展．格点规范理论的基

本思想是：把四维连续时空用离散的晶格代替，这

样，连续理论无法计算(发散的)物理量就可在格点

上作Monte Carlo数值模拟计算，最后还要用重整化

群方法把结果外推到连续极限．当然这个工作量是

非常巨大的，要用到高速并行计算机．以纯规范场

为例，把规范场定义连接格点间的链上

r

。

1

( ， )=exPl_ig 。 l dx A ( ，)1． (7)

。

。

格点之间的距离是o(晶格常数)．格点位置为 ，

为链的方向．其Wilson作用量为

s =一卢 P (8)

其中P =Tr(up+h．c．)／6， 绕在格点位置为 ，

方向为 和 的1×1方块的4条规范场链的乘积．

=

1／g ，而g是规范场的耦合常数．可以证明，在

连续极限0—0，以上作用量回复到Yang-Mills作用

量．

规范场的配分函数是

z=JI [dU]exp(一S )． (9)

像Ising模型那样，可以用重点抽样法(较常用的是

热浴法和杂化Monte Carlo法)来产生满足Boltzmann

分布exp(一S )／Z的规范链组态，然后对各种物理

量进行测量．

当然，格点QCD的方法有系统的误差．众所周

知，Wilson规范场作用量与Yang．Mills作用量之间有

O(o )的差别．只有。非常小时，这个误差才会变

小．为了消除有限体积效应，需要很大的格点体积

，

这无疑会耗费更多的计算资源．Symanzik改进方

案 的想法是：通过把次近邻的项加入Wilson作用量

来减小格距误差，例如，一个最简单的改进作用量为

s =一卢Σ(詈 ， 一 ， )， (10)

其中尺 对应于2×1和1×2的长方形的6条规范

场链的乘积的项．

与Tadpole改进 或非微扰改进n 结合起来，

Symanzik改进QCD方案已导致十分有意义的进展，

使在粗糙的格点上获得趋近物理连续极限的物理结

果成为可能．

一"

二二e= 一面“

S 一一e

Σ__

=

A

124 高能物理与核物理(HEP＆NP) 第28卷

4 并行算法的设计

并行算法是指一个过程可分配到不同的子进程

在不同的处理器上同时处理的计算方法．以上格点

系统十分容易并行化：即可把大格点划分成几块，分

别分配给不同的子进程执行组态更新，最后再把更

新了的数据归纳于一起进行测量．现在可以初步预

见，对于较大的格点，如果通信／计算时间比是较小

的，则可使用并行计算将大大提高计算时间．

主体算法描述如图2．以二维Ising模型为例，

根据其特点，可采用以下并行算法．先在主任务中

把二维格点划分成数条 ]，然后分配到各个处理器

或进程中去，每个处理器或进程再在自己的内存空

间中执行分配给它的子任务，每一子任务代表完成

对某一个分条块作一次子组态的更新，并且只测量

此条格点的数据，如磁场强度和内能．最后由MPI

函数MPI．ALLReduce将所有的处理器或进程中的数

据归纳，得到整个格点的总体数据．

循

图2 并行算法主要流程

由于并行程序的特殊性，其算法设计远不同于

串行程序，主要难点在于任务划分，进程间通信，此

外还涉及网络拓扑、采用的语言及支持库等．这里

采用了MPICH作为并行开发环境．以下来看几个

技术要点．

条分解(Strip Decomposition)是把主任务分配给

子任务的方法 ．假设有P个处理器(分别标为0，

1，2⋯，P一1)，就可把大格点分成L／P个子格点，如

图3所示，其中每块大小为L／P×L，分配给处理器

P 的处理单元为从第(L／P)×a行到[( ／P)×(a

+1)]一1行．图中 =l6，P=4，口=0，1，2，3．

0 2 4 6 8 10 l2 l4

J

图3 条分解

(1)式中的日描述整个格点上所有相邻自旋的

相互作用，要把日对所有子格点求和．在划分块之

后，每块子格点的头尾两行要从其他进程(处理器)

中的相应行取得数据，要通过进程间通信来实现．

此处采用进程间缓存(Interprocessor Caching)的技术

来解决．在子格点的头尾两行上再加两行缓存行

(Caching Row)，用来接收和发送数据给近邻的进程，

如图4所示．为了确切自身的进程号及其上下相邻

的进程号，每个进程需要设置以下变量：

int pid一自身的进程号；

im up—id=(pid+1)％ P一处理上一块子格点

的进程号；

int lw_id=(pid+P一1)％ P一处理下一块子格

点的进程号．

图4 缓存模型

因此，设置子格点的自旋数组大小为s[(L／P)

+2][ ]，多出来的两行用来放置缓存行．

5 编程

对于二维Ising模型，用C语言编程，并用MPI

把程序并行化．以下给出程序的主要思路，具体源

码在附录B，C，D，E给出．

程序由主函数main，子函数simulate，update组

成．其中主进程主要执行任务分发，最后将所有子

进程的数据集中起来进行下一步的数量分析．各子

进程在接受到所分配的子任务后将读取各自所需的

数据(readData函数)，然后初始化各自的格点自旋

初态，接着开始调用simulate子函数实行模拟过程．

第2期 刘军等：Ising模型的并行计算

其中涉及到对自旋组态的多次更新，这里通过多次

调用update子函数实现对格点自旋组态的更新．

Update子函数中调用lOW．update来实现Metropolis算

法．由于在扫描格点的首尾两行时需要相邻进程的

数据，所以这里需要用到上面所说的“进程间缓存”

技术来实现，具体算法可以参见附录．

各进程完成各自的数据处理后就由请进程将数

据归约起来进行数据分析，由主进程调用measure子

函数来实现．

6 结果分析

2000年初，本组建造了一套有20个CPU的PC

集群式高性能并行计算系统 引．其基本配置为

(1)硬件配置： 。

10个节点，每个节点有双CPU PIII 500MHz；

每个节点内存：128MB；

交换机速度：100MBbit／s，

(2)软件环境：

操作系统：RedHat Linux；

并行计算环境：MPICH．

对于k =J=1，B=0的二维Ising模型，测量

了系统的平均内能和平均比热等物理量．结果如图

5所示．从结果同准确解比较可以看出，计算机模拟

和准确解符合较好．(测量次数1000，组态间隔200，

格点大小60×60)．

一

0．8

—

1．2

0 l 2 3 4 5 6

r

图5

(a)内能随温度的变化；(b)比热随温度的变化

图6显示二维Ising模型计算机时与格点大小、

所用CPU数目的关系．由图可以看出，在格点总数

V=L×L不大的情况下，并行计算的优势并不能很

好地体现出来．如 =20，40和60的情况下，用两

个CPU计算只比用一个CPU计算的时间略快，加速

比在1．3—1．9附近，当增加CPU数时，计算时间并

没有多大的缩短，甚至有增长的情况出现(如 =20

时)．这是由于￡不大的时候，每个CPU的计算量

较少，通信时间占了总的计算时间的较大比例，所以

加速比不好．当体积增大到一定程度，例如上图中

的￡≥100时，并行计算的优势就能较好地体现出

来．如用两个CPU的加速比很接近于2．用多个

CPU计算时，时间也有较大的减少．这是因为当

CPU间的通信时间比每个CPU的计算时间小得多

时，并行计算的优势才能很好地体现出来加速比也

会上升．

并行计算结果

cPu数目

图6 Ising模型计算时间与CPU数目的关系

以上Ising模型的计算方法和结果对格点规范

理论的大规模数值模拟研究很有启发作用．格点

QCD的并行算法也是把大格点分成多个子格点在

不同的CPU运算．在子格点的边界，必然有信息交

换．幸好，这种信息交换量并不大．在本组的并行

计算系统上，用(10)式给出的最简单的改进纯规范

场作用量进行Monte Carlo模拟，格点体积V=10 ．测

表1 本组计算系统的改进格点QCD的运行结果．并与

其他计算系统比较

⋯

更新每条链所花机时：节点间的通信速度

机器 链

Mbits／s

∞ 如 ∞ 如 O

2 I l

∞＼匡营琳4{

6 O

l 2

126 高能物理与核物理(HEP＆NP) 第28卷

试了更新每条链所花的计算机时和节点间的通信速

度．表1给出本组的测试数据，并与其他系统比

较 ．还在本组的系统中，运行MILC的格点QCD

并行计算程序库 ，结果是令人鼓舞的¨ ．

我们已经利用本系统，用格点QCD计算了轻强

子¨ 、胶球n 和混杂态ll 的质量，这些结果和经

验，对利用国内更高性能并行计算系统，进行更大规

参

模的格点规范理论数值模拟研究，相信是十分有意

义的．

国内其他组对格点QCD的研究，综述见文献

[20]．胶球新算符的构造和相应的质量谱计算见文

献[21]．

感谢Eric B．Gregory的讨论和帮助．

附录A Metropolis抽样方法的伪码

for(i=1；i< ；i

系统状态为{Sf ；

文件的指针

char pname[MPL MAX_PROCESSOR_NAME]； ／／处理

++){随机从格点中选一格点，并设当前 器名

求此格点四周格点的自旋值；

定义此格点自旋翻转后的系统的状态为{Sf +，；

计算翻转此格点自旋后对系统能量改变的值dH =

日({Sf +，)一日({Sf )，

if(dH<0){／／跃迁几率为1，接受此自旋翻转{else{

／／根据随机数决定是否翻转

if($<exp(一dH／k T))／／$为0到1之间的随机数

自旋翻转

，，否则不翻转

f

{／／end for

附录B lsing模型的MPI主函数C程序

void main(int argc，char*argv[])

{

FILE *rid； ／／存储本进程的所测量物理量的输出

char filename[3]； ／／存储本进程的所测量物理量的

输出文件

int nprocs，namelen，i， ；

MPLInit(&arge，&argv)； ／／MPI初始化

MPL Comm-size(MPL COMM．_WORLD，&nprocs)； ／／取

通讯器中进程的数目

MPI—Comm-rank(MPL COMM～WORLD，&pid)； ／／取

当前的进程编号

MPL Get_processor_name(pname，&namelen)； ／／取当

前处理器的名字

／* 设置处理器(进程)的奇偶性*／

．f(P==1)

myparity=2； ／／进程数为1时

else if(p／d％ 2 1 =0)

myparity=1； ／／进程数为奇时

else

myparity=O； ／／进程数为偶时

1 2 3 4 5 6 7 8 9 0

第2期 刘军等：Ising模型的并行计算 l27

，* 近邻近程 *，

up_id=(pid+1)％ P；

1w_id=(pid+P一1)％ P；

，*各进程读取数据文件中的参数*，

readData()；

，，打开文件，供写入本进程的计算结果

脚=fopen(“dataIn．dat”，⋯’)；／／dataln．dat为输入

数据文件

／／格点自旋组态按均匀随机分布初始化为spin[i][J]=

1或一1

for(i=1；i(=LP；i++){

f0r( =0；J(L；J++){

spin[i][ ]=one_rand()；

}

}

／／设置随机数种子

srand((unsigned)(time((int*)0) pid))；

，，调用模拟程序

simulate(fid)；

MPL Finalize()；

}

附录C lsing模型的MPI子函数simulate的C

程序

void simulate(FILE*ofp){

int i，J；

，，开始时预热系统，扫描整个网格Ⅳ 次(N=mea$．

are—maln)

for(i=1；i(=measure—nlMn；i++){

，，每隔interva1次update才测量一次(以保持组

态的独立性)

for(j=1；J<=interval； ++){

update()；／／在up&re()子函数中将完成对

网格的并行化处理，

}

，，测量，结果记录到文件中

me~ure(ofp)；

}

}

附录D Ising模型的MPI子函数update的C程序

void update(void){

，，发送底行(第1行)给low_id，再接收数据到顶行

cache—row(1，1w~id)；

／／行扫描更新

for(i=LP；i>1；i一一){

，，对当前行执行Metropolis抽样

cache行

IOW

—

update(i)；

}

，，先发送顶行数据到up—id，再接收数据到底端

cache～mw(LP，up—id)；

／／更新最底行

rOW_update(1)；

附录E Ising模型的MPI子函数measure的C

程序

void measure(FILE * ofp)}

int i，J，im， ，jm，

spin；

double magnetization=0．0；

double interaction=0．0；

，

tIIis_spin， south_ spin，easL

／*平均磁场强度*，

，*扫描格点 *￡，测量磁场强度和内能*，

for(i=l；i(=LP；i++){，，每行

im = i一1：

ip= i+1：

for(j=O；J<￡； ++)}，／每列

jm=( 一1+￡)％ L；／／左行

jP=( +1)％ L；／／右行

tIlis_I spin=spin[i][J]；

souttLI spin=spin[／p][ ]；

e舾L spin=spin[i儿 ]；

magnetization+ =tIlis～spin：

interaction+=this_spin*(s0uttLI spin+easL spin)；

}

}

，，取正值

if(magnetization(0){

magnetization* = 一1：

}

magnetization：magnetization／(L*￡)；

interaction=interaction／(￡*￡)；

，*各进程数据归约*，

MPI—Allreduce (&magnetization。&totM， 1， MPL DOUBLE ，

MPI—SUM，MPI—COMM_WORLD)；

MPLAllreduee(&interaction，&totH，1，MPL DOUBLE ，MPL

SUM，MPLCOMM—WORLD )；

，*主进程记录数据到文件sO*，

if(pid==0){

fprinf(ofp，“％ If％ If＼n”，totM，totH)；

fllush(ofp)；

128 高能物理与核物理(HEP＆NP) 第28卷

Parallel Computing of the Ising M odel

LIU Jun SHEN Yang LUO Xiang Qian

(Department of Physics，Zhongshan University，GuangzhOU 510275，China)

Abstract We give an introduction to the design an d coding of parallel computing for Monte Carlo simulations of lattice

systems，by taking the Ising model as an example．The performance resuhs on our PC cluster are also provided．We be—

lieve that such information is useful for large scale simulation of lattice QCD．

Key words Monte Cado simulations，Ising mod el，lattice gauge theory，parallel algorithm ，hi gh pe rforman ce parallel

computers

Received 26 May 2003，Revised 27 July 2003

*Supported by NSFC(10235040)，Gu．angdong Ministry of Education，Foundation of Zhongshan University Advanced Research Center

1)Corresponding author，E—mall：stslxq@zsu．edu．ca