llo 標量勢差分計算的Co－array Fortran 解決方案及算例

姜 義，劉 瑩

(長春工程學院 理學院，長春130012)

0 引言

隨著計算技術的迅猛發(fā)展，幾乎所有工程技術應用領域都已進入定量化數(shù)值處理階段。近年來，64 位多處理器計算機的廣泛普及給存在定制方案/模型需求的非專業(yè)編程工作人員帶來了極大的便利，它提供了自主開發(fā)專業(yè)適應性更強的小型(微型)計算軟件最為基本的強力硬件支持。但是，由于專業(yè)、行業(yè)限制，如何簡便快捷地發(fā)揮先進計算硬件的哪怕是基本功能，也是多數(shù)應用工程技術人員難以面對，甚至是無法面對的。有鑒于此，本文詳細地介紹線性差分計算的并行處理方案并且給出詳盡的處理手續(xù)。

1 并行處理和Co－Array fortran 語言［1］

并行處理技術涉及廣域網(wǎng)、局域網(wǎng)服務器/工作站節(jié)點之間和單機(服務器或工作站)的多CPU、多核協(xié)同運算處理，以實現(xiàn)耗時冗長的大運算量任務。不論在何種硬件架構下，應用程序的并發(fā)執(zhí)行都無法回避進程/線程之間的數(shù)據(jù)互訪和同步問題。有鑒于此，專業(yè)軟件研發(fā)機構開發(fā)了多種并行計算接口或者環(huán)境，為非專業(yè)人士自主開發(fā)小型并行程序提供支持，其中較為流行的是OPENMP(由OpenMP Architecture Review Board 提出的并行處理方案)和MPI(Multi Point Interface)［2］。但是這兩種并行處理接口/環(huán)境在高級語言編程中仍然涉及頻繁的函數(shù)調用，這對于不具備軟件專業(yè)工作背景的人士而言仍然存在障礙。由于歷史沿襲和積淀，科學計算編程多采用FORTRAN 語言，一種淺顯易懂、便于操作的并行編程語言擴展應運而生，這就是Co－Array Fortran。

Co－Array Fortran 用于處理需要通過消息傳遞(message－passing)模型來處理的任務分解和并發(fā)處理模型，類似于Fortran，并且已經(jīng)融入Fortran 編程語言體系，使得在Fortran 程序編制過程中可以簡潔、直觀地進行并行處理。有了Co－Array 方案，只需具備Fortran 語言編程基礎，掌握有限幾條新的語法規(guī)則，就可以開發(fā)出穩(wěn)定、高效的并行程序。

從應用角度來看，本文不擬全面介紹Co－Array Fortran，只關注并行處理的三個關鍵步驟:任務分解，進程副本(Image)之間的互訪和同步問題。其中的任務分解將在下節(jié)通過具體模型加以詳解，本節(jié)只討論進程鏡像之間的數(shù)據(jù)互訪和同步。

為了實現(xiàn)數(shù)據(jù)互訪，Co－Array Fortran 提供了一條可以在Image 之間交換數(shù)據(jù)的變量/數(shù)組聲明語法，例如:

real，dimension(10)，codimension［*］::x，y，z

其含義是每一進程副本都定義名稱為x，y，z 三個實型、10 元素的一維數(shù)組，該數(shù)組可以被任意副本(Image)訪問，*表示遍及全部進程副本，x，y，z 即為Co－array。如果編譯環(huán)境中設定開出3 個進程副本，則每一進程副本進行本地運算時，可以直接操作變量符號x 進行數(shù)組變量運算，也可以操作或運算該數(shù)組變量的第k 個元素x(k)。一旦涉及進程副本之間的操作，則必須指明被操作的是哪個副本的哪個元素，x(2)［2］表示第2 個進程副本中Co－array x 的第二個元素，y(5)［3］表示第3 個進程副本中Co－array y 的第5 個元素，依此類推。

進程副本內的Co－array 操作一般無需給出Co－array 的歸屬，而進程副本間的Co－array 操作必須指明其歸屬才有效。例如當前進程副本為副本1，則x=y(:)［2］將把進程副本2 中Co－array y 的各元素依次賦給進程1 的Co－array x。

Co－array Fortran 中的兩個重要函數(shù)This_image()和Num_Images()分別表達進程副本序號和進程副本總數(shù)。其應用如下例:

