李 波，沈文楓，鄭衍衡，徐煒民

（上海大學(xué) 計(jì)算機(jī)工程和科學(xué)學(xué)院，上海200072）

0 引 言

心電仿真計(jì)算需要高性能計(jì)算環(huán)境，心電計(jì)算機(jī)仿真計(jì)算是心臟生物電現(xiàn)象研究中的一個(gè)重要課題［1］，隨著計(jì)算機(jī)計(jì)算能力的提高，心電活動(dòng)的計(jì)算機(jī)仿真也隨著迅速發(fā)展。但是對(duì)于心臟模型的仿真計(jì)算，要求在空間上達(dá)到高分辨率、在計(jì)算上達(dá)到高精度，需要巨大的計(jì)算量，必須采用并行計(jì)算［2］。到目前為止，絕大部分的這類計(jì)算是在高性能計(jì)算機(jī)系統(tǒng)中采用（message passing interface，MPI）。GPU用于高性能計(jì)算，使得在高校和研究機(jī)構(gòu)內(nèi)部大量的PC機(jī)在計(jì)算性能上有很大的提高。利用這些PC機(jī)的空閑時(shí)間進(jìn)行計(jì)算，可以代替昂貴的高性能計(jì)算機(jī)，廉價(jià)的、快速的實(shí)現(xiàn)心電仿真計(jì)算。

本文提出采用由PC機(jī)中的CPU和GPU組成異構(gòu)網(wǎng)絡(luò)計(jì)算的方案，利用空閑的PC資源，實(shí)現(xiàn)心電仿真計(jì)算，并且提出了面向短的響應(yīng)時(shí)間要求的調(diào)度算法，取得了很好的效果。

1 相關(guān)工作

Tusscher等進(jìn)行心室傳導(dǎo)系統(tǒng)的模擬模型，采用C＋＋和MPI，是在由10臺(tái)Dell650Precision服務(wù)器，每臺(tái)服務(wù)器用雙核Intel Xeon 2.66GHz CPU組成的Beowolf集群式系統(tǒng)上運(yùn)行的［3］。Reumann et al。在配置有512、1024、2048、4096個(gè)節(jié)點(diǎn)的IBM Blue Gene／L上，進(jìn)行了大規(guī)模的心臟模擬［4］，這兩篇文章表明心電仿真主要在昂貴的高性能計(jì)算機(jī)系統(tǒng)上進(jìn)行，在器件上都沒有采用GPU。本文研究在由CPU和GPU組成的網(wǎng)絡(luò)計(jì)算系統(tǒng)上運(yùn)行心電仿真。Wenfeng Shen et al.研究了在多核和GPU上進(jìn)行心電仿真的并行計(jì)算［1，5］，是在多核CPU和GPU組成的單機(jī)系統(tǒng)中的并行處理。本文研究的是CPU和GPU在多機(jī)系統(tǒng)中的并行處理。

在本文所采用的網(wǎng)絡(luò)計(jì)算系統(tǒng)的體系結(jié)構(gòu)上，和著名的Bonic系統(tǒng)是相近的，但是有兩個(gè)不同點(diǎn)：①本文所研究的系統(tǒng)對(duì)于響應(yīng)時(shí)間方面的要求比較高，這是和Bonic系統(tǒng)的主要差別。Bonic系統(tǒng)的調(diào)度程序以高吞吐率為出發(fā)點(diǎn)，追求公平性和高效率，在Bonic系統(tǒng)中，如果引入短的響應(yīng)時(shí)間的要求，會(huì)影響公平性和高效率［6］。我們研究的心電仿真計(jì)算的要求是較短的響應(yīng)時(shí)間，因此對(duì)于調(diào)度算法提出了不同的要求。②計(jì)算規(guī)模比較小，主要應(yīng)用在高校和研究單位內(nèi)部，對(duì)于提交的作業(yè)采取先提交先服務(wù)的策略，在運(yùn)行時(shí)，獨(dú)占全部可用資源，要求每一個(gè)作業(yè)的運(yùn)行時(shí)間盡量短，可以減少作業(yè)的平均等待時(shí)間，相當(dāng)于提高了系統(tǒng)的計(jì)算速度［7－8］。同時(shí)，計(jì)算資源來自部門內(nèi)部，可以獲得各個(gè)處理器的計(jì)算性能，為設(shè)計(jì)面向較短的響應(yīng)時(shí)間的調(diào)度程序提供了條件。

2 心電仿真計(jì)算原理及其算法

