丁磊，王武，姜金榮，趙蓮

1. 中國科學院計算機網絡信息中心，北京 100190

2. 中國科學院大學，北京 100049

引言

N 體問題[1]（N-body problem）是已知多個粒子的初始位置、速度和質量，求解出其在經典力學下的后續(xù)運動的問題，屬于高性能數值計算的七類典型應用之一[2]，在宇宙學模擬[3]、分子動力學研究[4]等諸多領域有十分重要的應用。

常用的N 體求解方法有直接法[5]（Particle-Particle, PP）、粒子-網格法[5]（Particle-Mesh, PM）、樹方法[6]（Barnes-Hut, BH）、快速多極子方法[7]（Fast Multipole Method, FMM）等。FMM 采用長程力近似計算方式，基于樹結構的多層遞歸盒子劃分策略，使用核心函數的級數展開，將遠程力的作用近似為粒子與盒子的作用及盒子間的作用，極大降低了場中粒子相互作用的計算復雜度。由于可以指定劃分的層數和展開階數，FMM 能有效平衡執(zhí)行時間與計算誤差[8]。與BH 相比，FMM 基于高階展開，因此精度可以達到更高，同時可以控制計算量。

近年來，國內外有許多FMM 相關研究與應用成果。如曹小林[9]等基于JASMIN 框架實現了FMM 的求解器；Cruz[10]等基于PETSc 框架實現了PetFMM 軟件；Winkel[11]等基于MPI 實現了N 體計算軟件PEPC。隨著異構系統的不斷發(fā)展，學者開始針對特定體系結構進行研究。Lashuk[12]等使用256塊GPU 在NCSA Lincoln cluster 上對KIFMM 進行加速；Yokota[13]等在TSUBAME 2.0 系統上使用4096塊GPU 來計算各向同性湍流；王武[14]等基于申威眾核處理器進行FMM 優(yōu)化。此外，還有基于FMMPM 方法在GPU 上進行N 體模擬實現[15]等。

Charm++是一種基于消息驅動的編程框架，具有過分解（Over-decomposition）、可遷移（Migratability）、異步消息驅動（Asynchronous Message-Driven Execution）等特點，可以提供運行時的動態(tài)負載均衡與容錯機制[16]。近年來，有許多工作利用Charm++的動態(tài)負載均衡特性，取得了不錯的效果。如中國科學院計算機網絡信息中心研發(fā)的并行框架SC_Tangram[17]、基于Charm++的非結構網格[18]等。本文中，我們將基于Charm++實現FMM 程序的并行化，并與MPI 版本的并行實現進行對比分析，結果表明整體性能有10%左右的提升。

1 并行FMM 介紹

1.1 FMM 的計算流程

FMM 是一種區(qū)域劃分的計算方法，其核心思路在于對位勢函數在遠場進行多極展開，再將多極展開轉化為近場的局部展開，繼而通過位勢計算出粒子在引力場所受到的作用力。

FMM 使用球坐標的多極展開來近似多個粒子形成的引力場對外界的作用，用局部展開來近似外界引力場對區(qū)域內粒子的作用，包括P2M、M2M、M2L、L2L、L2P、P2P 六個計算步驟。

下列公式中，α、β、ρ分別表示球坐標的極角、方位角與徑向距離，mi為粒子質量。P為展開階數，球諧展開中的系數m、n與P正相關。

（1）P2M (Particle to Multipole Expansion)

P2M 操作只發(fā)生在計算區(qū)域的最底層，也就是最小的盒子區(qū)域內。主要的思想是根據盒子里面的粒子位置，計算出這個盒子的多極展開系數。

P2M 的計算方法如下[19]：

其中，k 為盒子中的粒子數，Mnm是多極展開系數，Ynm(θ,φ)是球諧函數，計算公式為：

式中Pn(x)為Legendre 多項式。

（2）M2M (Multipole Expansion to Multipole Expansion)

M2M 操作的計算方法如下[19]：

其中，Mkj是轉移中心之后的多極展開系數，Okj是下一層盒子的多極展開系數。

M2M 是一個自下向上的過程，不斷地收集下層的多極展開系數，計算本層的多極展開系數，再到更高一層進行計算。