該代碼段將把其它各進程副本中Co－array y 累加存儲于當前副本中的Co－array x。

進程副本之間的同步通常直接調用sync all 即可。需要注意sync all 必須保證在全部進程副本中不被越過，否則將導致進程副本之間的不同步。以下是個極端的例子

顯然進程副本1 中的Sync all 無效，從而無法與其它副本同步。因此sync all 操作不應在進程副本之間存在差異的選擇、轉移和循環(huán)塊中調用，必須確保其不被繞開。或者說各進程副本中的Sync all 調用必須數(shù)目相同，因為每一次Sync all 調用都意味著各進程之間某個流程節(jié)點處“對一次鐘表”，以便共同出發(fā)，進行后續(xù)計算!

2 靜電勢差分計算模型的Co－Array fortran 處理方案

2.1 模型描述

依據(jù)唯一性定理，給定邊界條件下拉普拉斯方程的解是唯一的［3］。但是涉及多區(qū)域或者多對象、復雜邊界條件的靜電模型一般無法解析處理，需要借助數(shù)值計算才能對空間電勢、電場分布狀況進行準確估量。其中有限線性差分法是明智的選擇，這種用算術手段數(shù)值求解拉普拉斯方程的方法對于數(shù)值處理能力和硬件高度發(fā)達的今天是方便易行的。

算例為無限長圓柱狀平行相離雙導體模型，其一處于單位電勢，另一圓柱體懸浮。顯然，模型本質上是二維的，二維空間的有限線性差分法適用，并且模型可以抽象為圖1。圖中上下邊界接地，左右兩側采用周期性邊界條件，這樣做的理由是距離高電位導體充分遠處電場強度幾乎可以忽略，或者說對電工設備的安全性影響極小(這點在計算結果中可見一斑)。

計算模型網(wǎng)格橫縱向節(jié)點均為2N，理由是作者使用的服務器配置是雙路且每路4 核心，可以進行2 ～8 進程副本的有效并行，充分發(fā)揮服務器CPU 資源顯然8 副本是最為恰當?shù)摹Ｒ话愕姆掌?工作站的總核心數(shù)都是2N個，小節(jié)3.3 中的代碼可以直接使用。

圖1 差分計算原始模型

圖2 任務拆分和Co－array

2.2 任務拆分

下面以二拆分(2 進程副本協(xié)同處理模式)為例說明任務拆解規(guī)則。假定網(wǎng)格分布為512*512 格點陣，則每一進程副本需完成的差分格點空間為256*512。本算例采用4 最近鄰平均法刷新中心格點勢函數(shù)量值，所以進程副本image(1)格點空間左邊界格點的勢函數(shù)刷新必須使用image(2)的右邊界格點數(shù)據(jù)，相應地，image(2)的右邊界格點勢刷新也需要使用image(1)的左邊界格點數(shù)據(jù);image(1)的右邊界則與image(2)的左邊界相關聯(lián)。也就是說，每完成一次格點勢刷新，必須在兩個進程副本之間進行左右邊界格點數(shù)據(jù)的相互交換，該交換數(shù)據(jù)量為512 元素的列向量，顯然該向量就是一對Co－array 數(shù)組，如圖2 所示，且

任意進程副本的(i，j)內格點勢函數(shù)值的刷新公式為(注:因上下邊界接地，無需刷新，故而1＜j＜512)

