段紅英

摘 要：雖然Fortran常用來(lái)進(jìn)行科學(xué)計(jì)算，但是面對(duì)計(jì)算量大的程序仍然很耗時(shí)。通常人們用MPI進(jìn)行粗粒度的并行來(lái)提高程序的運(yùn)行效率，近年來(lái)隨著GPU計(jì)算能力的提高，將程序進(jìn)行細(xì)粒度GPU并行化成為一種趨勢(shì)。文章基于NVIDIA公司的CUDA框架，就Fortran程序向CUDA移植過(guò)程中的一些問(wèn)題進(jìn)行總結(jié)，并給出了相應(yīng)的解決方案。

關(guān)鍵詞：Fortran ；C；CUDA；Fortran程序

中圖分類號(hào)：TP311 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：2095-1302（2015）11-00-02

0 引 言

Fortran是常用的科學(xué)計(jì)算語(yǔ)言，其突出的特性就是能實(shí)現(xiàn)自然描述且描述接近數(shù)學(xué)公式，有較好的執(zhí)行效率，但是由于在計(jì)算流體力學(xué)、現(xiàn)代醫(yī)學(xué)影像、分子動(dòng)力學(xué)等領(lǐng)域的模擬中，存在大量的程序計(jì)算，仍然耗時(shí)很多，有的計(jì)算需要幾天甚至幾十天才能完成。為了提高程序的計(jì)算效率，我們將Fortran代碼并行化。通常人們用MPI進(jìn)行粗粒度的并行來(lái)提高程序的運(yùn)行效率，近年來(lái)隨著GPU計(jì)算能力的提高，將程序進(jìn)行細(xì)粒度GPU并行化成為一種趨勢(shì)。

CUDA是NVIDIA公司推出的一種用于 GPU 高性能計(jì)算的軟硬件架構(gòu)，它是對(duì)C語(yǔ)言的擴(kuò)展。在其編程模型中，CPU作為主機(jī)（Host）端調(diào)度整個(gè)程序，GPU作為計(jì)算設(shè)備（device）對(duì)計(jì)算量大、數(shù)據(jù)并行性強(qiáng)的程序并行處理。運(yùn)行在GPU上的并行計(jì)算程序稱為kernel，其必須通過(guò)__global__函數(shù)類型限定符定義，由host端程序調(diào)用啟動(dòng)。

Fortran程序CUDA并行化的完成一般分為編碼、編譯、測(cè)試、優(yōu)化幾個(gè)階段，以下從這幾個(gè)方面分別進(jìn)行總結(jié)闡述。

1 編 碼

在編碼階段，F(xiàn)ortran程序CUDA并行化即為Fortran→C→CUDA的過(guò)程。

1.1 Fortran→C的轉(zhuǎn)化過(guò)程

從Fortran到C的轉(zhuǎn)化過(guò)程，只需要在掌握二者語(yǔ)法的基礎(chǔ)上，逐行翻譯即可，但翻譯工作中需要注意以下細(xì)節(jié)。

（1）數(shù)組

C語(yǔ)言中數(shù)組的起始編號(hào)為0，而Fortran的默認(rèn)起始編號(hào)為1，但也可以用（idx1：idx2）的方式自己定義，這就需要我們把程序中的每個(gè)數(shù)組變量的定義弄明白，翻譯時(shí)對(duì)默認(rèn)定義的數(shù)組標(biāo)號(hào)減1，非默認(rèn)定義的，則用[i-idx1]來(lái)計(jì)算實(shí)際標(biāo)號(hào)。

其次是多維陣列。雖然C和Fortran中所謂的多維陣列都是一個(gè)連續(xù)的一維存儲(chǔ)空間，但是它們對(duì)于行列的分割卻相反。如圖1（a）和圖1（b）分別給出了C和Fortran對(duì)于數(shù)組a[3][2]各自的數(shù)組分割方式。由此，我們?cè)诜g過(guò)程中定義和使用多維數(shù)組時(shí)都須將數(shù)組的行列轉(zhuǎn)換。如a（3：2）變?yōu)镃時(shí)應(yīng)為a[2][3]，對(duì)應(yīng)圖1（b）。