（3）M2L (Multipole Expansion to Local Expansion)

P2M 操作與M2M 只進行多極展開系數的計算，目的是計算出盒子對空間的作用效果。而M2L 操作是計算相互作用。

M2L 是一個多極展開轉化為局部展開的過程，計算方法如下[19]：

其中，Lkj為局部展開系數，Onm為多極展開系數。

在計算M2L 的過程中，需要找到當前盒子的次相鄰盒子，即與當前盒子不相鄰，但是父盒子與當前盒子的父盒子相鄰的盒子。根據這些盒子的多極展開系數，計算出局部展開系數。

（4）L2L (Local Expansion to Local Expansion)

L2L 的計算公式如下[19]：

其中，Lkj為局部展開系數，Onm為局部展開系數。

L2L 是一個與M2M相反的過程，完成了相互作用的中心轉移。在M2L 中只考慮了本層的相互作用，那些來自父盒子的作用力是通過L2L 來計算。這樣，就可以計算出來自非鄰居盒子的作用。

（5）L2P (Local Expansion to Particle)

L2P 是P2M 的逆過程。經過L2L 操作后，已經在最底層盒子中獲取來自非鄰居盒子的作用，而L2P 就是把這些局部展開系數轉為盒子中粒子所受到的作用力。

L2P 的計算公式如下[19]：

（6）P2P (Particle to Particle)

由于近似計算需要兩個區(qū)域有一定的距離，無法使用近似的方式求解出最底層盒子來自鄰居盒子的作用力。因此，只能依照萬有引力公式，求解出鄰居盒子中粒子的作用力。

1.2 并行FMM 算法

FMM 的并行分解即為空間分解。如算法1 所示，每個并行單元負責一個區(qū)域，在區(qū)域中按照預定分布初始化全局粒子的位置，通過partition 操作將粒子劃分給處理區(qū)域，并依照希爾伯特編碼建立樹。此后執(zhí)行粒子到盒子以及盒子到盒子的P2M 和M2M 操作，獲得多極展開系數。

算法1 并行FMM

鑒于P2P 操作和M2L 操作需要了解相鄰和次相鄰盒子的信息，需要進行一次全局通信，將這些信息保存在本地關建樹（local essential tree）中。最終執(zhí)行樹的自上而下的遍歷，完成L2L 與L2P 操作。因此，在并行FMM 中需要進行兩次全局通信，分別發(fā)生在partition 和local essential tree 處。

2 Charm++介紹

Charm++是基于C++的面向對象并行編程框架，由UIUC 并行編程實驗室開發(fā)，具有過分解、可遷移、異步消息驅動的特點[20]。過分解指并行計算對象的數量多于實際的處理器數量?？蛇w移指并行計算對象可以在處理器間移動，完成動態(tài)負載均衡。異步消息驅動指每個處理器上的并行對象的執(zhí)行流由消息驅動，計算對象間的消息傳遞采用異步方式，不阻塞執(zhí)行流。

Chare 是Charm++的基本計算單元，其本質是分布在不同的處理器上的C++對象。代理（Proxy）是一種特殊的C++對象，用以表示一個遠程的Chare。Chare 間的通信通過調用代理的特殊成員方法來實現，由Runtime 負責將本次調用的參數封裝成消息，并發(fā)送到對應的遠程Chare，該過程通常被稱為遠程調用。

Charm++提供了Chare array、group、node group等并行數據結構。其中，Chare array 提供一個統一的代理，每個array 的成員Chare 都有一個thisIndex，可以通過proxy[someIndex]來調用編號為someIndex的Chare 上的方法。

PUP (Pack and Unpack)框架是Charm++提供的序列化框架，負責遠程調用和Chare 遷移時的序列化工作。

PE (Processing Element)和Node 是Charm++的抽象概念。其中PE 是Chare 的執(zhí)行實體，其上有許多Chare 和消息。Node 由多個PE 組成，是一種共享內存的抽象表達。Node 上的不同Chare 共享Node的地址空間。在實現上，PE 一般對應線程，而Node對應于進程。

Charm++的各種實體概念如圖1所示。