心臟的計(jì)算機(jī)模型一般是通過對(duì)心臟斷面圖像離散化并進(jìn)行空間插補(bǔ)，本論文所采用實(shí)驗(yàn)?zāi)Ｐ蜑槿毡疚菏先呐K模型［9－10］，該模型是由一個(gè)56＊56＊90的斜交坐標(biāo)系構(gòu)成。模型上部為心房，下部為心室，模型細(xì)胞在空間上是等距的，直徑為1.5ms。模型的實(shí)際細(xì)胞數(shù)據(jù)為5萬(wàn)個(gè)，通過使模型細(xì)胞具有方向性和具有和實(shí)際細(xì)胞相近的電生理特性，在解剖學(xué)和電生理學(xué)兩個(gè)方面對(duì)模型作進(jìn)一步的修飾使得模型更加接近實(shí)際的心臟。

2.1 心電仿真算法

心電計(jì)算機(jī)仿真的核心便是心臟興奮序列的計(jì)算機(jī)仿真，所用算法是其關(guān)鍵所在，直接決定了計(jì)算的速度和仿真的效果，并影響到模型的實(shí)用性和仿真結(jié)果。在心電興奮傳播的過程中，算法的目標(biāo)是在每一個(gè)計(jì)算過程中確定未興奮的心肌細(xì)胞時(shí)候接受臨近的已興奮心肌細(xì)胞的興奮傳播而開始其自身的興奮過程。

心電計(jì)算的仿真算法主要分為三步：第一步是以惠更斯原理模型為基礎(chǔ)，仿真心電興奮傳播的過程。第二步是根據(jù)傳播過程中記錄的偶極子數(shù)目，計(jì)算心電圖的容積導(dǎo)體和體表電勢(shì)。第三步是保存計(jì)算結(jié)果。其中第二步計(jì)算容積導(dǎo)體和體表電勢(shì)占用了整個(gè)計(jì)算機(jī)仿真過程的絕大部分時(shí)間。在體表軀干表面，心內(nèi)膜和心外膜等四邊邊界大約3800三角形元素的電勢(shì)計(jì)算，產(chǎn)生大概50000個(gè)電偶極子源。容積導(dǎo)體中每個(gè)三角元素是異勢(shì)的。根據(jù)特定方程，從當(dāng)前偶極子源采用邊界元方法計(jì)算的體表電位，心外膜和心內(nèi)膜。如圖1所示。

電勢(shì)計(jì)算方程如下

圖1 人體軀干

式中：（rk，j）——在第k個(gè)表面上的第j個(gè)三角的體表電勢(shì)值，rk，j——中央位置向量。－號(hào)和＋號(hào)分別代表內(nèi)部和外部信號(hào)。M則代表細(xì)胞模型的當(dāng)前偶極子數(shù)目，J代表偶極子源。在整個(gè)容積導(dǎo)體中體表、體外膜表面、右心內(nèi)膜表面和左心內(nèi)膜表面分別被劃分為344，1002，307，278個(gè)節(jié)點(diǎn)和684，2002，610，552個(gè)三角形。心電仿真計(jì)算主要是對(duì)方程（1）進(jìn)行計(jì)算機(jī)分解計(jì)算，求得計(jì)算結(jié)果。

2.2 心電仿真計(jì)算流程圖

心電仿真計(jì)算流程如圖2所示。

圖2 心電仿真計(jì)算流程

通過實(shí)驗(yàn)驗(yàn)證，對(duì)算法中每一步計(jì)算作時(shí)間分析，得到如下結(jié)果：在第一步的傳播過程中，所占時(shí)間為整個(gè)計(jì)算機(jī)仿真過程的0.1%，第二步的電勢(shì)計(jì)算占總仿真時(shí)間的99.7%，保存結(jié)果大概占總時(shí)間的0.2%。因此主要從第二步計(jì)算電勢(shì)來考慮，尋求最優(yōu)并行化。第二步是根據(jù)傳播過程中記錄的偶極子數(shù)目，計(jì)算心電圖的容積導(dǎo)體和體表電位。由于心電計(jì)算的特殊性質(zhì)，在每3ms的時(shí)間周期內(nèi)所計(jì)算的偶極子數(shù)目不盡相同。每個(gè)周期的電勢(shì)計(jì)算只和這個(gè)周期產(chǎn)生的偶極子數(shù)目相關(guān)，一般情況下產(chǎn)生的偶極子數(shù)目越多，電勢(shì)計(jì)算所花費(fèi)時(shí)間就越多。而偶極子數(shù)目在第一階段的興奮性傳播過程中便可以統(tǒng)計(jì)出來，通過合理的分析可以預(yù)測(cè)計(jì)算負(fù)載，為調(diào)度算法創(chuàng)造了條件。