（a） C語(yǔ)言數(shù)組分割圖 （b）Fortran數(shù)組分割圖

圖1 不同數(shù)組的分割圖

（2）函數(shù)參數(shù)傳遞

Fortran中函數(shù)調(diào)用時(shí)一般傳遞的是參數(shù)的內(nèi)存地址，而C既可以進(jìn)行值傳遞，也可以進(jìn)行地址傳遞，一般需要返回多個(gè)參數(shù)值時(shí)用地址傳遞。在翻譯中，為了方便，所有函數(shù)都采用地址傳遞。

（3）函數(shù)重載

在Fortran中為了共享數(shù)據(jù)的方便一般會(huì)用common，如下例所示，F(xiàn)ortran代碼為：

Integer：：I，j，k

Integer：：kk（3）

COMMON/test/I，j，k

用C表達(dá)，需要用以下的方式：

int I，j，k；

int kk[3]；

int*test[3]={&I，&J，&K}；

即把變量的地址連續(xù)地存儲(chǔ)到一個(gè)數(shù)組中。在函數(shù)參數(shù)傳遞時(shí)，在Fortran中調(diào)用addkt函數(shù)就可以傳遞數(shù)組kk，也可以傳遞test，代碼如下：

subroutineaddkt（ kd， kt ）

Integer， intent（inout） ：： kd（3）

但是在C語(yǔ)言中要傳遞以上兩種參數(shù)就出現(xiàn)函數(shù)重載問(wèn)題（一個(gè)為一維數(shù)組，另一個(gè)為二維數(shù)組），但對(duì)于過(guò)程化語(yǔ)言C則沒(méi)有該功能，我們只能把函數(shù)addt定義為addkt1（int *kd，intkt）和addkt2（int *kd[]，intkt）兩個(gè)函數(shù)。

1.2 C→CUDA的轉(zhuǎn)化過(guò)程

CUDA是一種數(shù)據(jù)并行性而非功能并行性的并行計(jì)算解決方案。在C到CUDA的轉(zhuǎn)化過(guò)程中，最關(guān)鍵的就是分析整個(gè)程序，找到最耗時(shí)的代碼部分，分析整個(gè)的可并行性，在對(duì)整個(gè)物理過(guò)程理解的基礎(chǔ)上，進(jìn)行算法設(shè)計(jì)，然后并行化。

以核物理中的蒙卡輸運(yùn)程序?yàn)槔?，蒙特卡羅（MC）方法采用隨機(jī)方法模擬物理過(guò)程，應(yīng)用數(shù)理統(tǒng)計(jì)獲取計(jì)算結(jié)果的計(jì)算方法。蒙卡的整個(gè)輸運(yùn)框架如圖2所示，其中，n為粒子編號(hào)，N為粒子總數(shù)。由于每個(gè)粒子輸運(yùn)過(guò)程相對(duì)獨(dú)立、粒子間通信量少、循環(huán)次數(shù)多，因此，可以一個(gè)粒子對(duì)應(yīng)一個(gè)線程來(lái)并行。

圖2 蒙卡輸運(yùn)框架

當(dāng)然，若是有對(duì)程序足夠的理解，并且Fortran和C都很精通的情況下，則可以直接將Fortran程序CUDA并行化。

2 測(cè) 試

我們借助GDB調(diào)試工具，將測(cè)試過(guò)程分為由下到上，和由上到下兩步。首先，由下到上的對(duì)單個(gè)程序逐級(jí)測(cè)試；然后，根據(jù)程序?qū)懗龆喾N輸入?yún)?shù)，由上往下整體測(cè)試。

判定程序正確的方法就是此程序有相同的輸入和輸出，有隨機(jī)數(shù)的程序會(huì)給我們的測(cè)試帶來(lái)很大的困難。如上面提到的蒙卡輸運(yùn)程序，我們既要保證程序中輸運(yùn)過(guò)程的隨機(jī)性，又要通過(guò)測(cè)試保證程序的正確性。一般大家會(huì)想到產(chǎn)生一個(gè)很大的隨機(jī)數(shù)文件，分別讀入到Fortran和C程序中，此方法可行，但是在粒子數(shù)很多的情況下，效率很低。文章就此問(wèn)題提出了很好的解決方案，此處不再詳細(xì)解述。