Charm++的項目由ci (Charm++ interface)文件和C++源碼構成。其中ci 文件定義了Charm++中的Chare、message、module 等實體。預處理系統解析ci 文件后生成對應頭文件和實現文件，使用者需要在C++源碼中引入這些文件。

圖1 Charm++ 概念示意圖Fig.1 Charm++ concept diagram

Charm++不支持全局可變變量，但支持全局只讀變量。只讀變量在ci 文件中使用readonly 關鍵字引入，在mainChare 的構造函數中進行初始化。

每個PE 上都有一個消息隊列。以PE 的視角看，有一個執(zhí)行流不斷從隊列中取出消息。如果消息是種子，將創(chuàng)建新的Chare。如果消息是遠程調用，則找到PE 上對應的Chare，將消息解包并調用。在調用過程中可能會產生新的種子和遠程調用，由Runtime 決定新消息對應的PE。

因此，Chare 間的消息傳遞是異步的方式。調用者只知道消息在隊列里，但不清楚消息是否已被執(zhí)行，即調用者無法立即獲得執(zhí)行的結果。使用者可以定義一個遠程方法供接收方調用，并在這個遠程方法中獲取通信結果，執(zhí)行后續(xù)操作。這種風格使代碼雜亂且易錯。為了方便表示執(zhí)行流的依賴關系，Charm++提供了SDAG[21](structure dagger)表示方法。

SDAG 的語法采用 { } 模式，其中有when、serial 等。serial 用于表示一段完整執(zhí)行的不被搶占的代碼。when 表示等待消息，如 when method (arg){serial {do_something (arg)}} 表示當前執(zhí)行流等待遠程方法method 被調用，并在該條件觸發(fā)后執(zhí)行do_something 操作。

Charm++預處理器會分析ci 文件中的SDAG 語法，并生成對應的代碼，將SDAG 語句塊拆分成許多continuation 函數，即多個階段。以when 語句為例，Runtime 會在method 被觸發(fā)后，調用下一個階段對應的continuation。這樣使用者無需自行編寫等待、判斷以及跳轉的遠程方法。

3 并行實現

3.1 任務分解與預處理

在本文中，我們定義一個名為master 的mainChare作為并行任務的發(fā)起者和規(guī)約消息的接收者，定義一個名為slave 的一維Chare array，作為基本的并行計算單元。

鑒于Chare 需要在PE 上執(zhí)行，而且存在運行時遷移的情況，我們將原始的全局變量保存在slave 的成員中，并以參數的方式傳遞進去。

對于計算開始前已知且只讀的全局變量，使用Charm++提供的readonly 變量表示，并在master 中進行初始化。

3.2 異步通信轉化

Charm++采用異步通信的方式，這使得通信的發(fā)起方無法了解到消息是否被接收。

在MPI 中，我們可以在一個函數中發(fā)起多次通信。對于阻塞通信，操作系統會將函數的實參和局部變量保存在棧中，等到通信結束后再進行調度；對于非阻塞通信，程序中有循環(huán)代碼對通信結束條件進行判斷，控制流始終在該函數中。在Charm++中，用戶只能通過消息驅動來觸發(fā)后續(xù)操作，這意味著通信發(fā)起時，我們就失去了局部變量和實參的內存分配。為了恢復狀態(tài)，需要將跨越通信的實參和局部變量保存到Chare 的成員中。

針對通信函數XXmethod，我們引入結構體XXmethod_param 來保存這些信息?？紤]到實參的初始化有傳值和傳引用的方式，我們只在結構體里保存變量地址。對于傳值的方式，需要引入init_XX_by_value 來動態(tài)分配內存并保存；對于傳引用的方式，引入init_XX_by_ref，來保存已有變量的值，并提供一個release 方法用于釋放值初始化時分配的空間。

另一方面，某個通信函數可能被調用多次，結束后跳轉到不同的執(zhí)行流。為了保證通信結束后函數的正常驅動，我們需要保存通信結束后的下一步執(zhí)行的位置。

本文中，我們將需要通信的函數設置為slave的成員函數，并根據通信將程序劃分為多個state。在slave 的成員變量中，使用一個棧來保存通信結束后下一步執(zhí)行的函數指針以及state 值。通過set_continuation 來將函數指針和state 壓棧，do_continuation 來將函數指針和state 出棧，并調用該函數。在每個state 開始時，我們取出函數中需要用到的局部變量和實參。在第一個state 的時候對局部變量進行初始化。實參的初始化發(fā)生在函數的調用處。