φ(i，j)=(φ(i－1，j)+φ(i+1，j)+φ(i，j－1)+φ(i，j+1)/4

Image1 的左邊界格點刷新公式為

φ(1，j)=(x(j)+φ(2，j)+φ(1，j－1)+φ(1，j+1)/4

右邊界格點

φ(256，j)=(φ(255，j)+y(j)+φ(256，j－1)+φ(256，j+1)/4

Image2 的左邊界格點刷新公式為

φ(1，j)=(x(j)+φ(2，j)+φ(1，j－1)+φ(1，j+1)/4

右邊界格點

φ(256，j)=(φ(255，j)+y(j)+φ(256，j－1)+φ(256，j+1)/4

每刷新格點勢函數(shù)一遍，將邊界格點數(shù)據(jù)集中互傳一次

x=φ(1，j)，y=φ(256，j)

z=x，x(:)［1］=y(:)［2］，y(:)［2］=z(:)［1］

x(:)［2］=y(:)(1)，y(:)［1］=z(:)［2］

調用sync all 實現(xiàn)副本之間的同步，然后進行下一次格點刷新，直至計算結果滿足精度要求。

2.3 代碼及解釋

下屬代碼為單機運行通用程序，由于版面限制，省略了其中的差分循環(huán)、誤差估計和輸出等常規(guī)處理模塊。變量聲明模塊中的maxGridx、maxGridy、maxStep 和localGridx 分別為模型橫向格點總數(shù)、縱向格點總數(shù)、最大差分循環(huán)步數(shù)和每副本內的橫向格點數(shù)。

3 算例結果和討論

3.1 并行效率

在給定硬件資源的前提下，影響并行效率的因素仍然很多，最主要的是大任務拆解的合理性和硬件資源的利用率。本算例運行的硬件環(huán)境為聯(lián)想萬全T350 G7 64 位服務器，雙CPU，每CPU 有4 核心，核心主頻為2.4GHz，服務器內存32G。分別在512/1024/2048 差分網(wǎng)格密度下針對前述計算模型進行了2/4/8 image的并行測試，測試結果用機時消耗表達且以兩image 作為參考，列于下表。在任務很小的情形下，如果拆分過細，會導致進程、線程或者進程副本之間的數(shù)據(jù)傳遞消耗過多機時，從而降低效率。表中數(shù)據(jù)表明，512 差分網(wǎng)格密度4 拆分(4 Images 并行處理)已經(jīng)足夠，進一步拆分效率提升并不顯著。而2048 網(wǎng)格密度所對應的計算量很大，所以從4 核心運算(4 拆分)到用滿8 核心(8 拆分)效率提升仍然十分顯著。

表1 不同任務拆解下的并行效率對比

3.2 物理模型討論

圖3 勢函數(shù)分布圖

由于累積誤差的影響，格點之差分刷新方向不僅影響勢函數(shù)的計算結果，而且會影響收斂進程(需要的總差分循環(huán)步數(shù))。為了消除格點刷新方向性的影響，實際計算過程中，每一差分循環(huán)包含自左向右和自右向左一對操作，或者說對應一對相向的差分循環(huán)進程。實踐證明，這樣做可以有效抑制格點刷新方向性對收斂進度和收斂品質的影響，否則會出現(xiàn)循環(huán)步數(shù)不足導致的懸浮導體表面不等勢。圖3 為40000 步差分循環(huán)的計算結果，其中圖3(a)是稀疏等位線密度下的數(shù)值標注圖式，便于了解勢函數(shù)的定量分布，而圖3(b)則為加密等位線圖式，更有助于定性考察電場強度之空間分布。圖式結果清楚表明，懸浮導體的引入使空間等勢線受到擠壓，形成一個等勢區(qū)域。高電位導體表面電場強度遠大于懸浮導體表面場強，如果環(huán)境介質均勻，則可能性的介質電離應始于高電位導體表面附近，可能的介質擊穿通道一定在高電位導體和懸浮導體之間。

圖4 計算誤差統(tǒng)計分布

差分計算雖然原理簡單，操作方便，但是其中的誤差和傳遞不可回避。圖4 為40000 差分循環(huán)末尾100步逐格點標準誤差計算結果分布。從圖中可見，懸浮導體附近標準誤差要較其它區(qū)域高一個數(shù)量級。究其原因，應當來源于靜電感應過程和懸浮導體電荷遷移進程。也就是說，這里的差分循環(huán)迭代過程不僅僅是求解拉普拉斯方程，還是懸浮導體在外電場中的靜電感應過程模擬(導體的靜電極化!)。盡管實際靜電感應過程在10－9秒的“瞬間”內就可以完成，數(shù)值模擬過程中的電荷移動必須逐步推進，而且越是接近靜電平衡，該進程越發(fā)緩慢。

［1］ Numrich，R W，Reid J K.Co－arrays in the next Fortran Standard［J］.Scientific Programming，2006(14):1－18.

［2］ 周偉明.多核計算與程序設計［M］.武漢:華中科技大學出版社，2009:75－76.

［3］ 王家禮，朱滿座，等.電磁場與電磁波［M］.西安:西安電子科技大學出版社，2012:87－88.