3 異構(gòu)網(wǎng)絡(luò)計(jì)算的結(jié)構(gòu)

3.1 心電仿真計(jì)算在計(jì)算結(jié)構(gòu)上的特點(diǎn)分析

心電仿真計(jì)算在計(jì)算結(jié)構(gòu)上有以下的特點(diǎn)：高度并行，占總仿真時(shí)間99.7%的計(jì)算是高度并行的，計(jì)算模塊之間的相互通信量非常少，各個(gè)模塊完全可以獨(dú)立計(jì)算。非常適合利用網(wǎng)絡(luò)計(jì)算。各個(gè)模塊的計(jì)算量是可以預(yù)測(cè)的。為在網(wǎng)絡(luò)計(jì)算中實(shí)現(xiàn)定量調(diào)度提供了條件。

3.2 異構(gòu)網(wǎng)絡(luò)結(jié)構(gòu)介紹

異構(gòu)網(wǎng)絡(luò)結(jié)構(gòu)如圖3所示。

圖3 異構(gòu)網(wǎng)絡(luò)結(jié)構(gòu)

整個(gè)網(wǎng)絡(luò)由多臺(tái)的四核處理器、八核處理器、GPU處理器組成。這個(gè)網(wǎng)絡(luò)環(huán)境在高校的實(shí)驗(yàn)室和研究所的研究室還是有比較大的代表性。

4 調(diào)度算法

調(diào)度算法的目標(biāo)是達(dá)到按照各個(gè)服務(wù)器處理速度的比例分配任務(wù)，達(dá)到這個(gè)目標(biāo)的方案很多，最簡(jiǎn)單的方案是依照處理器的計(jì)算能力按比例進(jìn)行分配，但是GPU用作高性能計(jì)算時(shí)，它的性能是隨著被處理的程序的并行程度變化的，因此在調(diào)度程序中要盡量使得計(jì)算量大（并行度高）的模塊分配給GPU計(jì)算。圖4是Geforce440GPU的處理速度隨著處理模塊的偶極子數(shù)量變化曲線。

調(diào)度算法的執(zhí)行步驟如下：

圖4 Geforce440GPU的處理速度與處理模塊偶極子數(shù)量的關(guān)系

設(shè)不帶有GPU的處理器有m臺(tái)，帶有GPU的處理器有n臺(tái)。

（1）把總計(jì)算量按照CPU和GPU的計(jì)算能力分二大類，以保證異構(gòu)系統(tǒng)中GPU的能力得到充分發(fā)揮。然后，在分類內(nèi)部按照處理器能力分配計(jì)算量。步驟如下：

1）求出需要計(jì)算偶極子的總量，用Q總計(jì)算量表示。

2）對(duì)于各個(gè)CPU和GPU種類，用實(shí)驗(yàn)方法獲得每類的計(jì)算能力（單位時(shí)間內(nèi)可以計(jì)算的偶極子數(shù)量），CPUi計(jì)算能力和GPUi計(jì)算能力 表示分別表示第i個(gè)CPU和GPU的計(jì)算能力。按照降序排序，i＝1，表示計(jì)算能力最強(qiáng)。

3）把異構(gòu)系統(tǒng)中的計(jì)算機(jī)處理能力，按照CPU和GPU分為兩類，對(duì)兩類處理器能力分別進(jìn)行統(tǒng)計(jì)，CPU總計(jì)算能力和GPU總計(jì)算能力分別代表CPU處理器和GPU處理器總的計(jì)算能力。

4）把計(jì)算總量為Q總計(jì)算量分配給CPU處理器和GPU處理器，用CPU總計(jì)算量和GPU總計(jì)算量表示CPU處理器和GPU處理器所承擔(dān)的總計(jì)算量。則GPU總計(jì)算量：CPU總計(jì)算量＝GPU總計(jì)算能力：CPU總計(jì)算能力，GPU總計(jì)算量＋CPU總計(jì)算量＝Q總計(jì)算量

5）按照計(jì)算能力的大小計(jì)算每個(gè)處理器所分配的計(jì)算量，用CPUi計(jì)算量表示分配給第i類CPUi的計(jì)算量，i＝1，2…m。GPUi計(jì)算量表示分配給第i類GPUi的計(jì)算量，i＝1，2…n。則CPUi計(jì)算量＝CPU總計(jì)算量＊CPUi計(jì)算能力／CPU總計(jì)算能力。同樣，GPUi計(jì)算量＝GPU總計(jì)算量＊GPUi計(jì)算能力／GPU總計(jì)算量。