圖2 MPI 風格示例通信函數Fig.2 Example of communication in MPI style

圖2是一個MPI 風格的通信函數示例，在Some-Method 函數中，有全局通信global_communication。其中有實參parameter，有局部變量local_variable，這些變量都要跨越通信的。

圖3是對應的Charm++表達方式。全局通信將SomeMethod 分成了兩個部分。state 是Chare array 的成員變量，用于標識當前處理的狀態(tài)。CharmSome-Method_param 是我們針對通信函數CharmSome-Method 引入的結構體，其中保存了該函數執(zhí)行時需要的實參parameter、局部變量local_variable。

在state 為0 的時候，我們調用set_continuation 函數來指定在本次通信后需要執(zhí)行的函數為Charm-SomeMethod 的下一個state，并調用init_XX_by_value來為局部變量分配內存。實參parameter 的初始化發(fā)生在CharmSomeMethod 被調用處。此后，我們通過CharmSomeMethod_param 獲得parameter 和local_variable，并執(zhí)行串行操作do_some_serial_job_1。鑒于通信函數global communication 需要傳入local_variable 的引用，這里我們對global_communication_param 使用init_XX_by_ref 進行初始化，記錄local_variable 的地址。在執(zhí)行完global_communication 通信后，執(zhí)行流到了CharmSomeMethod 的state1，再一次從CharmSomeMethod_param 中獲得parameter 和local_variable 變量，執(zhí)行串行操作do_some_serial_job_2。由于該函數已經執(zhí)行到最后，可以釋放CharmSome-Method_param 中的空間，并調用do_continuation，獲取之前保存的下一步執(zhí)行的函數和state。global_communication 內部實現與CharmSomeMethod 類似。

圖3 Charm++風格示例通信函數Fig.3 Example of communication in Charm++ style

本文只對最基本的通信函數使用SDAG。這樣對于調用層數較深處發(fā)生的通信，只需要提供基本通信的實現，無需將所有包含通信的函數都改寫成SDAG 風格，可以較大程度減少代碼修改量。至此，除基本通信函數外，所有的MPI 風格的同步通信邏輯都可以轉化為Charm++異步通信的表示，且用戶修改部分較少。

基本通信函數的實現將在下一節(jié)中介紹。

3.3 基本通信的實現

并行FMM 的通信發(fā)生在partition 和local essential tree 操作中，包括基本通信Alltoallv 與Allreduce。

Alltoallv 需要每個并行計算單元向其它并行計算單元發(fā)送不等長的消息，同時從其它并行計算單元中接收這些不等長的消息。Allreduce 需要對所有并行單元的某個參數進行規(guī)約。

圖4 message 類型的定義Fig.4 Defination of message

在Charm++中，使用message 比參數序列化快[22]，而且用戶無需過多了解序列化知識。我們針對不同的類型，定義了對應的message 類型，如圖4所示。其中，data 保存信息的本體，size 表示消息的大小，from 表示消息的發(fā)送者。

偽代碼 1 Alltoallv 的 Charm++實現

serial {sendbuf, sendcount, sdispl = get_data_from (param)for i=0 to Chare_size do message = init_message (sendbuf, sendcount,sdispl)thisProxy[i].send_message (message)end for recv_cnt = Chare_size}while (recv_cnt --){when send_message (msg)serial{recvbuf, recvcount, rdispl = get_data_from (param)set_variable_from (recvbuf, recvcount, rdispl, msg)}}serial {do_continuation ()}

在實現Alltoallv 時，我們使用結構體來保存整個通信過程需要用到的變量。在serial 語句塊中，每個Chare 構造消息并發(fā)送，同時觸發(fā)等待消息的when 語句塊，從消息中取出數據，放到緩沖區(qū)中。Alltoallv 的實現如偽代碼1 所示，其中param 保存了alltoalv 調用時的參數，thisProxy 指Chare array 的代理，通過這個代理向自己發(fā)送消息。

偽代碼2 Allreduce 在slave 上的實現

