llo 標(biāo)量勢(shì)差分計(jì)算的Co－array Fortran 解決方案及算例

姜 義，劉 瑩

(長(zhǎng)春工程學(xué)院 理學(xué)院，長(zhǎng)春130012)

0 引言

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

1 并行處理和Co－Array fortran 語(yǔ)言［1］

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

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

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

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

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

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

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

Co－array Fortran 中的兩個(gè)重要函數(shù)This_image()和Num_Images()分別表達(dá)進(jìn)程副本序號(hào)和進(jìn)程副本總數(shù)。其應(yīng)用如下例:

該代碼段將把其它各進(jìn)程副本中Co－array y 累加存儲(chǔ)于當(dāng)前副本中的Co－array x。

進(jìn)程副本之間的同步通常直接調(diào)用sync all 即可。需要注意sync all 必須保證在全部進(jìn)程副本中不被越過(guò)，否則將導(dǎo)致進(jìn)程副本之間的不同步。以下是個(gè)極端的例子

顯然進(jìn)程副本1 中的Sync all 無(wú)效，從而無(wú)法與其它副本同步。因此sync all 操作不應(yīng)在進(jìn)程副本之間存在差異的選擇、轉(zhuǎn)移和循環(huán)塊中調(diào)用，必須確保其不被繞開(kāi)?；蛘哒f(shuō)各進(jìn)程副本中的Sync all 調(diào)用必須數(shù)目相同，因?yàn)槊恳淮蜸ync all 調(diào)用都意味著各進(jìn)程之間某個(gè)流程節(jié)點(diǎn)處“對(duì)一次鐘表”，以便共同出發(fā)，進(jìn)行后續(xù)計(jì)算!

2 靜電勢(shì)差分計(jì)算模型的Co－Array fortran 處理方案

2.1 模型描述

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

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

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

圖1 差分計(jì)算原始模型

圖2 任務(wù)拆分和Co－array

2.2 任務(wù)拆分

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

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

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

Image1 的左邊界格點(diǎn)刷新公式為

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

右邊界格點(diǎn)

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

Image2 的左邊界格點(diǎn)刷新公式為

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

右邊界格點(diǎn)

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

每刷新格點(diǎn)勢(shì)函數(shù)一遍，將邊界格點(diǎn)數(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］

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

2.3 代碼及解釋

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

3 算例結(jié)果和討論

3.1 并行效率

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

表1 不同任務(wù)拆解下的并行效率對(duì)比

3.2 物理模型討論

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

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

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

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

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

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

［3］ 王家禮，朱滿(mǎn)座，等.電磁場(chǎng)與電磁波［M］.西安:西安電子科技大學(xué)出版社，2012:87－88.