6）對(duì)偶極子計(jì)算任務(wù)按照計(jì)算量進(jìn)行降序排序，并且劃分成Q1，Q2…QN，其中N＝n（GPU處理器數(shù)）＋m（CPU處理器數(shù)）。使得 Qi＝GPUi計(jì)算量，i＝1，2…n，Qn＋I(xiàn)＝CPUi計(jì)算量，I＝1，2…m。

（2）對(duì)于GPU1處理器，按照偶極子計(jì)算量的降序逐個(gè)模塊提交，對(duì)于CPUm處理器，按照偶極子計(jì)算量的升序逐個(gè)模塊提交。對(duì)于其它的GPUi處理器（i＝2，3…n）和CPUi處理器（i＝1，2…m－1）按照偶極子計(jì)算量的中間優(yōu)先，逐步向二邊交叉延伸。這樣的提交次序，便于在計(jì)算的最后階段，把剩余計(jì)算的任務(wù)向鄰居節(jié)點(diǎn)轉(zhuǎn)移。

（3）網(wǎng)絡(luò)計(jì)算的特點(diǎn)之一，是在運(yùn)行過程中會(huì)發(fā)生增加或者減少處理器的情況，此時(shí)，按照以下步驟進(jìn)行：

1）計(jì)算出還沒有提交的任務(wù)（包括在減少的處理機(jī)中被中斷并且取消的任務(wù)）的剩余總計(jì)算量，設(shè)定為Q剩余總計(jì)算量。并且按照新的系統(tǒng)配置情況計(jì)算出到全部任務(wù)執(zhí)行完成還需要多少時(shí)間，設(shè)定為T剩余總計(jì)算時(shí)間。

2）計(jì)算出正在正常運(yùn)行的各個(gè)處理器i當(dāng)前的任務(wù)結(jié)束還需要多少時(shí)間，設(shè)定為T i剩余計(jì)算時(shí)間，修改處理器i的計(jì)算能力為CPUi新計(jì)算能力＝CPUi計(jì)算能力×（1－T i剩余計(jì)算時(shí)間／T剩余總計(jì)算時(shí)間），對(duì)于GPU處理器也同樣處理。

3）按照新的系統(tǒng)配置情況，使用調(diào)度程序計(jì)算各個(gè)處理器需要承擔(dān)的任務(wù)。由于各個(gè)處理器的剩余計(jì)算量是不同的，因此各個(gè)處理器的新計(jì)算能力也是各不相同的，原有的分類將取消，每個(gè)處理器都單獨(dú)成為一個(gè)類。

（4）當(dāng)某一處理器k結(jié)束自身計(jì)算任務(wù)，同時(shí)k－1或者k＋1處理器還有多個(gè)任務(wù)等待，可以把相應(yīng)的計(jì)算任務(wù)遷移到處理器k上。直至整個(gè)系統(tǒng)計(jì)算任務(wù)結(jié)束為止。

（5）在開始計(jì)算過程以后，需要密切注意GPUn計(jì)算速度，如果發(fā)現(xiàn)GPUn速度低于期望值很多時(shí)，則表示在GPUn中，計(jì)算了大量的小任務(wù)，使得GPU計(jì)算速度下降。從我們分析偶極子計(jì)算量的分布情況來看，70－80%的偶極子模塊計(jì)算時(shí)，GPU處理器的計(jì)算速度基本上保持著最佳狀態(tài)，也就是說，只要在網(wǎng)絡(luò)計(jì)算系統(tǒng)中，CPU處理器的計(jì)算量超過20%，我們所提出的調(diào)度算法是可以獲得最佳結(jié)果的。我們認(rèn)為在高性能計(jì)算中，這種情況是有典型性的。當(dāng)然如果在極端情況下，計(jì)算任務(wù)中包含了大量的小任務(wù)，迫使GPU處理器承擔(dān)了很多小任務(wù)，由于我們的調(diào)度算法在分配GPU處理器和CPU處理器承擔(dān)的計(jì)算量時(shí)，是按照GPU處理器在處理大任務(wù)情況下的性能。當(dāng)GPU處理器需要處理大量小任務(wù)時(shí)，完成所分配的任務(wù)需要的時(shí)間將會(huì)超過預(yù)定時(shí)間，適當(dāng)把分配給GPU處理器的任務(wù)遷移到CPU處理器。在遷移的過程中，可以參照步驟（3）的方法。

5 實(shí)驗(yàn)結(jié)果