Allreduce 的實現包括規(guī)約和廣播。首先在slave上獲取到所有的規(guī)約數據，指定master 的處理函數來做規(guī)約。在規(guī)約結束后，master 發(fā)布廣播消息，觸發(fā)slave 的when 語句，master 和slave 的實現過程見偽代碼2 與偽代碼3。

偽代碼3 Allreduce 在master 上的實現

3.4 負載均衡與遷移策略

大規(guī)模分布式并行系統上動態(tài)負載均衡，可以分成分布式、集中式和層次式三種[23]。

分布式策略只考慮相鄰處理器，如Cybenko[24]等給出了擴散算法的一般化形式，處理器從高負載處理器中獲取負載，并把負載給低負載處理器。Hui[25]等提出了基于水動力學的負載均衡，利用水的連通性特性來控制負載均衡。

集中式策略考慮全局負載，將信息集中到一個處理器上進行計算。如全局隨機指派[26]，貪心方法[26]等。

層次式負載均衡策略根據拓撲結構建立層次樹。如HBM[27]根據網絡結構將處理器劃分成多個組，在層次結構的每一層都執(zhí)行負載均衡遷移，直到根節(jié)點為止。

Charm++提供了在運行時移動Chare 的方法，使得處理器上的負載可以進行動態(tài)調節(jié)。本文中，我們利用MigrateMe 函數，使用集中式負載均衡策略，在運行時對slave 的Chare 進行遷移。

圖5是負載均衡的實現時序圖。slave 在某個階段將自己包含的粒子信息發(fā)送到master 上。master負責對這些信息進行分析，經由負載均衡策略給出一個Chare 的重新分布，并依次通知相應的Chare。這些Chare 調用MigrateMe 函數來移動到新的PE 上，并向master 匯報遷移的結束。master 收到所有的遷移結束消息后，向slave 發(fā)起一個廣播，使其繼續(xù)執(zhí)行計算。

圖5 負載均衡時序圖Fig.5 Sequence Diagram of load balance

本文中，負載均衡的執(zhí)行時期為partition 結束時，此時所有的Chare 都已經拿到自己負責處理的粒子。為了簡化負載均衡的實現邏輯，我們假設所有的處理器都是相同性能，僅考慮PE 的計算負載，并用PE 上的總粒子數來近似。

該負載均衡問題被稱為相同并行機調度（Identical Parallel-Machine Scheduling）問題，屬于NP-hard 問題[28]。本文中，我們采用最長處理時優(yōu)先策略（Longest Processing Time，LPT）來求解，簡述如下：

(1)將Chare 按照粒子數量從大到小排序

(2)依次將Chare 指派到當前粒子總數最少的PE 上。

該策略可以達到4/3 近似[28]，即對于該問題的所有情況，近似策略都能給出最優(yōu)解4/3 倍以內的可行解。

3.5 實現技術

上述message、Alltoallv、Allreduce 等需要針對不同類型實現。此外，對于每個通信函數，都需實現一個XX_param 的結構體，其成員函數需要覆蓋到在通信中使用到的變量，且應當區(qū)分by_value 和by_ref。在通信函數的最開始部分，都需要獲得全部的XX_param 結構體成員。這些代碼都是相似的，可以使用腳本語言生成。

本文將需要生成的類型參數寫成json 的格式，通過Python 腳本解析，基于模板自動生成了message、基礎通信函數、XX_param 結構體的C++代碼片段和宏定義。在并行FMM 源碼中只需在對應的地方引入這些宏定義即可，無需大量修改源碼。

主要的實現流程包括：

(1)全局變量局部化；

(2)基于通信劃分state；

(3)Python 腳本生成C++代碼片段；

(4)調整通信函數內部風格；

(5)在通信函數中引入XX_param 等代碼。

4 數值結果

以下測試均在中國科學院計算機網絡信息中心超級計算機“元”上進行測試。每臺刀片計算節(jié)點配置 2 顆 Intel E5-2680 V2（Ivy Bridge | 10C | 2.8GHz）處理器，64 GB DDR3 ECC 1866MHz 內存。編譯環(huán)境為Intel Composer XE 2013 編譯器，Charm++ 版本為6.6.1，編譯時使用-O3 優(yōu)化選項和C++0x 標準。