2.1 Fortan→C測(cè)試

對(duì)于變量少的程序，我們只需要手動(dòng)打印出需要檢測(cè)的變量值進(jìn)行對(duì)比，但是對(duì)于有幾千個(gè)全局變量的計(jì)算程序，逐一手動(dòng)輸出效率會(huì)很低。因此我們首先要找到程序中用到的全局變量，然后根據(jù)這些變量書寫main函數(shù)測(cè)試。這個(gè)問(wèn)題我們可以借助Perl、Shell等腳本語(yǔ)言完成。其步驟如下：

（1）人工找到程序中所有的全局變量，其類型，維數(shù)，每一維的長(zhǎng)度按照某種格式寫到文件A中；

（2）人工找出程序中的所有程序、子程序及函數(shù)的名字，寫到B中；

（3）書寫腳本，根據(jù)文件A，在所要測(cè)試的Fortran程序的初始化部分，打印出所有全局變量的值，作為C程序的輸入；在所要測(cè)試的Fortran程序結(jié)束之前，打印出所有全局變量的值，作為判斷C程序的標(biāo)準(zhǔn)；

（4）書寫腳本，實(shí)現(xiàn)初始化函數(shù)，即給C語(yǔ)言的變量初始化；

（5）書寫腳本，實(shí)現(xiàn)讀入Fortran的輸出值，判斷C程序的正確性。

在具體的腳本實(shí)現(xiàn)過(guò)程中，需要對(duì)Fortran語(yǔ)法詳細(xì)分解，如一個(gè)程序可能有多處結(jié)束，而在每個(gè)結(jié)束前都需要輸出打印。為了方便，同時(shí)打印出變量的類型、維數(shù)及每一維的長(zhǎng)度。

2.2 Fortan→CUDA測(cè)試

CUDA程序的測(cè)試并不像C程序的測(cè)試那樣簡(jiǎn)單，因?yàn)橛捎谒惴ㄔ?，在功能不變的情況下整個(gè)程序可能會(huì)被打亂，甚至對(duì)變量數(shù)組進(jìn)行排序，由于原來(lái)的數(shù)組會(huì)打亂，致使不能用GDB調(diào)試工具進(jìn)行對(duì)比，并且無(wú)法定位錯(cuò)誤或者判斷計(jì)算中的正確性。雖然程序及數(shù)組順序亂了，但是程序的功能不會(huì)變，并且在相同功能的地方這些數(shù)組之和不會(huì)變，所以我們可以在功能相同的地方對(duì)數(shù)據(jù)之和進(jìn)行對(duì)比。

3 結(jié) 語(yǔ)

面對(duì)科學(xué)計(jì)算中串行程序的瓶頸，我們需要應(yīng)用并行化方法來(lái)解決，本文就Fortran向基于CUDA架構(gòu)的GPU移植過(guò)程中所遇到的一些問(wèn)題進(jìn)行總結(jié)。當(dāng)然為了更高的計(jì)算效率，可以對(duì)初步的程序優(yōu)化。

參考文獻(xiàn)

[1]龔春葉.面向異構(gòu)體系結(jié)構(gòu)的粒子輸運(yùn)并行算法研究[D].長(zhǎng)沙：國(guó)防科學(xué)技術(shù)大學(xué)，2011.

[2] NVIDIA.CUDA C PROGRAMMING GUIDE[Z].2012.

[3] Stephen J.Chapman. Fortran 95/2003程序設(shè)計(jì)[M].劉瑾，龐巖梅，趙越，等，譯.北京：中國(guó)電力出版社，2009.

[4] Jason Sanders，EdwardKandrot.GPU高性能編程CUDA實(shí)戰(zhàn)[M]. 聶雪軍，等，譯.北京：機(jī)械工業(yè)出版社，2011：43-68.

[5] Randal L.Scbwartz，TomPboenix brain d foy. Perl語(yǔ)言入門[M].盛春，蔣永清，王暉，譯.南京：東南大學(xué)出版社，2007.