實(shí)驗(yàn)環(huán)境：

CPU處理器：Intel（R）Core（TM）i7CPU 920＠2.67GHz雙核。

GPU處理器：NVIDIA GeForce 9800GT。

實(shí)驗(yàn)內(nèi)容：

（1）采用3種不同配置，如下：

配置1：CPU處理器10臺(tái)，GPU處理器30臺(tái)。

配置2：CPU處理器20臺(tái)，GPU處理器20臺(tái)。

配置3：CPU處理器30臺(tái)，GPU處理器10臺(tái)。

（2）在每個(gè)配置上，采用二種不同算法，如下：

算法1：本調(diào)度算法。

算法2：先到先服務(wù)算法，有多個(gè)空閑處理器時(shí)，選用性能最高的處理器。

（3）輸入負(fù)荷：室上性心動(dòng)過速心電仿真的4.8s偶極子數(shù)。

實(shí)驗(yàn)結(jié)果：

如圖5所示，對(duì)于三證不同的配置平臺(tái)，分別使用算法1和算法2。

圖5 不同配置下的實(shí)驗(yàn)結(jié)果

6 結(jié)束語(yǔ)

本論文在分析了心電仿真計(jì)算特點(diǎn)的基礎(chǔ)上，提出了在網(wǎng)絡(luò)計(jì)算系統(tǒng)中的一個(gè)新的調(diào)度算法，在系統(tǒng)的響應(yīng)時(shí)間方面取得了很好的效果。

在本文心電仿真計(jì)算中，除了在并行性方面適合網(wǎng)絡(luò)計(jì)算的要求以外，其中的主要前提是計(jì)算工作量的可預(yù)測(cè)性，這個(gè)要求在高性能計(jì)算方面還是有一定的可行性的，例如，本調(diào)度算法可以應(yīng)用到N－Body問題Fmm算法中，因?yàn)樵诿總€(gè)分塊中近距離粒子作用力的計(jì)算量和粒子數(shù)的平方成正比，它的計(jì)算量也是可以預(yù)測(cè)的。

進(jìn)一步工作，將在總結(jié)其它類型的高性能計(jì)算的算法特點(diǎn)的基礎(chǔ)上，研究動(dòng)態(tài)調(diào)度算法。

［1］Shen W，Wei D，Xu W，et al.Parallelized computation for computer simulation of electrocardiograms using personal computers with multi－core CPU and general－purpose GPU［J］.Computer Methods and Programs in Biomedicine，2010，100（1）：87－96.

［2］McFarlane R，Biktasheva I V.High performance computing for the simulation of cardiac electrophysiology［C］／／Sliema：Proceedings of the Third International Conference on Software Engineering Advances，2008：13－18.

［3］Tusscher K，Panfilov A.Modelling of the ventricular conduction system［J］.Progress in Biophysics and Molecular Biology，2008，96（1－3）：152－170.

［4］Reumann M，F(xiàn)itch BG，Rayshubskiy A，et al.Large scale cardiac modeling on the Blue Gene supercomputer［C］／／Vancouver，BC：Proceedings of IEEE Engineering in Medicine and Biology Society，2008：577－580.

［5］Reumann M，F(xiàn)itch B G，Rayshubskiy A，et al.Strong scaling and speedup to 16，384processors in cardiac electro－mechanical simulations［C］／／Minneapolis，MN：Annual International Conference of the IEEE Engineering in Medicine and Biology Society，2009：2795－2798.

［6］Donassolo B，Legrand A，Geyer C.Non－cooperative scheduling considered harmful in collaborative volunteer computing environments cluster，cloud and grid computing（CCGrid）［C］／／Newport Beach，CA：11th IEEE／ACM International Symposium on，2011：144－153.

［7］Sato D，Xie Y，Weiss J N，et al.Acceleration of cardiac tissue simulation with graphic processing units［J］.Medical and Biological Engineering and Computing，2009，47（9）：1011－1015.

［8］Zhu H，Sun Y，Rajagopal G，et al.Facilitating arrhythmia simulation：The method of quantitative cellular automata modeling and parallel running［OL］.Biomedical Engineering Online，2004－03－29.

［9］Zhu Xin，Wei Daming，Wang Hui.Simulation of intracardial potentials with anisotropic computer heart models［C］／／Shanghai：Proceedings of the IEEE Engineering in Medicine and Biology，2005：1071－1074.

［10］Zhu Xin，Wei Daming.Computer simulation of intracardiac potential with whole－h(huán)eart model［J］.International Journal of Bioinformatics Research and Applications，2007，3（1）：100－122.