4.1 正確性檢驗

為了驗證FMM 的Charm++實現正確性，我們選擇了不同分布與參數進行測試，比較了Charm++與MPI 下每個粒子的受力值，結果完全一致，這表明我們的計算精度和已有的MPI 實現的精度相同。

4.2 性能檢驗

以16384 萬個粒子為例，采用均勻的立方體分布作為粒子的初值分布，分別對MPI 和Charm++實現進行測試，保證Chare 數量和PE 數量相等（即不使用過分解），結果如圖6所示?？梢钥闯觯S著處理器規(guī)模從32 核增加到1024 核，兩種實現的速度基本接近，Charm++實現相比MPI 實現在千核規(guī)模上有提升。

圖6 16384 萬粒子的執(zhí)行時間Fig.6 Execution time of 163,840,000 particles

在上述測試中，以32 核規(guī)模的執(zhí)行時間作為基準，計算了Charm++和MPI 版本的加速比，結果如圖7所示?？梢钥闯?，與理想情況相比，FMM 的Charm++實現在256 核以下的強擴展性能較好。隨著核數增加到千核規(guī)模，全局通信會擴大，擴展性有一些下降。與MPI 實現相比，Charm++實現在1024 核規(guī)模的加速效果較好，但在64 核至512 核規(guī)模的加速效果較為遜色。這一現象是因為千核規(guī)模下，通信占比增大，而MPI 實現使用全局通信，在partition 階段的通信時間遠大于使用點對點通信的Charm++實現。

圖7 并行 FMM 的加速比Fig.7 Speedup ratio of different implementations

為了驗證過分解對負載不均衡情況的優(yōu)化效果，本文針對65536 萬粒子的球殼狀不均勻分布，固定處理器數量為64，調整每個處理器上的Chare 數量，從64 增加到512，進行了測試。我們統計了每個處理器上的總粒子數Ni，并計算了所有處理器上的極差(Nmax-Nmin)和標準差

表1 過分解測試結果Table 1 Result of over-decomposition

測試結果如表1所示，在不使用過分解的情況下，可以看出粒子分布的標準差和極差較大。隨著Chare 數量增加，LPT 負載均衡策略降低了標準差和極差的數量級，也減少了計算時間。最優(yōu)情況在使用256 個Chare 時取得，與不采用過分解的情況相比，減少了10%的計算時間。上述結果表明，過分解和負載均衡策略對粒子分布不均勻下的FMM 計算有一定優(yōu)化效果，可以平衡不同處理器上的計算負載。

5 總結與展望

本文基于Charm++并行編程框架，實現了并行FMM 程序。首先分析了FMM 的并行實現關鍵和主要通信，介紹了Charm++并行編程環(huán)境。在并行實現的過程中，本文討論了MPI 編程模型和Charm++編程模型的異同，引入了用異步模型表示同步計算的一般性方法，給出了一個通用的MPI 并行程序向Charm++轉化的流程，并針對基本通信函數給出了SDAG 的實現方案。最終，本文在“元”超級計算機上對FMM 的Charm++實現進行了檢驗與測試。結果表明，基于Charm++實現的FMM 具有較好的擴展性，針對負載不均衡的情況，過分解和LPT 負載均衡策略可以減少10%的執(zhí)行時間。

文本中實現中只利用了Charm++提供的Chare array 類型，即把每個Chare 看作一個獨立且內存不共享的個體。如何使用Charm++提供的邏輯節(jié)點與本地存儲來減少通信量，提高并行效率，將是未來的一個嘗試方向。另一方面，本文在負載均衡時僅考慮了PE 上的總粒子數，且假定計算負載與粒子數成正相關，不一定反映了真實的負載情況。因此，使用更貼近FMM 的負載均衡指標也是一個可以考慮的方向。

針對目前高性能計算系統廣泛使用異構計算的情況[29]，本文提供的Charm++實現也具有一定適應性。只需聲明一組綁定PE 的Charm group，在其中執(zhí)行耗時的P2P 與M2L 操作。異構加速與具體體系結構有關，涉及內容較多，故本文僅提供接口。

利益沖突聲明

所有作者聲明不存在利益沖突關系