mpif_2003

Transcription

1 Fortran 語言 MPI 平行計算程式設計 編著 : 鄭守成 期 間 : 民國 92 年 2 月 12 日 電話 : (03) x c00tch00@nchc.gov.tw 1

2 Fortran 語言...1 MPI 平行計算程式設計...1 第一章 前言 MPI 平行計算軟體 國家高速網路與計算中心的平行計算環境 在 IBM 電腦系統上如何使用 MPI IBM 電腦系統的 MPI Fortran 程式編譯指令 IBM 電腦系統的 Job command file IBM 電腦系統的平行程式的執行指令 在 P luster 上如何使用 MPI P luster 上的 MPI Fortran 程式編譯指令 P luster 上的 Job command file P luster 上的平行程式執行指令...15 第二章 無邊界資料交換的平行程式 MPI 基本指令 mpif.h include file MPI_INIT, MPI_FINALIZE MPI_OMM_SIZE, MPI_OMM_RANK MPI_SEND, MPI_REV 無邊界資料交換的循序程式 T2SEQ 資料不切割的平行程式 T2P MPI_SATTER,MPI_GATHER,MPI_REDUE MPI_SATTER,MPI_GATHER MPI_REDUE, MPI_ALLREDUE 資料切割的平行程式 T2DP...32 第三章 需要邊界資料交換的平行程式 MPI_SENDREV, MPI_BAST MPI_SENDREV MPI_BAST 邊界資料交換的循序程式 T3SEQ 資料不切割的邊界資料交換平行程式 T3P 資料切割的邊界資料交換平行程式 ( 一 ) T3DP_ 資料切割的邊界資料交換平行程式 ( 二 ) T3DP_ 第四章 格點數不能整除的平行程式 格點數不能整除的循序程式 T4SEQ MPI_SATTERV MPI_GATHERV MPI_PAK UNPAK BARRIER WTIME

3 4.3.1 MPI_PAK MPI_UNPAK MPI_BARRIER MPI_WTIME 資料切割的平行程式 T4DP...67 第五章多維陣列的平行程式 多維陣列的循序程式 T5SEQ 多維陣列資料不切割的平行程式 T5P 多維陣列末維資料切割的平行程式 T5DP 與二維切割有關的 MPI 副程式 垂直座標圖示法則 (artesian Topology) 界定二維切割的 MPI 副程式 MPI_ART_REATE...93 MPI_ART_OORDS MPI_ART_SHIFT 定義固定間隔資料的 MPI 副程式...96 MPI_TYPE_VETOR MPI_TYPE_OMMIT 多維陣列末二維切割的平行程式 T5_2D...98 第六章 MPI 程式的效率提昇 Nonblocking 資料傳送 資料傳送的合併 以邊界資料計算取代邊界資料交換 輸出入資料的安排 事先切割輸入資料 事後收集切割過的輸出資料 第七章導出的資料類別 導出的資料類別 陣列的轉換 兩方迴歸與管線法 第八章多方依賴及 SOR 解法 四方依賴及 SOR 解法 黑白點間隔 SOR 解法 斑馬線 SOR 解法 八方依賴與四色點間隔 SOR 解法 第九章有限元素法程式 有限元素法的循序程式 有限元素法的平行程式 附錄一撰寫 語言的 MPI 程式 參考書目 Parallel Processing without Partition of 1-D Arrays Parallel Processing with Partition of 1-D Arrays Parallel on the 2 nd Dimension of 2-D Arrays without Partition

4 Parallel on the 2 nd Dimension of 2-D Arrays with Partition Partition on the 3rd dimension of 3-D Arrays

5 第一章前言 本章將介紹 MPI 平行計算軟體 國家高速網路與計算中心現有的平行計算環境 以及在各型機器上如何使用 MPI 第一節簡單介紹 MPI 平行計算軟體 第二節介紹國家高速網路與計算中心現有的平行計算環境 第三節介紹如何在 IBM 電腦系統上使用 MPI, 包括路徑的設定 平行程式的編譯 及平行程式的執行等 第四節介紹如何在 P cluster 上使用 MPI, 包括路徑的設定 平行程式的編譯 及平行程式的執行等 5

6 1.1 MPI 平行計算軟體 MPI (Message Passing Interface) 是第一個標準化的 Message Passing 平行語言 可以使用在 Fortran ++ 等語言撰寫的程式上 MPI 平行程式可以在分散式記憶體平行系統上執行, 也可以在共用記憶體平行系統上執行 目前系統廠商所提供的 MPI 軟體是屬於 MPI1.2 版 它提供了一百多個副程式, 讓程式人員來選用 MPI 協會在 1998 年公布了 MPI 2.0 版的規格, 數年之後就會有 MPI 2.0 版的軟體可用了 日前美國的 Argonne National Lab 已經公布了 MPIH 版的整套軟體, 該版含有 MPI 2.0 版的部份功能 有興趣的讀者可以免費自網路下載該軟體, 其網址是 也可以用 anonymous ftp 下載該軟體, 其網址是 ftp.mcs.anl.gov 其下目錄 (directory) pub/mpi 裏檔名 mpich tar.z 或 還有許多與 MPI 相關的資訊可供參考 mpich tar.gz, 在該目錄之下 6

7 1.2 國家高速網路與計算中心的平行計算環境 目前國家高速網路與計算中心的 IBM SP2 IBM SP2 SMP IBM p690 HP SPP2200 SGI Origin2000 和 Fujitsu VPP300 等系統上均有該公司自備的 MPI 平行軟體,P cluster 上是裝用 MPIH 公用平行軟體, 也都有能力執行平行程式 但是到目前為止, 只有 P cluster IBM SP2 IBM SP2 SMP 和 IBM p690 設有一個 PU 只執行一個程式的平行環境, 其他機器上則無此種設定 例如, 若有一個用戶要用四個 PU 來執行其平行程式, 他在 IBM SP2 上取得四個 PU 之後, 這四個 PU 就僅只執行這個平行程式直到它執行完畢為止, 不會有其他程式進來跟他搶 PU 時間 但他在其他機器 ( 如 HP SPP2000) 上取得四個 PU 之後, 如果所有使用者對 PU 的需求數量超過該系統的 PU 總數時, 他所取得四個 PU 之中的每一個 PU, 都有可能要跟其他程式以分時方式 (time sharing) 共用一個 PU HP SPP2000 和 SGI ORIGIN2000 為共用記憶體平行系統, 這種電腦系統是 16 顆 PU 共用一組記憶體 SP2 和 VPP300 是屬於分散式記憶體平行系統, 每一個 PU 備有它獨用的記憶體 IBM SP2 SMP 及 IBM p690 是共用記憶體及分散式記憶體混合的平行系統,SP2 SMP 每一個 node 備有 4 顆 PU 共用一組記憶體, 目前備有 42 個 node 的 SMP cluster p690 每一個 node 備有 32 顆 PU 共用一組記憶體, 目前備有 8 個 node 的 SMP cluster SP2 SP2 SMP 和 p690 是採用該系統專屬的工作排程軟體 (job scheduler) LoadLeveler 來安排用戶的批次工作 (batch job) 使用者必須備妥 LoadLeveler 的 job command file, 使用 llsubmit 指令把該批次工作交給該系統來執行 SPP2000 ORIGIN2000 和 VPP300 是採用 NQS (Network Queue System) 工作排程軟體來安排用戶的批次工作 使用者必須備妥 NQS 的 job command file, 使用 qsub 指令把該批次工作交給各該系統來執行 P cluster 是採用 DQS (Distributed Queue System) 工作排程軟體來安排用戶的批次工作, 其使用方式類似 NQS 7

8 1.3 在 IBM 電腦系統上如何使用 MPI 首先, shell 用戶要在自己 home directory 的.cshrc 檔裏加入下列路徑, 這樣才能夠抓得到 include file (mpif.h mpif90.h mpi.h) 編譯指令 (mpxlf mpxlf90 mpcc) MPI library 和 LoadLeveler 指令 (llsubmit llq llstatus llcancel) set lpath=(. ~ /usr/lpp/ppe.poe/include /usr/lpp/ppe.poe/lib) set lpath=($lpath /usr/lpp/ppe.poe/bin /home/loadl/bin ) set path=($path $lpath) 加好上述路徑之後, 將.cshrc 存檔, 再執行 source.cshrc 指令, 即可進行平行程式的編譯與執 行 簽退 (logout) 後再簽到 (login) 之後就不必再執行 source.cshrc 指令 IBM 電腦系統的 MPI Fortran 程式編譯指令 使用 MPI 的 Fortran 77 平行程式, 其編譯器 (compiler) 一般叫做 mpif77, 但是在 IBM 電腦系統上卻叫做 mpxlf mpxlf 常用的編譯選項如下 : mpxlf -O3 -qarch=auto -qstrict -o file.x file.f 其中選項 -O3 是作最高級的最佳化 (level 3 Optimization), 可使程式的計算速度加快數倍 -qarch=auto 是通知編譯器該程式要在同型機器上執行 -qstrict 是通知編譯器不要改變計算的順序 -o file.x 是指定執行檔名為 file.x, 不指定時其內定 (default) 檔名為 a.out IBM 電腦系統的 Job command file 要在 IBM SP2(ivy) 上執行平行程式, 使用者必須備妥 LoadLeveler 的 job command file 例如, 下面這個 job command file 叫做 jobp4, 它要在四個 PU 上執行平行程式 file.x 8

9 #!/bin/csh executable = /usr/bin/poe #@ arguments = /your_working_directory/file.x #@ output = outp4 #@ error = outp4 #@ job_type = parallel #@ class = medium #@ min_processors = 4 #@ max_processors = 4 #@ requirements = (Adapter == "hps_user") #@ wall_clock_limit = 20 #@ queue euilib us 其中 executable = /usr/bin/poe 是固定不變,poe 是指 Parallel Operating Environment arguments = 執行檔所在之全路徑及檔名 output = 標準輸出檔名 (stdout) error = 錯誤訊息 (error message) 輸出檔名 class = SP2 PU 的分組別, 使用 llclass 指令可以看到分組別 : short (PU 時間上限為 12 小時, 共有 10 顆 120MHz PU) medium (PU 時間上限為 24 小時, 共有 64 顆 160MHz PU) long (PU 時間上限為 96 小時, 共有 24 顆 120MHz PU) min_processors = 最少的 PU 數目 max_processors = 最多的 PU 數目 requirements = (Adapter == "hps_user") 是固定不變 wall_clock_limit = 該 job 最多需要的時間, 單位為分鐘 queue 是固定不變 平行計算可以使用的 PU 數目,short class 最多 4 個 PU,medium class 最多 32 個 PU, long class 最多 8 個 PU 由於 MPI 1.2 版不具備取得 PU 控制 PU 和歸還 PU 的功能, 所以 min_processors 和 max_processors 要填相同的數字 要在 IBM SP2 SMP (ivory) 上執行平行程式, 使用者必須備妥 LoadLeveler 的 job command file 例如, 下面這個 job command file 叫做 jobp4, 它要在四個 PU 上執行平行程式 file.x 9

10 #!/bin/csh network.mpi= css0,shared,us executable = /usr/bin/poe #@ arguments = /your_working_directory/file.x #@ output = outp4 #@ error = outp4 #@ job_type = parallel #@ class = medium #@ tasks_per_node = 4 #@ node = 1 #@ queue euilib us 由於 IBM SP2 SMP 每個 Node 含有四棵 375MHz PU 共用 4GB 或 8GB 的記憶體 class = SP2 SMP PU 的分組別, 使用 llclass 指令可以看到分組別 : short (PU 時間上限為 12 小時, 3 個 Node 共有 6 顆 PU) medium (PU 時間上限為 24 小時,32 個 Node 共有 128 顆 PU) bigmem (PU 時間上限為 48 小時, 4 個 Node 共有 16 顆 PU) 這個 class 一個 Node 備有 8GB 的共用記憶體 tasks_per_node=4 是說明一個 Node 選用四棵 PU node=1 是說明要用一個 Node, 一共四棵 PU 平行計算可以使用的 PU 數目 medium class 是 8 個 Node 一共 32 顆 PU 要在 IBM p690 上執行平行程式, 使用者必須備妥 LoadLeveler 的 job command file 例如, 下面這個 job command file 叫做 jobp8, 它要在 8 個 PU 上執行平行程式 file.x #!/bin/csh #@ executable = /usr/bin/poe #@ network.mpi= csss,shared,us #@ arguments=/your-working-directory/file.x #@ output = outp8 #@ error = outp8 #@ job_type = parallel #@ class = 8cpu #@ tasks_per_node = 8 #@ node = 1 #@ queue 10

11 由於 IBM p690 每個 Node 含有 32 棵 1.3 GHz PU 共用 128GB 的記憶體 class = p690 PU 的分組別, 使用 llclass 指令可以看到分組別 : 32cpu (PU 時間上限為 36 小時, 2 個 Node 共有 64 顆 PU) 16cpu (PU 時間上限為 72 小時, 1 個 Node 共有 32 顆 PU) 8cpu (PU 時間上限為 72 小時, 6 個 Node 共有 188 顆 PU) 4cpu (PU 時間上限為 96 小時, 6 個 Node 共有 188 顆 PU) serial (PU 時間上限為 168 小時, 1 個 Node 共有 32 顆 PU) tasks_per_node=8 是說明一個 Node 選用 8 棵 PU node=1 是說明要用一個 Node, 一共 8 棵 PU IBM 電腦系統的平行程式的執行指令 要在 IBM 電腦系統上執行平行程式, 使用者在備妥 LoadLeveler 的 job command file 之後, 就可以使用 llsubmit 指令將該 job command file 交給該系統排隊等候執行 例如上一節的 job command file 例子 jobp4 即可用下述指令交付執行 : llsubmit jobp4 工作交付之後, 該工作執行的情形可用 llq 指令查詢 要縮小查詢的範圍可在 llq 指令之後加 上 grep 指令敘明要查詢的 class 或 user id 例如上一個例子 jobp4 所選用的分組別為 medium, 就可用下述指令進行查詢 : llq grep medium llq 顯示之內容有下列事項 : job_id user_id submitted status priority class running on ivy u43ycc00 8/13 11:24 R 50 medium ivy39 ivy u50pao00 8/13 20:12 R 50 short ivy35 其中 job_id user_id submitted status 是 LoadLeveler 給交付的工作編定的工作代號是使用者的 login name 是交付工作的時刻, 月 / 日時 : 分是工作執行的情形 R 表 Running I 表 Idle (=waiting in queue) 11

12 Priority lass Running on ST 表 Start execution NQ 表 Not Queued, 還在隊伍之外是交付工作的優先次序, 不用更動它是 PU 分組別是執行交付工作的第一個 PU 代號 工作交付執行之後, 如果要中止該工作的執行可用 llcancel 指令殺掉該工作 llcancel job_id 此處的 job_id 就是使用 llq 指令所顯示之使用者交付工作的工作代號 執行過 llcancel 指令 之後, 再使用 llq 指令就可以看出該工作已經消失不見了 12

13 1.4 在 P luster 上如何使用 MPI 首先, 使用 MPIH 的 shell 用戶要在自己 home directory 的.cshrc 檔裏加入下列路徑, 這樣才能夠抓得到 include file (mpif.h mpi.h) 編譯指令 (mpif77 mpicc) MPI library 和 DQS 指令 不同的 P luster 這些存放的路徑可能不同, 要向該系統的管理人詢問 其路徑設定如下 : setenv PGI /usr/local/pgi set path = (. ~ /usr/local/pgi/linux86/bin $path) set path = ( /package/dqs _hpcserv2/bin $path) set path = ( /package/mpich _hpcserv2/bin $path) 其中第一行是 PGI 公司 (Portland Group Inc.) 軟體存放的路徑, 第二行是 PGI 公司 Fortran77 編譯器 pgf77 存放的路徑, 第三行是 DQS 批次工作排程軟体存放的路徑, 第四行是 MPIH 編 譯系統存放的路徑 沒有購用 PGI 公司的軟體時前面兩行可以省略 P luster 上的 MPI Fortran 程式編譯指令 MPIH 的 Fortran77 平行程式編譯器叫做 mpif77, 其底層是使用 GNU 的 g77 來編譯, 因此可以使用 g77 的調適選項 舉例如下 : mpif77 -O3 -o file.x file.f 其中選項 -O3 是選用 g77 最高層次的調適選項 -o file.x 是指定編譯產生的執行檔為 file.x 沒有指定時, 內定的執行檔為 a.out file.f 是 Fortran77 平行程式 如果選用 PGI 公司的 MPI 平行程式編譯器 mpif77, 其底層是使用該公司的 pgf77 來編譯, 因此可以使用 pgf77 及 pgf90 的調適選項 其 makefile 舉例如下 : 13

14 OBJ = file.o EXE = file.x MPI = /home/package/mpich_pgi LIB = $(MPI)/lib/LINUX/ch_p4 LFLAG = -L$(LIB) -lfmpich MPIF77 = $(MPI)/bin/mpif77 OPT = -O2 -I$(MPI)/include $(EXE) : $(OBJ) $(MPIF77) $(LFLAG) -o $(EXE) $(OBJ) $(LIB).f.o : $(MPIF77) $(OPT) -c $< 備妥 makefile 之後, 只要下 make 指令就會開始程式的編譯工作 P luster 上的 Job command file 如果該 P cluster 是採用 DQS 排程軟體來安排批次工作時, 要在其上執行平行程式, 使用者必須備妥 DQS 的 job command file 例如, 下面這個 job command file 叫做 jobp4, 它要在四個 PU 上執行平行程式 hubksp : #!/bin/csh #$ -l qty.eq.4 #$ -N HUP4 #$ -A user_id #$ -cwd #$ -j y cat $HOSTS_FILE > MPI_HOST mpirun -np 4 -machinefile MPI_HOST hubksp >& outp4 其中 #!/bin/csh 是說明這是個 shell script #$ -l qty.eq.4 是向 DQS 要求四個 PU,qty 是數量 (quantity) #$ -N HUP4 是說明這個工作的名字 (Name) 叫做 HUP4 #$ -A user_id 是說明付費帳號 (Account) 就是使用者帳號 #$ -cwd 是說明要在現在這個路徑 (working directory) 上執行程式 內定的路徑是 home directory #$ -j y 是說明錯誤訊息要輸出到標準輸出檔 $HOST_FILE 是 DQS 安排給這項工作的 node list -np 4 hubksp 是告訴 mpirun 要在四個 PU 上執行平行程式 hubksp >& outp4 是要把標準輸出檔寫入 outp4 14

15 1.4.3 P luster 上的平行程式執行指令 要在 P cluster 上執行平行程式, 使用者在備妥 DQS 的 job command file 之後, 就可以使用 qsub32 指令將該 job command file 交給 P cluster 排隊等候執行 例如上一節的 job command file 例子 jobp4 即可用下述指令交付執行 : qsub32 jobp4 工作交付之後, 可以使用 qstat32 指令 ( 不加參數 ) 查詢整個 cluster 交付工作執行的情形, 使用 qstat32 -f 指令查詢整個 cluster 各個 node 的狀況 上述指令 qsub32 jobp4 之後使用 qstat32 指令顯示的內容如下 : c00tch00 HUP4 hpcs :1 r RUNNING 02/26/99 10:51:23 c00tch00 HUP4 hpcs :1 r RUNNING 02/26/99 10:51:23 c00tch00 HUP4 hpcs :1 r RUNNING 02/26/99 10:51:23 c00tch00 HUP4 hpcs :1 r RUNNING 02/26/99 10:51: Pending Jobs c00tch00 RAD5 70 0:2 QUEUED 02/26/99 19:24:32 第一欄是 user_id, 第二欄是交付工作的名稱, 第三欄是 PU 代號, 第四欄是 DQS 替交付的工作編定的工作編號 job_id(62), 第五欄 0:1 的 0 是交付工作的優先序號, 0:1 的 1 是該用戶交付的第一個工作, 第六欄的 r 和第七欄的 RUNNING 表示該工作正在執行中, 最後是該工作交付時的時刻, 月 / 日 / 年時 : 分 : 秒 排隊等待執行的工作則出現在 Pending Jobs 之列, 對應 RUNNING 的欄位則為 QUEUED 工作交付執行之後, 如果要中止該工作的執行可用 qdel32 指令殺掉該工作 Qdel32 job_id 此處的 job_id 就是使用 qstat32 指令所顯示之第四欄 執行過 qdel32 指令之後, 再使用 qstat32 指令就可以看出該工作已經消失不見了 15

16 第二章無邊界資料交換的平行程式 最簡單的平行程式就是無邊界資料交換的平行程式 本章將利用一個很簡單的循序程式 (sequential program) 使用 MPI 指令加以平行化, 並比較其計算結果以資驗證 2.1 節介紹六個 MPI 基本指令 MPI_INIT MPI_FINALIZE MPI_OMM_SIZE MPI_OMM_RANK MPI_SEND MPI_REV 2.2 節介紹無邊界資料交換的循序程式 T2SEQ 2.3 節說明使用這六個 MPI 基本指令平行化循序程式 T2SEQ 而成為平行程式 T2P 2.4 節介紹另外四個常用的 MPI 指令 MPI_SATTER MPI_GATHER MPI_REDUE MPI_ALLREDUE 2.5 節是使用這些指令平行化循序程式 T2SEQ 而成為平行程式 T2DP 16

17 2.1 MPI 基本指令 MPI 的基本指令有下列六個, 將於本節分段加以介紹 MPI_INIT, MPI_FINALIZE, MPI_OMM_SIZE, MPI_OMM_RANK, MPI_SEND, MPI_REV mpif.h include file 使用 MPI 撰寫 Fortran 平行程式時, 必須在每一個程式 ( 包括主程式和副程式 ) 的宣告段落裏加上 INLUDE 'mpif.h' 陳述 (statement) mpif.h 檔案裏含有編譯 MPI 平行程式所必須的 MPI 字彙與 MPI 常數 (constant) 例如 : PROGRAM DRIVER IMPLIIT REAL*8... INLUDE 'mpif.h'... ALL HEF(... )... STOP END SUBROUTINE HEF(... ) IMPLIIT REAL*8... INLUDE 'mpif.h' RETURN END 讀者可以在 MPI 軟體所在之路徑裏查看 mpif.h 的內容 不同廠商設定的 MPI 常數也許不盡相同, 但是所使用的 MPI 字彙則是完全一致 MPI_INIT, MPI_FINALIZE 在叫用 (ALL) 其他 MPI 函數或副程式之前必須先叫用 MPI_INIT 副程式, 來啟動該程式在多個 PU 上的平行計算工作 在程式結束 (STOP) 之前必須叫用 MPI_FINALIZE 副程式, 以結束平行計算工作 所以 MPI_INIT 和 MPI_FINALIZE 在主程式裏只要叫用一次就夠了, 例如 : 17

18 PROGRAM T2P PARAMETER (.. ) INLUDE 'mpif.h' REAL*8... INTEGER... ALL MPI_INIT(IERR)... ALL MPI_FINALIZE(IERR) STOP END 所有 MPI 副程式引數 (argument) 的資料類別除了資料名稱之外, 其餘的都是整數 (integer), 傳回的引數 IERR 其值為零時是正常結束, 否則就有錯誤發生 MPI_OMM_SIZE, MPI_OMM_RANK 通常在叫用過 MPI_INIT 之後, 就必須叫用 MPI_OMM_SIZE 以得知參與平行計算的 PU 個數 (NPRO), 及叫用 MPI_OMM_RANK 以得知我是第幾個 PU (MYID), 第幾個 PU 是從 0 開始起算 所以第一個 PU 的 MYID 值為零, 第二個 PU 的 MYID 值為 1, 第三個 PU 的 MYID 值為 2, 餘類推 通常要在幾個 PU 上作平行計算是在下執行命令時決定的, 而不是在程式裏事先設定 當然, 使用者也可以在程式裏事先設定要在幾個 PU 上作平行計算, 其意義只供程式人員做參考, 實際上使用幾個 PU 作平行計算是根據 job command file 裏 min_processors 和 max_processors 的設定值, 或 -np 的設定值 MPI_OMM_SIZE 和 MPI_OMM_RANK 的叫用格式如下 : ALL MPI_OMM_SIZE (MPI_OMM_WORLD, NPRO, IERR) ALL MPI_OMM_RANK (MPI_OMM_WORLD, MYID, IERR) 引數 MPI_OMM_WORLD 是 MPI 內定的 (default) communicator, 參與該程式平行計算的全部 PU 都是屬於同一個 communicator 屬於同一個 communicator 的各個 PU 之間才可以傳送資料 MPI 1.2 版不具備 PU 的取得與控制功能, 參與平行計算的 PU 顆數從程式開始執行到程式結束都是固定不變的 因此, 這兩個 MPI 副程式在一個程式裏只要叫用一次就可以了 例如 : PROGRAM T2P PARAMETER (.. ) INLUDE 'mpif.h' 18

19 REAL*8... INTEGER NPRO, MYID ALL MPI_INIT (IERR) ALL MPI_OMM_SIZE (MPI_OMM_WORLD, NPRO, IERR) ALL MPI_OMM_RANK (MPI_OMM_WORLD, MYID, IERR)... ALL MPI_FINALIZE (IERR) STOP END MPI_SEND, MPI_REV 參與平行計算的各個 PU 之間的資料傳送方式有兩種, 一種叫做 ' 點對點通訊 ' (point to point communication), 另外一種叫做 ' 集體通訊 ' (collective communication) 此處先介紹 ' 點對點通訊 ' 類的 MPI_SEND 和 MPI_REV, 其他常用的 ' 點對點通訊 ' 及 ' 集體通訊 ' 指令容後再介紹 一個 PU 與另外一個 PU 之間的資料傳送屬於 ' 點對點通訊 ', 送出資料的 PU 要叫用 MPI_SEND 來送資料, 而收受資料的 PU 要叫用 MPI_REV 來收資料 一個 MPI_SEND 必須要有一個對應的 MPI_REV 與之配合, 才能完成一份資料的傳送工作 MPI_SEND 的叫用格式如下 : & ALL MPI_SEND (DATA, IOUNT, DATA_TYPE, IDEST, ITAG, MPI_OMM_WORLD, IERR) 引數 DATA 要送出去的資料起點, 可以是純量 (scalar) 或陣列 (array) 資料 IOUNT 要送出去的資料數量, 當 IOUNT 的值大於一時,DATA 必須是陣列 DATA_TYPE 是要送出去的資料類別,MPI 內定的資料類別如表 1.1 IDEST 是收受資料的 PU id ITAG 要送出去的資料標籤 MPI data types MPI_HARATER MPI_LOGIAL MPI_INTEGER MPI_REAL, MPI_REAL4 MPI_REAL8, MPI_DOUBLE_PREISION MPI_OMPLEX, MPI_OMPLEX8 MPI_OMPLEX16, MPI_DOUBLE_OMPLEX Fortran data types HARATER LOGIAL INTEGER REAL, REAL*4 REAL*8, DOUBLE PREISION OMPLEX, OMPLEX*8 OMPLEX*16 表 1.1 MPI data types 19

20 MPI_REV 的叫用格式如下 : & ALL MPI_REV (DATA, IOUNT, DATA_TYPE, ISR, ITAG, MPI_OMM_WORLD, ISTATUS, IERR) 引數 DATA IOUNT DATA_TYPE ISR ITAG ISTATUS 是要收受的資料起點是要收受的資料數量是要收受的資料類別是送出資料的 PU id 是要收受的資料標籤是執行 MPI_REV 副程式之後的狀況 ISTATUS 為一整數陣列, 該陣列的長度為在 mpif.h 裏已經設定的常數 MPI_STATUS_SIZE, 寫法如下 INTEGER ISTATUS(MPI_STATUS_SIZE) 一個 PU 同時要收受多個 PU 送來的資料時, 若不依照特定的順序, 而是先到先收, 則其 指令為 ALL MPI_REV( BUFF, IOUNT, DATA_TYPE, MPI_ANY_SOURE, ITAG, 1 MPI_OMM_WORLD, ISTATUS, IERR) 若要判別送出該資料的 PU id 時就要用到 ISTATUS 變數如下 ISR = ISTATUS( MPI_SOURE ) MPI 在傳送資料 (MPI_SEND MPI_REV) 時, 是以下列四項構成其 ' 信封 ' (envelope), 用 以識別一件訊息 (message) 1. 送出資料的 PU id 2. 收受資料的 PU id 3. 資料標籤 4. communicator 所以一個 PU 送給另外一個 PU 多種資料時, 不同的資料要用不同的資料標籤, 以資識別 20

21 2.2 無邊界資料交換的循序程式 T2SEQ T2SEQ 是個無邊界資料交換的循序程式, 在 test data generation 段落裏設定陣列 B D 的值, 然後把這些陣列寫到磁檔上 其目的是便利往後的範例程式可以讀入同一組資料作平行計算, 用來驗證其計算的結果是否正確 這個程式的計算部份只有一個 DO loop, 而且該 loop 裏只有兩個計算陳述, 其目的是方便往後說明如何將這一類 DO loop 平行化 實際的計算程式也許有數百個或數千個 DO loop, 但是其平行化的方法是一樣的 PROGRAM T2SEQ PARAMETER (NTOTAL=200) REAL*8 A(NTOTAL), B(NTOTAL), (NTOTAL), D(NTOTAL), SUMA test data generation and write to file 'input.dat' DO I=1,NTOTAL B(I)=3.D0/DFLOAT(I)+1.0 (I)=2.D0/DFLOAT(I)+1.0 D(I)=1.D0/DFLOAT(I)+1.0 OPEN(7,FILE='input.dat',FORM='UNFORMATTED') WRITE(7) B WRITE(7) WRITE(7) D LOSE(7) read 'input.dat', compute and write out the result OPEN(7,FILE='input.dat',STATUS='OLD',FORM='UNFORMATTED') READ (7) B READ (7) READ (7) D SUMA=0.0 DO I=1,NTOTAL A(I)=B(I)+(I)*D(I) SUMA=SUMA+A(I) 21

22 WRITE(*,101) (A(I),I=1,NTOTAL,5) 101 FORMAT(10F8.3) WRITE(*,102) SUMA 102 FORMAT('SUM of array A =', E15.5) STOP END 循序程式 T2SEQ 的測試結果如下 : SUM of array A =.43855E+03 22

23 2.3 資料不切割的平行程式 T2P 平行程式的切割 (decomposition/partition) 方式有兩種 一種是計算切割而資料不切割, 另外一種是計算和資料都切割 前一種切割方式不能夠節省記憶體的使用量是其缺點, 但是陣列的描述與循序版 (sequential version) 完全相同, 程式容易閱讀也容易維護是其優點 後一種切割方式能夠節省記憶體的使用量是其最大優點, 但是陣列的描述與循序版差異較大, 程式的閱讀與維護比較困難是其缺點 如何將循序程式 T2SEQ 平行化呢? 這一節先介紹 ' 計算切割而資料不切割 ' 的方法,2.5 節再介紹 ' 計算及資料同時切割 ' 的方法 假如 T2SEQ 程式要在四個 PU 上平行計算而資料不切割時, 就把一維陣列 A B D 均分為四段, 各個 PU 負責計算其中的一段, 分工合作完成整個計算工作 此處是利用一個 STARTEND 副程式來計算各個 PU 負責計算段落的起迄 index 它是把第一段分給 PU0, 第二段分給 PU1, 第三段分給 PU2, 餘類推 如圖 2.1 所示 : computing partition without data partition cpu0 cpu1 cpu2 cpu3 istart iend ntotal istart iend ntotal istart iend ntotal istart iend array element inside territory array element outside territory 圖 2.1 計算切割而資料不切割的示意圖 圖 2.1 裏符號 代表轄區內的陣列元素, 符號 代表轄區外的陣列元素, 各個 PU 負責計算的範圍是從該 PU 的 istart 到 iend 由於 MPI 1.2 版不具備平行輸出入 (Parallel I/O) 的功能, 所以輸入資料由 PU0(MYID 值為零 ) 讀入後, 利用一個 DO loop 分段傳送 (MPI_SEND) 給其他三個 PU, 而其他三個 PU 23

24 (MYID 值大於零 ) 則接收由 PU0 送來給該 PU 的陣列片段 請留意, 傳送不同的陣列片段要使用不同的資料標籤 (ITAG), 每一個 MPI_SEND 一定有一個對應的 MPI_REV 每一個 PU 算完自己負責的段落後, 把計算的結果 ( 陣列 A 的一部份 ) 傳送給 PU0,PU0 利用一個 DO loop 把其他三個 PU 送來的陣列片段逐一接收下來 然後由 PU0 單獨計算整個 A 陣列各個元素的和 SUMA, 再把 A 陣列和 SUMA 列印出來 PROGRAM T2P computing partition without data partition INLUDE 'mpif.h' PARAMETER (NTOTAL=200) REAL*8 A(NTOTAL), B(NTOTAL), (NTOTAL), D(NTOTAL), SUMA INTEGER NPRO, MYID, ISTATUS(MPI_STATUS_SIZE), ISTART, IEND, 1 OMM, GOUNT(0:31), GSTART(0:31), GEND(0:31) ALL MPI_INIT (IERR) ALL MPI_OMM_SIZE (MPI_OMM_WORLD, NPRO, IERR) ALL MPI_OMM_RANK (MPI_OMM_WORLD, MYID, IERR) ALL STARTEND ( NPRO, 1, NTOTAL, GSTART, GEND, GOUNT) ISTART=GSTART(MYID) IEND=GEND(MYID) PRINT *,' NPRO,MYID,ISTART,IEND=', NPRO,MYID,ISTART,IEND OMM=MPI_OMM_WORLD Read input data and distribute data to each processor IF (MYID.EQ.0) THEN OPEN(7,FILE='input.dat',STATUS='OLD',FORM='UNFORMATTED') READ (7) B READ (7) READ (7) D DO IDEST=1,NPRO-1 ISTART1=GSTART(IDEST) KOUNT1=GOUNT(IDEST) ALL MPI_SEND (B(ISTART1), KOUNT1, MPI_REAL8, IDEST, & 10, OMM, IERR) ALL MPI_SEND ((ISTART1), KOUNT1, MPI_REAL8, IDEST, 24

25 & 20, OMM, IERR) ALL MPI_SEND (D(ISTART1), KOUNT1, MPI_REAL8, IDEST, & 30, OMM, IERR) ELSE KOUNT=GOUNT(MYID) ISR=0 ALL MPI_REV (B(ISTART), KOUNT, MPI_REAL8, ISR, 10, 1 OMM, ISTATUS, IERR) ALL MPI_REV ((ISTART), KOUNT, MPI_REAL8, ISR, 20, 1 OMM, ISTATUS, IERR) ALL MPI_REV (D(ISTART), KOUNT, MPI_REAL8, ISR, 30, 1 OMM, ISTATUS, IERR) ENDIF ompute, collect computed result and write out the result DO I=1,NTOTAL DO I=ISTART, IEND A(I)=B(I)+(I)*D(I) ITAG=110 IF (MYID.NE.0) THEN KOUNT=GOUNT(MYID) IDEST=0 ALL MPI_SEND (A(ISTART), KOUNT, MPI_REAL8, IDEST, ITAG, & OMM, IERR) ELSE DO ISR=1,NPRO-1 ISTART1=GSTART(ISR) KOUNT1=GOUNT(ISR) ALL MPI_REV (A(ISTART1), KOUNT1, MPI_REAL8, ISR, ITAG, & OMM, ISTATUS, IERR) ENDIF IF(MYID.EQ.0) THEN WRITE(*,101) (A(I),I=1,NTOTAL,5) SUMA=0.0 25

26 DO I=1,NTOTAL SUMA=SUMA+A(I) WRITE(*,102) SUMA ENDIF 101 FORMAT(10F8.3) 102 FORMAT('SUM of array A =', E15.5) ALL MPI_FINALIZE (IERR) STOP END SUBROUTINE STARTEND(NPRO,IS1,IS2,GSTART,GEND,GOUNT) INTEGER ID,NPRO,IS1,IS2,GSTART(0:31),GEND(0:31),GOUNT(0:31) LENG=IS2-IS1+1 IBLOK=LENG/NPRO IR=LENG-IBLOK*NPRO DO I=0,NPRO-1 IF(I.LT.IR) THEN GSTART(I)=IS1+I*(IBLOK+1) GEND(I)=GSTART(I)+IBLOK ELSE GSTART(I)=IS1+I*IBLOK+IR GEND(I)=GSTART(I)+IBLOK-1 ENDIF IF(LENG.LT.1) THEN GSTART(I)=1 GEND(I)=0 ENDIF GOUNT(I)=GEND(I)-GSTART(I)+1 END 資料不切割平行程式 T2P 的測試結果如下 : ATTENTION: nodes allocated by LoadLeveler, continuing... NPRO,MYID,ISTART,IEND= NPRO,MYID,ISTART,IEND= NPRO,MYID,ISTART,IEND= NPRO,MYID,ISTART,IEND=

27 SUM of array A =.43855E+03 這一種平行程式的寫法叫做 SPMD (Single Program Multiple Data) 程式 這個程式在多顆 PU 上平行時, 每一顆 PU 都是執行這一個程式, 有些地方是用 rank (MYID) 來判斷要執行條件陳述 (IF statement) 裏的那一個區段, 有些地方是用不同的 index 起迄位置來執行 DO loop 的某一段落 沒有使用 index 或 rank 來區分執行的段落部份, 則每一顆 PU 都要執行 所以 T2P 程式裏, 每一顆 PU 都會執行該程式一開始的下列陳述 : ALL MPI_INIT (IERR) ALL MPI_OMM_SIZE (MPI_OMM_WORLD, NPRO, IERR) ALL MPI_OMM_RANK (MPI_OMM_WORLD, MYID, IERR) ALL STARTEND (NPRO, 1, NTOTAL, GSTART,GEND,GOUNT) ISTART=GSTART(MYID) IEND=GEND(MYID) PRINT *,' NPRO,MYID,ISTART,IEND=', NPRO,MYID,ISTART,IEND 雖然執行的陳述完全一樣, 但是每一個 PU 得到的 MYID 值卻不一樣, 所得到的 ISTART 和 IEND 也就隨之而異, 從所列印出來的結果可以得到印證 這樣一個 SPMD 平行程式每一個 PU 都有工作做 PU0 會多做資料的讀和寫的工作 但是一般計算程式資料輸出入所佔的時間比例很小, 所以各個 PU 的負載相當均勻, 也沒有 master PU 和 slave PU 的區別 27

28 2.4 MPI_SATTER,MPI_GATHER,MPI_REDUE MPI_SATTER MPI_GATHER MPI_ALLGATHER MPI_REDUE MPI_ALLREDUE 都是屬於 ' 集體通訊 ' 類副程式 這一類的資料傳輸, 凡是屬於同一個 communicator 的每一個 PU 都要參與運作 所以使用這一種指令時, 每一個 PU 都必須叫用同一個副程式 MPI_SATTER,MPI_GATHER MPI_SATTER 的叫用格式如下 : IROOT = 0 ALL MPI_SATTER (T, N, MPI_REAL8, B, N, MPI_REAL8, & IROOT, MPI_OMM_WORLD, IERR) MPI_SATTER 是 IROOT PU 把一個陣列 T 等分為 NPRO 段,(NPRO= 參與平行計算的 PU 數量 ), 每一段資料長度為 N, 依 PU id 的順序分送給每一個 PU ( 包括 IROOT PU 在內 ) 它是把第一段分給 PU0, 第二段分給 PU1, 第三段分給 PU2, 餘類推 如圖 2.2 所示 : PU0 T T1 T2 T3 T4 PU0 B T1 PU1 SATTER PU1 B T2 > PU2 PU2 B T3 PU3 PU3 B T4 圖 2.2 MPI_SATTER 示意圖 請注意每一段資料必須等長 其引數依序為 T 是待送出陣列的起點 N 是送給每一個 PU 的資料數量 MPI_REAL8 是待送出資料的類別 B N MPI_REAL8 IROOT 是接收資料存放的起點, 如果 N 值大於一時,B 必須是個陣列是接收資料的數量是接收資料的類別是送出資料的 PU id 28

29 MPI_GATHER 的叫用格式如下 : IDEST = 0 ALL MPI_GATHER (A, N, MPI_REAL8, T, N, MPI_REAL8, & IDEST, MPI_OMM_WORLD, IERR) MPI_GATHER 與 MPI_SATTER 的動作剛好相反, 是 IDEST PU 收集每一個 PU 送給它的陣列 A, 依 PU id 的順序存入陣列 T 裏頭 也就是從 PU0 收到的 N 個陣列元素存入 T 陣列的第一段, 從 PU1 收到的 N 個陣列元素存入 T 陣列的第二段, 從 PU2 收到的 N 個陣列元素存入 T 陣列的第三段, 餘類推 如圖 2.3 所示 : PU0 T T1 T2 T3 T4 PU0 A T1 PU1 GATHER PU1 A T2 < PU2 PU2 A T3 PU3 PU3 A T4 圖 2.3 MPI_GATHER 示意圖 請注意每一段資料必須等長 MPI_GATHER 的引數依序為 A 是待送出的資料起點, 如果 N 值大於一時,A 必須是個陣列 N 是待送出資料的數量 MPI_REAL8 是待送出資料的類別 T 是接收資料存放的陣列起點 N 是接收來自各個 PU 的資料數量 MPI_REAL8 是接收資料的類別 IDEST 是收集資料的 PU id MPI_ALLGATHER 的叫用格式如下 : ALL MPI_ALLGATHER (A, N, MPI_REAL8, T, N, MPI_REAL8, & MPI_OMM_WORLD, IERR) 29

30 MPI_ALLGATHER 與 MPI_GATHER 的運作功能相似,MPI_GATHER 是把運作的結果存入 指定的一個 PU, 而 MPI_ALLGATHER 則是把運作的結果存入每一個 PU PU0 T T1 T2 T3 T4 PU0 A T1 PU1 T T1 T2 T3 T4 ALLGATHER PU1 A T2 PU2 T T1 T2 T3 T4 PU2 A T3 PU3 T T1 T2 T3 T4 PU3 A T4 圖 2.4 MPI_ALLGATHER 示意圖 MPI_REDUE, MPI_ALLREDUE 另外一種集體資料傳輸功能叫作 ' 縮減運作 ' (reduction operation), 例如把各個 PU 算出來的部份和 (partial sum) 加總, 或找出各個 PU 上某一個變數的最大值或最小值 MPI_REDUE 運作的結果只存放在指定的 PU (IROOT) 裏,MPI_ALLREDUE 則是把運作的結果存放在每一個 PU 裏 MPI_REDUE 和 MPI_ALLREDUE 叫用格式如下 : IROOT = 0 ALL MPI_REDUE ( SUMA, SUMALL, KOUNT, MPI_REAL8, MPI_SUM, IROOT, & MPI_OMM_WORLD, IERR) & ALL MPI_ALLREDUE( SUMA, SUMALL, KOUNT, MPI_REAL8, MPI_SUM, MPI_OMM_WORLD, IERR) 引數 SUMA 是待運作 ( 累加 ) 的變數 SUMALL 是存放運作 ( 累加 ) 後的結果 ( 把各個 PU 上的 SUMA 加總 ) KOUNT 是待運作 ( 累加 ) 的資料個數 MPI_REAL8 是 SUMA 和 SUMALL 的資料類別 MPI_SUM 是運作涵數, 可以選用的涵數如表 2.1 IROOT 是存放運作結果的 PU_id 30

31 MPI 指令 Operation Data type MPI_SUM sum 累加 MPI_INTEGER, MPI_REAL, MPI_REAL8, MPI_PROD product 乘積 MPI_OMPLEX, MPI_OMPLEX16 MPI_MAX MPI_MIN MPI_MAXLO MPI_MINLO MPI_LAND MPI_LOR MPI_LXOR MPI_BAND MPI_BOR MPI_BXOR maximum 最大值 minimum 最小值 max value and location min value and location logical AND logical OR logical exclusive OR binary AND binary OR binary exclusive OR MPI_INTEGER, MPI_REAL, MPI_REAL8, MPI_DOUBLE_PREISION MPI_2INTEGER, MPI_2REAL, MPI_2 DOUBLE_PREISION MPI_LOGIAL MPI_INTEGER, MPI_BYTE 表 2.1 MPI Reduction Function MPI_REDUE 的運作方式如圖 2.5 所示,MPI_ALLREDUE 的運作方式如圖 2.6 所示 PU0 SUMA PU0 SUMALL PU1 SUMA REDUE PU1 > PU2 SUMA (MPI_SUM) PU2 PU3 SUMA PU3 圖 2.5 MPI_REDUE 示意圖 PU0 SUMA PU0 SUMALL PU1 SUMA ALLREDUE PU1 SUMALL > PU2 SUMA (MPI_SUM) PU2 SUMALL PU3 SUMA PU3 SUMALL 圖 2.6 MPI_ALLREDUE 示意圖 此處 SUMALL 等於各個 PU 上 SUMA 對應項目之和 31

32 2.5 資料切割的平行程式 T2DP 程式 T2DP 是一個計算及資料都切割的平行程式 這個程式在 NP 個 PU 上平行計算時, 陣列 A B D 的長度只需要原來長度 NTOTAL 的 NP 分之一 但是輸入資料 B D 陣列的長度是 NTOTAL, 所以要再設定一個長度為 NTOTAL 的暫用陣列 T 來輪流存放讀入資料陣列 B D 及存放輸出資料整個陣列 A PU0 每讀入一個輸入陣列就利用 MPI_SATTER 分段分送給每一個 PU IROOT=0 ALL MPI_SATTER (T, N, MPI_REAL8, B, N, MPI_REAL8, & IROOT, MPI_OMM_WORLD, IERR) 各個 PU 把分配給它的陣列段落計算完畢之後, 就利用 MPI_GATHER 把計算的結果 A 陣列 片段送回給 PU0 IDEST=0 ALL MPI_GATHER (A, N, MPI_REAL8, T, N, MPI_REAL8, & IDEST, MPI_OMM_WORLD, IERR) 比較 T2P 和 T2DP 兩個程式就可以看出利用 MPI_SATTER,MPI_GATHER 的程式要簡 潔得多 但是使用這種指令時分送給每一個 PU 的陣列長度必須相等, 而使用 MPI_SEND MPI_REV 指令時就沒有這種限制 由於 Fortran 77 的 DIMENSION 必須是常數, 所以必須在程式裏先設定 PU 的個數 NP, 切割後的陣列長度 N 等於 NTOTAL / NP, 使用 PARAMETER 陳述設定如下 : PARAMETER (NTOTAL=200, NP=4, N=NTOTAL/NP) 此處 NTOTAL 必須能被 NP 整除 於是, 陣列的 dimension 由 NTOTAL 改為 N : REAL*8 A(N), B(N), (N), D(N), T(NTOTAL) 現在, 每一個 PU 執行的 DO loop 範圍是從 1 到 N, 所以算出來的 SUMA 是 A 陣列 NTOTAL 個元素裏 NP 分之一個陣列元素的和, 是為部份和 (partial sum) 所以此處叫用 MPI_REDUE 把各個 PU 的部份和 SUMA 加總, 存放在 PU0 的 SUMALL 裏 32

33 IROOT=0 ALL MPI_REDUE (SUMA, SUMALL, 1, MPI_REAL8, MPI_SUM, & IROOT, MPI_OMM_WORLD, IERR) 整個計算及資料都切割的平行程式 T2DP 如下 : PROGRAM T2DP Data & omputational Partition Using MPI_SATTER, MPI_GATHER NP=4 must be modified when run on other than 4 processors index n n n n o o o o o o o o o o o o o o o o o o o o o o o o P0 P1 P2 P3 PARAMETER (NTOTAL=200, NP=4, N=NTOTAL/NP) INLUDE 'mpif.h' REAL*8 A(N), B(N), (N), D(N), T(NTOTAL), SUMA, SUMALL INTEGER NPRO, MYID, ISTATUS(MPI_STATUS_SIZE), ISTART, IEND IF (MOD(NTOTAL,NP).NE.0) THEN PRINT *,' NTOTAL IS NOT DIVISIBLE BY NP, PROGRAM WILL STOP' PRINT *,' NTOTAL,NP=', NTOTAL,NP STOP ENDIF ALL MPI_INIT (IERR) ALL MPI_OMM_SIZE ( MPI_OMM_WORLD, NPRO, IERR) ALL MPI_OMM_RANK ( MPI_OMM_WORLD, MYID, IERR) PRINT *,' NPRO,MYID=', NPRO, MYID Read and distribute input data IF (MYID.EQ.0) THEN OPEN(7,FILE='input.dat',STATUS='OLD',FORM='UNFORMATTED') READ (7) T! read array B ENDIF 33

34 IROOT=0 ALL MPI_SATTER (T, N, MPI_REAL8, B, N, MPI_REAL8, & IROOT, MPI_OMM_WORLD, IERR) IF(MYID.EQ.0) THEN READ (7) T! read array ENDIF ALL MPI_SATTER (T, N, MPI_REAL8,, N, MPI_REAL8, & IROOT, MPI_OMM_WORLD, IERR) IF(MYID.EQ.0) THEN READ (7) T! read array D ENDIF ALL MPI_SATTER (T, N, MPI_REAL8, D, N, MPI_REAL8, & IROOT, MPI_OMM_WORLD, IERR) ompute, GATHER computed data, and write out the result SUMA=0.0 DO I=1,NTOTAL DO I=1,N A(I)=B(I)+(I)*D(I) SUMA=SUMA+A(I) IDEST=0 ALL MPI_GATHER (A, N, MPI_REAL8, T, N, MPI_REAL8, & IDEST, MPI_OMM_WORLD, IERR) ALL MPI_REDUE (SUMA, SUMALL, 1, MPI_REAL8, MPI_SUM, & IDEST, MPI_OMM_WORLD, IERR) IF (MYID.EQ.0) THEN WRITE(*,101) (T(I),I=1,NTOTAL,5) WRITE(*,102) SUMALL ENDIF 101 FORMAT(10F8.3) 102 FORMAT('SUM of array A =', E15.5) ALL MPI_FINALIZE (MPI_ERR) STOP END 計算及資料都切割的平行程式 T2DP 測試結果如下 : 34

35 ATTENTION: nodes allocated by LoadLeveler, continuing... NPRO,MYID= 4 0 NPRO,MYID= 4 1 NPRO,MYID= 4 2 NPRO,MYID= SUM of array A =.43855E+03 35

36 第三章需要邊界資料交換的平行程式 最常見的平行程式就是在計算的過程中需要邊界資料交換的平行程式 這一章將利用一個需要邊界資料交換的循序程式使用 MPI 指令加以平行化, 並比較其計算結果以資驗證 3.1 節將介紹兩個 MPI 指令 MPI_SENDREV MPI_BAST MPI_SENDREV 是用來交換邊界資料, 而 MPI_BAST 則是用來分送輸入資料 3.2 節介紹需要邊界資料交換的循序程式 T3SEQ 3.3 節說明使用 MPI_SENDREV MPI_SEND MPI_REV 指令平行化循序程式 T3SEQ 而成為平行程式 T3P 然後說明使用 MPI_BAST 取代 T3P 所使用的 MPI_SEND 和 MPI_REV 的方法 T3P 是計算切割而資料不切割的平行程式 3.4 節說明交換一個陣列元素的資料切割平行程式 T3DP_1 3.5 節說明交換兩個陣列元素的資料切割平行程式 T3DP_2 36

37 3.1 MPI_SENDREV, MPI_BAST MPI_SENDREV 是屬於 ' 點對點通訊 ' 類副程式, 而 MPI_BAST 是屬於 ' 集體通訊 ' 類副程式 MPI_SENDREV 甲 PU 送出一些資料給乙 PU, 又要接受丙 PU 送來另外一些資料時, 可以叫用一個 MPI_SENDREV 副程式來做這兩件工作 其作用等於一個 MPI_SEND 加一個 MPI_REV 下面這個指令是把轄區裏最後一個陣列元素送給右鄰 PU, 同時自左鄰 PU 接受轄區外的前一個陣列元素 ITAG=110 ALL MPI_SENDREV (B(IEND), IOUNT, DATA_TYPE, R_NBR, ITAG, 1 B(ISTARTM1), IOUNT, DATA_TYPE, L_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) 引數 B(IEND) 是要送出去的資料起點 IOUNT 是要送出去的資料數量 DATA_TYPE 是要送出去的資料類別 R_NBR 是要送出去的目的地 PU id ( 右鄰 ) ITAG 是要送出去的資料標籤 B(ISTARTM1) 是要接收的陣列起點 IOUNT 是要接收的資料數量 DATA_TYPE 是要接收的資料類別 L_NBR 是要接收的來源地 PU id ( 左鄰 ) ITAG 是要接收的資料標籤 ISTATUS 是叫用這個副程式執行後的狀況 MPI_BAST MPI_BAST 是屬於 ' 集體通訊 ' 類副程式,BAST 是 Broadcast 的縮寫 當你要把同一項資料傳送給屬於同一個 communicator 其他各個 PU 時就可以叫用這一個副程式 既然是 ' 集體通訊 ' 類副程式, 參與平行計算的每一個 PU 都要叫用這一個副程式, 不允許只有少數幾個 PU 叫用它, 而其他 PU 不叫用它 MPI_BAST 的叫用格式如下 : 37

38 IROOT=0 ALL MPI_BAST ( B, IOUNT, DATA_TYPE, IROOT, & MPI_OMM_WORLD, IERR) 引數 B IOUNT DATA_TYPE IROOT 是要送出去的資料起點, 簡單變數或陣列名稱是要送出去的資料數量是要送出去的資料類別是要送出資料的 PU id MPI_BAST 的運作方式如圖 3.1 所示 : PU0 B B1 B2 B3 B4 PU0 B B1 B2 B3 B4 PU1 MPI_BAST PU1 B B1 B2 B3 B4 PU2 PU2 B B1 B2 B3 B4 PU3 PU3 B B1 B2 B3 B4 圖 3.1 MPI_BAST 示意圖 38

39 3.2 邊界資料交換的循序程式 T3SEQ T3SEQ 程式的 DO loop 裏, 在算出等號左邊的 A(I) 時, 要用到等號右邊的 (I) D(I) 和 B(I-1) B(I) B(I+1) 這在循序程式裏是奚鬆平常的事情, 但是在計算分割的平行程式裏卻是一件大事 因為要用到該 PU 轄區外 (outside territory) 的資料, 這時就牽涉到邊界資料交換 (boundary data exchange) 的問題 另外還要找出 A 陣列各個元素的最大值 AMAX PROGRAM T3SEQ Boundary Data Exchange Program - Sequential Version PARAMETER (NTOTAL=200) REAL*8 A(NTOTAL), B(NTOTAL), (NTOTAL), D(NTOTAL), AMAX OPEN(7,FILE='input.dat',STATUS='OLD',FORM='UNFORMATTED') READ (7) B READ (7) READ (7) D AMAX=-1.D12 DO I=2,NTOTAL-1 A(I)=(I)*D(I)+(B(I-1)+2.0*B(I)+B(I+1))*0.25 AMAX=MAX(AMAX,A(I)) WRITE(*,101) (A(I),I=1,NTOTAL,5) WRITE(*,102) AMAX 101 FORMAT(10F8.3) 102 FORMAT('MAXIMUM VALUE OF A ARRAY is', E15.5) STOP END 循序程式 T3SEQ 的測試結果如下 : MAXIMUM VALUE OF ARRAY A is.57500e+01 39

40 3.3 資料不切割的邊界資料交換平行程式 T3P 如何將循序程式 T3SEQ 平行化呢? 這一節先介紹計算切割而資料不切割的方法,3.4 節及 3.5 節再介紹計算及資料同時切割的方法 程式 T3P_1 也是利用副程式 STARTEND 來切割陣列, 算出各個 PU 轄區的 index 起迄點, 它是把陣列的第一段分給 PU0, 第二段分給 PU1, 餘類推 如圖 3.2 所示 : 左 left mpi_proc_null iend+1 iend1 iend ntotal cpu iend+1 istart2 iend1 istart iend ntotal cpu istart iend istart2 iend1 istart-1 cpu istart istart2 istart -1 mpi_proc_null 右 right is owned data is exchanged data 圖 3.2 資料不切割的邊界資料交換平行計算示意圖 圖 3.2 裏符號 代表轄區內的陣列元素, 符號. 代表轄區外的陣列元素, 符號 代表從左鄰 右舍傳送過來的陣列元素, 箭頭代表陣列元素傳送的方向 各個 PU 負責計算的範圍是從該 PU 的 istart 到 iend 程式 T3SEQ 的計算部份如下 : AMAX=-1.D12 DO I=2,NTOTAL-1 A(I)=(I)*D(I)+(B(I-1)+2.0*B(I)+B(I+1))*0.25 AMAX=MAX(AMAX,A(I)) 40

41 index I 是從 2 到 NTOTAL-1 在切割後, 只有 PU0 是從 2 開始, 其他的 PU 都是從 istart 開始, 所以必須設一個變數 istart2 來解決這個問題 : ISTART2=ISTART IF(MYID.EQ.0) ISTART2=2 而 loop 的終點只有最後一個 PU 是到 NTOTAL-1, 也就是 iend-1, 其他的 PU 都是到 iend, 所以必須再設一個變數 iend1 來解決這個問題 : IEND1=IEND IF(MYID.EQ.NPRO-1) IEND1=IEND-1 當 A(I) 的 I 等於 istart 時, 要用到 B(ISTART-1) 而當 A(I) 的 I 等於 iend 時, 要用到 B(IEND+1), 所以還必須設一個變數 istartm1 (istart minus 1) 和 iendp1 (iend plus 1) 來解決這 B 變數的 index 問題 : ISTARTM1=ISTART-1 IENDP1=IEND+1 在需要域外資料 (I-1 或 I+1 等 ) 的 DO loop 之前就要叫用 MPI_SENDREV 來取得該項資料 不過在這之前先要知道該 PU 的左鄰右舍是那一個 PU, 副程式 STARTEND 在切割陣列的 index 時是把第一段分給 PU0, 第二段分給 PU1, 第三段分給 PU2, 餘類推 所以一個 PU 的左鄰 L_NBR 就是該 PU 的 PU id 減一, 而其右鄰 R_NBR 就是該 PU 的 PU id 加一 只有第一個和最後一個 PU 是例外, 第一個 PU 沒有左鄰, 而最後一個 PU 沒有右鄰, 這時的左鄰右舍就要給他一個特定的名子叫做 MPI_PRO_NULL 這個 MPI_PRO_RULL 是在 mpif.h 檔裏已經設定的常數 L_NBR=MYID-1 R_NBR=MYID+1 IF(MYID.EQ.0) IF(MYID.EQ.NPRO-1) L_NBR=MPI_PRO_NULL R_NBR=MPI_PRO_NULL 現在來解決 B(I-1) 和 B(I+1) 的邊界資料交換問題, 這需要兩個 MPI_SENDREV, 一個解決 B(I-1) 的資料交換, 另外一個解決 B(I+1) 的資料交換 先來解決 B(I-1) 的邊界資料交換問題 從圖 3.2 中的 PU1 來看, 它要送 B(IEND) 給右鄰 " 當作右鄰的 B(ISTARTM1)", 又要自左鄰取得 " 左鄰的 B(IEND)" 做為它自己的 B(ISTARTM1) 如果要傳送的對象是 MPI_PRO_NULL 時, 是沒有傳送動作發生的 每一個 PU 都做同樣的動作, 就解決了 B(ISTARTM1) 的邊界資料交換問題 也就是 : 41

42 ITAG=110 ALL MPI_SENDREV (B(IEND), 1, MPI_REAL8, R_NBR, ITAG, 1 B(ISTARTM1), 1, MPI_REAL8, L_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) 再來解決 B(I+1) 的邊界資料交換問題 從圖 3.2 中的 PU1 來看, 它要送 B(ISTART) 給左鄰 " 當作左鄰的 B(IENDP1)", 又要自右鄰取得 " 右鄰的 B(ISTART)" 做為它自己的 B(IENDP1) 如果要傳送的對象是 MPI_PRO_NULL 時, 是沒有傳送動作發生的 每一個 PU 都做同樣的動作, 就解決了 B(IENDP1) 的邊界資料交換問題 也就是 : ITAG=120 ALL MPI_SENDREV (B(ISTART), 1, MPI_REAL8, L_NBR, ITAG, 1 B(IENDP1), 1, MPI_REAL8, R_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) 現在, 每一個 PU 執行的 DO loop 範圍是從 1 到 N, 所以算出來的 AMAX 是 A 陣列 NTOTAL 個元素裏 NP 分之一個陣列元素的最大值, 只是部份資料的最大值 所以此處叫用 MPI_ALLREDUE 找出各個 PU 轄區裏的最大值 AMAX 之中的最大值 GMAX (global maximum), 存放在每一個 PU 裏 ALL MPI_ALLREDUE ( AMAX, GMAX, 1, MPI_REAL8, MPI_MAX, 1 MPI_OMM_WORLD, IERR ) 至於甚麼時候使用 reduce, 甚麼時候使用 allreduce, 視需要而定 每一個 PU 都要用到 reduce 的結果時就要用 MPI_ALLREDUE, 只有一個 PU 需要用到 reduce 的結果時只要用 MPI_REDUE 就可以了, 因為 MPI_ALLREDUE 比較耗時間 所以完整的邊界資料交換平行程式如下 : 42

43 PROGRAM T3P Boundary data exchange without data partition Using MPI_SEND, MPI_REV to distribute input data PARAMETER (NTOTAL=200) INLUDE 'mpif.h' REAL*8 A(NTOTAL), B(NTOTAL), (NTOTAL), D(NTOTAL), AMAX, GMAX INTEGER NPRO, MYID, ISTATUS(MPI_STATUS_SIZE), ISTART, IEND, 1 L_NBR, R_NBR, 2 GSTART(0:31), GEND(0:31), GOUNT(0:31) ALL MPI_INIT (IERR) ALL MPI_OMM_SIZE (MPI_OMM_WORLD, NPRO, IERR) ALL MPI_OMM_RANK (MPI_OMM_WORLD, MYID, IERR) ALL STARTEND (NPRO, 1, NTOTAL, GSTART, GEND,GOUNT) ISTART=GSTART(MYID) IEND=GEND(MYID) PRINT *,' NPRO,MYID,ISTART,IEND=', NPRO,MYID,ISTART,IEND ISTARTM1=ISTART-1 IENDP1=IEND+1 ISTART2=ISTART IF(MYID.EQ.0) ISTART2=2 IEND1=IEND IF(MYID.EQ.NPRO-1) IEND1=IEND-1 L_NBR=MYID-1 R_NBR=MYID+1 IF(MYID.EQ.0) L_NBR=MPI_PRO_NULL IF(MYID.EQ.NPRO-1) R_NBR=MPI_PRO_NULL IF(MYID.EQ.0) THEN OPEN(7,FILE='input.dat',STATUS='OLD',FORM='UNFORMATTED') READ (7) B READ (7) READ (7) D 43

44 DO IDEST=1,NPRO-1 ISTART1=GSTART(IDEST) KOUNT1=GOUNT(IDEST) ITAG=10 ALL MPI_SEND ( B(ISTART1), KOUNT1, MPI_REAL8, IDEST, 1 ITAG, MPI_OMM_WORLD, IERR) ITAG=20 ALL MPI_SEND ((ISTART1), KOUNT1, MPI_REAL8, IDEST, 1 ITAG, MPI_OMM_WORLD, IERR) ITAG=30 ALL MPI_SEND (D(ISTART1), KOUNT1, MPI_REAL8, IDEST, 1 ITAG, MPI_OMM_WORLD, IERR) ELSE KOUNT=(IEND-ISTART+1) ISR=0 ITAG=10 ALL MPI_REV (B(ISTART), KOUNT, MPI_REAL8, ISR, ITAG, 1 MPI_OMM_WORLD, ISTATUS, IERR) ITAG=20 ALL MPI_REV ((ISTART), KOUNT, MPI_REAL8, SR, ITAG, 1 MPI_OMM_WORLD, ISTATUS, IERR) ITAG=30 ALL MPI_REV (D(ISTART), KOUNT, MPI_REAL8, ISR, ITAG, 1 MPI_OMM_WORLD, ISTATUS, IERR) ENDIF Exchange data outside the territory ITAG=110 ALL MPI_SENDREV (B(IEND), 1, MPI_REAL8, R_NBR, ITAG, 1 B(ISTARTM1), 1, MPI_REAL8, L_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) ITAG=120 ALL MPI_SENDREV (B(ISTART), 1, MPI_REAL8, L_NBR, ITAG, 1 B(IENDP1), 1, MPI_REAL8, R_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) 44

45 omputate, gather and write out the computed result AMAX=-1.0D12 DO I=2,NTOTAL-1 DO I=ISTART2,IEND1 A(I)=(I)*D(I)+(B(I-1)+2.0*B(I)+B(I+1))*0.25 AMAX=MAX(AMAX,A(I)) ITAG=110 IF(MYID.NE.0) THEN KOUNT=(IEND-ISTART+1) IDEST=0 ALL MPI_SEND (A(ISTART), KOUNT, MPI_REAL8, IDEST, ITAG, 1 MPI_OMM_WORLD, IERR) ELSE DO ISR=1,NPRO-1 ISTART1=GSTART(ISR) KOUNT1=GOUNT(ISR) ALL MPI_REV (A(ISTART1), KOUNT1, MPI_REAL8, ISR, ITAG, 1 MPI_OMM_WORLD, ISTATUS, IERR) ENDIF ALL MPI_ALLREDUE (AMAX, GMAX, 1, MPI_REAL8, MPI_MAX, 1 MPI_OMM_WORLD, IERR) IF(MYID.EQ.0) THEN WRITE(*,101) (A(I),I=1,NTOTAL,5) WRITE(*,102) GMAX ENDIF 101 FORMAT(10F8.3) 102 FORMAT('MAXIMUM VALUE of ARRAY A is', E15.5) ALL MPI_FINALIZE(IERR) STOP END 45

46 計算切割而資料不切割的邊界資料交換平行程式 T3P_1 的測試結果如下 : ATTENTION: nodes allocated by LoadLeveler, continuing... NPRO,MYID,ISTART,IEND= NPRO,MYID,ISTART,IEND= NPRO,MYID,ISTART,IEND= NPRO,MYID,ISTART,IEND= MAXIMUM VALUE OF ARRAY A is.57500e+01 在陣列不切割的平行程式裏,PU0 讀入資料之後也可以使用 MPI_BAST 把整個陣列傳送給各個 PU 這樣做程式寫起來比較簡單, 但是卻傳送了許多額外的資料 利敝得失端看兩種傳送方式所花費總時數的多寡而定, 通常一次傳送大量資料比分多次傳送小量資料來得省時 如何取得平行程式執行時的時間資訊將在第四章裏介紹 下列程式片段就是使用 MPI_BAST 取代 T3P 所使用的 MPI_SEND 和 MPI_REV 來分送輸入資料的結果 Read and distribute input data IF(MYID.EQ.0) THEN OPEN(7,FILE='input.dat',STATUS='OLD',FORM='UNFORMATTED') READ (7) B READ (7) READ (7) D ENDIF IROOT=0 ALL MPI_BAST ( B, NTOTAL, MPI_REAL8, IROOT, 1 MPI_OMM_WORLD, IERR) ALL MPI_BAST(, NTOTAL, MPI_REAL8, IROOT, 1 MPI_OMM_WORLD, IERR) ALL MPI_BAST( D, NTOTAL, MPI_REAL8, IROOT, 1 MPI_OMM_WORLD, IERR) 46

47 3.4 資料切割的邊界資料交換平行程式 ( 一 ) T3DP_1 計算及資料同時切割時, 如果是在 NP 個 PU 上平行計算, 則陣列的長度 N 只需要原來長度 NTOTAL 的 NP 分之一 此時,NTOTAL 必須能被 NP 整除 再加上前後各保留一個界外陣列元素的存放位置, 則所需的陣列長度為 N+2, 本文主張其 DIMENSION 採用 (0:N+1) 的方式來撰寫, 則其轄區內資料是存放在 index 1 到 N 於是陣列長度的設定為 : REAL*8 A(0:N+1), B(0:N+1), (0:N+1), D(0:N+1), T(NTOTAL), AMAX, GMAX 如圖 3.3 所示 : left cpu0 cpu1 cpu2 cpu3 mpi_proc_null n + index n 1 n + index n 1 n + index n 1 n + index n 1 right is owned data is exchanged data mpi_proc_null 圖 3.3 資料切割且交換一個邊界資料的平行計算示意圖 這樣一來, 每一個 PU 的 DO loop 起點為 1, 終點為 N, 於是 : ISTART=1 IEND=N 第一個 PU 的 DO loop 起點為 2, 其他 PU 的 DO loop 起點為 1, 最後一個 PU 的 DO loop 終點為 N-1, 其他 PU 的 DO loop 終點為 N, 也就是 : 47

48 ISTART2=ISTART IF(MYID.EQ.0) ISTART2=2 IEND1=IEND IF(MYID.EQ.NPRO-1) IEND1=IEND-1 每一個 PU 要把轄區內最後一個陣列元素送給右邊的 PU, 這個陣列元素的位置是 iend, 同時接受來自左邊 PU 的陣列元素是放在 istart-1 所以 : ISTARTM1=ISTART-1 ITAG=110 ALL MPI_SENDREV (B(IEND), 1, MPI_REAL8, R_NBR, ITAG, 1 B(ISTARTM1), 1, MPI_REAL8, L_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) 每一個 PU 也要把轄區內第一個陣列元素送給左邊的 PU, 這個陣列元素的位置是 istart, 同時接受來自右邊 PU 的一個陣列元素是放在 iend+1 所以 : IENDP1=IEND+1 ITAG=120 ALL MPI_SENDREV (B(ISTART), 1, MPI_REAL8, L_NBR, ITAG, 1 B(IENDP1), 1, MPI_REAL8, R_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) NTOTAL 個元素的 B D 陣列資料讀入陣列 T 之後, 在用 MPI_SATTER 分送給各個 PU 時,B D 陣列的 DIMENSION 是從零開始, 而輸入資料是從 1 開始存放 所以在 MPI_SATTER 裏必需敘明要從 B(1) (1) D(1) 開始存放 MPI_GATHER 時亦然 IROOT=0 ALL MPI_SATTER (T, N, MPI_REAL8, B(1), N, MPI_REAL8, 1 IROOT, MPI_OMM_WORLD, MPI_ERR) ALL MPI_GATHER (A(1), N, MPI_REAL8, T, N, MPI_REAL8, 1 IROOT, MPI_OMM_WORLD, MPI_ERR) 計算及資料都切割的邊界資料交換平行程式 T3DP_1 如下 : 48

49 PROGRAM T3DP_1 Data Partition Using MPI_SATTER, MPI_GATHER with -1,+1 data exchange PARAMETER (NTOTAL=200, NP=4, N=NTOTAL/NP) INLUDE 'mpif.h' REAL*8 A(0:N+1), B(0:N+1), (0:N+1), D(0:N+1), T(NTOTAL), & AMAX, GMAX INTEGER NPRO, MYID, ISTATUS(MPI_STATUS_SIZE), ISTART, IEND, 1 L_NBR, R_NBR ALL MPI_INIT (IERR) ALL MPI_OMM_SIZE (MPI_OMM_WORLD, NPRO, IERR) ALL MPI_OMM_RANK (MPI_OMM_WORLD, MYID, IERR) IF (NP.NE.NPRO) THEN PRINT *,' NP NOT EQUAL TO NPRO, PROGRAM WILL STOP' PRINT *,' NP,NPRO=', NP,NPRO STOP ENDIF ISTART=1 IEND=N ISTARTM1=ISTART-1 IENDP1=IEND+1 ISTART2=1 IF(MYID.EQ.0) ISTART2=2 IEND1=N IF(MYID.EQ.NPRO-1) IEND1=N-1 L_NBR=MYID-1 R_NBR=MYID+1 IF(MYID.EQ.0) L_NBR=MPI_PRO_NULL IF(MYID.EQ.NPRO-1) R_NBR=MPI_PRO_NULL Read & distribute input data IF(MYID.EQ.0) THEN 49

50 OPEN(7,FILE='input.dat',STATUS='OLD',FORM='UNFORMATTED') READ (7) T! read array B ENDIF IROOT=0 ALL MPI_SATTER (T, N, MPI_REAL8, B(1), N, MPI_REAL8, 1 IROOT, MPI_OMM_WORLD, IERR) IF(MYID.EQ.0) THEN READ (7) T! read array ENDIF ALL MPI_SATTER (T, N, MPI_REAL8, (1), N, MPI_REAL8, 1 IROOT, MPI_OMM_WORLD, IERR) IF(MYID.EQ.0) THEN READ (7) T! read array D ENDIF ALL MPI_SATTER (T, N, MPI_REAL8, D(1), N, MPI_REAL8, 1 IROOT, MPI_OMM_WORLD, IERR) Boundary data exchange ITAG=110 ALL MPI_SENDREV (B(IEND), 1, MPI_REAL8, R_NBR, ITAG, 1 B(ISTARTM1), 1, MPI_REAL8, L_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) ITAG=120 ALL MPI_SENDREV (B(ISTART), 1, MPI_REAL8, L_NBR, ITAG, 1 B(IENDP1), 1, MPI_REAL8, R_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) ompute and output result AMAX=-1.0D12 DO I=2,NTOTAL-1 DO I=ISTART2,IEND1 A(I)=(I)*D(I)+(B(I-1)+2.0*B(I)+B(I+1))*0.25 AMAX=MAX(AMAX,A(I)) ALL MPI_GATHER (A(1), N, MPI_REAL8, T, N, MPI_REAL8, 50

51 1 IROOT, MPI_OMM_WORLD, IERR) ALL MPI_ALLREDUE (AMAX, GMAX, 1, MPI_REAL8, MPI_MAX, 1 MPI_OMM_WORLD, IERR) IF(MYID.EQ.0) THEN WRITE(*,101) (T(I),I=1,NTOTAL,5) WRITE(*,102) GMAX ENDIF 101 FORMAT(10F8.3) 102 FORMAT('MAXIMUM VALUE OF ARRAY A is', E15.5) ALL MPI_FINALIZE (MPI_ERR) STOP END 計算及資料都切割的邊界資料交換平行程式 T3DP_1 的測試結果如下 : ATTENTION: nodes allocated by LoadLeveler, continuing MAXIMUM VALUE OF ARRAY A is.57500e+01 51

52 3.5 資料切割的邊界資料交換平行程式 ( 二 ) T3DP_2 如果 3.4 節裏程式 T3DP_1 的 DO loop 需要用到兩個鄰近陣列元素時, 如下列所示 : DO I=3,NTOTAL-2 A(I)=(I)*D(I)+( B(I-2)+2.0*B(I-1)+2.0*B(I) & +2.0*B(I+1)+B(I+2) )*0.125 其平行化方法與程式 T3DP_1 相同 只要把切割後陣列的 DIMENSION 前後各加兩個陣列 元素, 把 DIMENSION 改為 (-1:N+2), 其轄區內資料仍然是存放在 index 1 到 N & REAL*8 A(-1:N+2), B(-1:N+2), (-1:N+2), D(-1:N+2), T(NTOTAL), AMAX, GMAX ISTART=1 IEND=N 如圖 3.4 所示 : LEFT cpu0 cpu1 cpu2 cpu3 mpi_proc_null n n n 1 2 n n + + index n 1 2 n n + + index n 1 2 n n + + index n 1 2 is exchanged data mpi_proc_null RIGHT is owned data 圖 3.4 資料切割且交換兩個邊界資料的平行計算示意圖 52

53 DO loop 的起迄 index 也要加以修改 第一個 PU 的起點是 3, 最後一個 PU 的終點是 NTOTAL-2, 所以 : ISTART3=1 IF(MYID.EQ.0) ISTART3=3 IEND2=IEND IF(MYID.EQ.NPRO-1) IEND2=IEND-2 當然, 邊界資料的交換量也要由一個改為兩個 每一個 PU 要把轄區內最後兩個陣列元素送 給右邊的 PU, 這兩個陣列元素的起點是 iend-1, 同時要接受來自左邊 PU 的兩個陣列元素, 從 istart-2 放起 所以 : IENDM1=IEND-1 ISTARTM2=ISTART-2 ITAG=110 ALL MPI_SENDREV (B(IENDM1), 2, MPI_REAL8, R_NBR, ITAG, 1 B(ISTARTM2), 2, MPI_REAL8, L_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) 每一個 PU 也要把轄區內最前面兩個陣列元素送給左邊的 PU, 這兩個陣列元素的起點是 istart, 同時要接受來自右邊 PU 的兩個陣列元素, 從 iend+1 放起 所以 : IENDP1=IEND+1 ITAG=120 ALL MPI_SENDREV (B(ISTART), 2, MPI_REAL8, L_NBR, ITAG, 1 B(IENDP1), 2, MPI_REAL8, R_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) 交換邊界兩個陣列元素的平行程式如下 : 53

54 PROGRAM T3DP_2 Balanced Data Partition Using MPI_SATTER, MPI_GATHER with -1,-2,+1,+2 data exchange PARAMETER (NTOTAL=200,NP=4,N=NTOTAL/NP) INLUDE 'mpif.h' REAL*8 A(-1:N+2), B(-1:N+2), (-1:N+2), D(-1:N+2), T(NTOTAL), & AMAX, GMAX INTEGER NPRO, MYID, ISTATUS(MPI_STATUS_SIZE), ISTART, IEND, 1 L_NBR, R_NBR ALL MPI_INIT (IERR) ALL MPI_OMM_SIZE (MPI_OMM_WORLD, NPRO, IERR) ALL MPI_OMM_RANK (MPI_OMM_WORLD, MYID, IERR) IF (NP.NE.NPRO) THEN PRINT *,' NP NOT EQUAL TO NPRO, PROGRAM WILL STOP' PRINT *,' NP,NPRO=', NP,NPRO STOP ENDIF ISTART=1 ISTARTM2=ISTART-2 ISTART3=1 IF(MYID.EQ.0) ISTART3=3 IEND=N IEND2=IEND IF(MYID.EQ.NPRO-1) IEND2=IEND-2 IENDM1=IEND-1 IENDP1=IEND+1 IEND1=N IF(MYID.EQ.NPRO-1) IEND1=N-1 L_NBR=MYID-1 R_NBR=MYID+1 IF(MYID.EQ.0) L_NBR=MPI_PRO_NULL IF(MYID.EQ.NPRO-1) R_NBR=MPI_PRO_NULL 54

55 Read & distribute input data IF(MYID.EQ.0) THEN OPEN(7,FILE='input.dat',STATUS='OLD',FORM='UNFORMATTED') READ (7) T! read array B ENDIF IROOT=0 ALL MPI_SATTER (T, N, MPI_REAL8, B(1), N, MPI_REAL8, 1 IROOT, MPI_OMM_WORLD, IERR) IF(MYID.EQ.0) THEN READ (7) T! read array ENDIF ALL MPI_SATTER (T, N, MPI_REAL8, (1), N, MPI_REAL8, 1 IROOT, MPI_OMM_WORLD, IERR) IF(MYID.EQ.0) THEN READ (7) T! read array D ENDIF ALL MPI_SATTER (T, N, MPI_REAL8, D(1), N, MPI_REAL8, 1 IROOT, MPI_OMM_WORLD, IERR) Boundary data exchange ITAG=110 ALL MPI_SENDREV (B(IENDM1), 2, MPI_REAL8, R_NBR, ITAG, 1 B(ISTARTM2), 2, MPI_REAL8, L_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) ITAG=120 ALL MPI_SENDREV (B(ISTART), 2, MPI_REAL8, L_NBR, ITAG, 1 B(IENDP1), 2, MPI_REAL8, R_NBR, ITAG, 2 MPI_OMM_WORLD, ISTATUS, IERR) ompute and output result AMAX=-1.0D12 DO I=3,NTOTAL-2 DO I=ISTART3,IEND2 A(I)=(I)*D(I)+( B(I-2)+2.0*B(I-1)+2.0*B(I) 55

56 & +2.0*B(I+1)+B(I+2) )*0.125 AMAX=MAX(AMAX,A(I)) IDEST=0 ALL MPI_GATHER (A(1), N, MPI_REAL8, T, N, MPI_REAL8, 1 IDEST, MPI_OMM_WORLD, IERR) ALL MPI_ALLREDUE (AMAX, GMAX, 1, MPI_REAL8, MPI_MAX, 1 MPI_OMM_WORLD, IERR) AMAX=GMAX IF(MYID.EQ.0) THEN WRITE(*,101) (T(I),I=1,NTOTAL,5) WRITE(*,102) AMAX ENDIF 101 FORMAT(10F8.3) 102 FORMAT('MAXIMUM VALUE OF ARRAY A is ', E15.5) ALL MPI_FINALIZE (IERR) STOP END 交換邊界兩個陣列元素的平行程式 T3DP_2 的測試結果如下 : ATTENTION: nodes allocated by LoadLeveler, continuing MAXIMUM VALUE OF ARRAY A is.44847e+01 56

57 第四章格點數不能整除的平行程式 本章將探討格點數 (grid points) 不能被參與平行計算的 PU 數目整除時的處理方法 所謂格點數就是陣列的 dimension 4.1 節循序程式 T4SEQ 裏陣列的 dimension 是 161, 只能夠被 7 和 23 整除 4.2 節介紹 MPI_SATTERV 和 MPI_GATHERV 兩個 ' 集體通訊 ' 類副程式, 其功能與 MPI_SATTER 和 MPI_GATHER 相似, 但是分送和來自各個 PU 的資料不必等長 4.3 節介紹 MPI_PAK 和 MPI_UNPAK 兩個資料集結與分解指令, 以及同步指令 MPI_BARRIER 和取得時鐘時刻指令 MPI_WTIME 4.4 節說明使用這些 MPI 指令來把循序程式 T4SEQ 改寫成平行程式 T4DP 57

58 4.1 格點數不能整除的循序程式 T4SEQ 循序程式 T4SEQ 裏陣列的 dimension 是 161, 只能夠被 7 和 23 整除 該程式除了陣列資料 A B D 之外, 還用到三個純量資料 (scalar data) P Q R 在設定各個變數的初值 (initial value) 之後, 也寫入磁檔, 便利平行化時的驗證 PROGRAM T4SEQ odd-dimensioned array with -1, +1 data reference PARAMETER (NTOTAL=161) REAL*8 A(NTOTAL), B(NTOTAL), (NTOTAL), D(NTOTAL), P, Q, R DATA P, Q, R/1.45, 2.62, 0.5/ Data generation DO I=1,NTOTAL B(I)=3.D0/DFLOAT(I)+1.D0 (I)=2.D0/DFLOAT(I)+1.D0 D(I)=1.D0/DFLOAT(I)+1.D0 OPEN(1,FILE='input.dat',FORM='UNFORMATTED') WRITE(1) B WRITE(1) WRITE(1) D WRITE(1) P,Q,R LOSE(1) Read input data & distribute input data OPEN(1,FILE='input.dat',STATUS='OLD',FORM='UNFORMATTED') READ (1) B READ (1) READ (1) D READ (1) P,Q,R ompute & output computed result 58

59 DO I=2,NTOTAL-1 A(I)=(I)*D(I)*P+(B(I-1)+2.0*B(I)+B(I+1))*Q+R WRITE(*,101) (A(I),I=1,NTOTAL,5) 101 FORMAT(10F8.3) STOP END 格點數不能整除的循序程式 T4SEQ 的測試結果如下 :

60 4.2. MPI_SATTERV MPI_GATHERV 如果要把 4.1 節循序程式 T4SEQ 平行化, 而且陣列也要切割時, 其 dimension 是 161, 只能夠被 7 和 23 整除 不能夠被 等參與平行計算的 PU 數目整除, 因此在 個 PU 上平行計算時, 就不能夠使用 MPI_SATTER 來分送輸入資料, 及使用 MPI_GATHER 來收集各個 PU 上的資料 當然, 此時還是可以使用 MPI_SEND 和 MPI_REV 來完成陣列資料的分送與與收集 不過,MPI 還備有 MPI_SATTERV 和 MPI_GATHERV 可以用來分送與與收集不等長陣列資料 MPI_SATTERV 和 MPI_SATTER 的功能相似,MPI_GATHERV 和 MPI_GATHER 的功能相似, 但是 MPI_SATTER 是用來分送等量資料給每一個 PU,MPI_GATHER 是收集來自每一個 PU 的等量資料,MPI_SATTERV 和 MPI_GATHERV 則不受這種限制 MPI_SATTERV 的叫用格式如下 : ALL MPI_SATTERV( T, GOUNT, GDISP, MPI_REAL8, 1 (1), MYOUNT, MPI_REAL8, 2 IROOT, MPI_OMM_WORLD, IERR) MPI_SATTERV 是 IROOT PU 把一個陣列 T 依 PU id 的順序分段分送給每一個 PU, 包括 IROOT PU 在內 由於分送給各個 PU 的資料量可以不一樣多, 因此必須使用一個陣列 GOUNT 來存放分送給各個 PU 的資料數量, 再另外使用一個陣列 GDISP 來存放送出資料在 T 陣列上的相對位置 其引數依序為 : T GOUNT GDISP MPI_REAL8 是待送出的資料陣列起點整數陣列, 是存放要送給各個 PU 的資料數量整數陣列, 是存放要送給各個 PU 資料起點在 T 陣列上的相對位置是待送出資料的類別 (1) MYKOUNT MPI_REAL8 IROOT 是接收資料存放的起點是接收資料的數量是接收資料的類別是送出資料陣的 PU id 此處 GOUNT 及 GDISP 是叫用 STARTEND 副程式算出來各個 PU 上經切割後陣列片段的長度及各該陣列片段在未切割前的起點 START index 再減一 為了配合 PU id 是從零起算, GOUNT GDISP 的 dimension 也是從零開始 所以 : REAL*8 A(0:N+1), B(0:N+1), (0:N+1), D(0:N+1), T(NTOTAL) 60

61 INTEGER NPRO, MYID, ISTATUS(MPI_STATUS_SIZE), ISTART, IEND, 1 L_NBR, R_NBR, 2 GSTART(0:31), GEND(0:31), GOUNT(0:31), GDISP(0:31) ALL STARTEND (NPRO, 1, NTOTAL, GSTART, GEND, GOUNT) DO I=0,NPRO-1 GDISP(I)=GSTART(I)-1 MYOUNT= GOUNT (MYID) 另外, 分割後的陣列 A B D 的 dimension 是從零開始, 而讀入資料是從一開始存放, 所以叫用 MPI_SATTERV 時要標明從一開始 B(1) (1) D(1) 還有 T 陣列上的相對位置 (displacement 或 offset) 是從零起算而不是從一起算, 所以其值由 STARTEND 副程式算出來的 GSTART 要再減掉一 MPI_GATHERV 的叫用格式如下 : ALL MPI_GATHERV (A(1), MYKOUNT, MPI_REAL8, 1 T, GOUNT, GDISP, MPI_REAL8, 2 IROOT, MPI_OMM_WORLD, IERR) MPI_GATHERV 與 MPI_SATTERV 的動作剛好相反, 是 IROOT PU 收集每一個 PU ( 包括 IROOT PU) 送給它的陣列 A 依 PU id 的順序存入陣列 T 裏頭 由於來自各個 PU 的資料量不一定一樣多, 因此必須用一個陣列 GOUNT 來存放來自各個 PU 的資料數量, 用另外一個陣列 GDISP 來存放要存入陣列 T 上的相對位置 其引數依序為 : A(1) MYKOUNT MPI_REAL8 是待送出的資料陣列起點 是待送出資料的數量 是待送出資料的類別 T GOUNT GDISP MPI_REAL8 IROOT 是接收資料存放的陣列位址整數陣列, 是存放來自各個 PU 的資料數量整數陣列, 是存放來自各個 PU 資料要存入 T 陣列的相對位置是接收資料的類別是接收資料的 PU id 61

62 4.3 MPI_PAK UNPAK BARRIER WTIME 循序程式 T4SEQ 裏除了陣列資料 A B D 之外, 還用到三個純量資料 P Q R 在陣列資料切割與平行化之後, 可以使用 MPI_SATTERV 來分送輸入之陣列資料, 使用 MPI_GATHERV 來收集由各個 PU 計算出來的結果 至於三個純量資料 P Q R 在讀入之後, 可以使用三個 MPI_BAST 分別傳送給各個 PU 不過, 到現在為止,PU 之間的資料傳送速度, 比記憶體內資料移動的速度要慢上很多倍 所以在 PU 之間一次傳送大量資料, 比分多次傳送而每次只傳送少量資料來得有效率 但是只有陣列資料才可以一次傳送多個資料, 還好 MPI 備有資料集結副程式 MPI_PAK 可以把非連續位址資料 (noncontiguous data) 移入記憶體的連續位址裏 (contiguous memory locations), 俗稱緩衝區 (buffer area), 緩衝區為一個字符陣列 (character array) 然後把該字符陣列傳送給對方, 對方收到該字符陣列後, 再依存入的順序使用 MPI_UNPAK 取出各個變數的數值 另外, 還可以叫用 ' 同步 ' 副程式 MPI_BARRIER 使各個 PU 達到同步的效果 叫用 MPI 涵數 (Fortran function) MPI_WTIME 來取得電腦時鐘時刻 (wall clock time) MPI_PAK MPI_UNPAK MPI_PAK 是要把非連續位址資料, 像循序程式 T4SEQ 裏數個純量資料 P Q R 或數個陣列片段移入一個字符陣列裏 這些資料可以是同一種資料類別, 也可以是不同的資料類別 送出資料的 PU 把集結 (pack) 好了的字符陣列傳送給對方, 對方收到該字符陣列後, 再依存入的順序使用 MPI_UNPAK 取出各個變數的數值 要使用 MPI_PAK MPI_UNPAK, 首先要在該程式裏設定一個字符陣列 例如要存放三 個倍準數 P Q R 的緩衝區, 一個倍準數需要 8 個字符, 三個倍準就需要 24 個字符 所以 設定一個字符陣列 BUF1 其長度為 24 個字符 如下 : INTEGER BUFSIZE PARAMETER (BUFSIZE=24) HARATER BUF1(BUFSIZE) MPI_PAK 的叫用格式如下 : & ALL MPI_PAK (P, 1, MPI_REAL8, BUF1, BUFSIZE, IPOS, MPI_OMM_WORLD, IERR) 62

63 其中引數 P 是要移入 BUF1 的變數起點 1 是要移入 BUF1 的資料數量 MPI_REAL8 是要移入 BUF1 的資料類別 BUF1 BUFSIZE IPOS 是存放集結資料的字符陣列起點 是字符陣列 BUF1 的字符長度 是要存入 BUF1 的字符位置, 從零開始, 自動更新 IPOS 在移入第一筆資料之前要設為零, 每次移入一筆資料, 其值會自動更新, 指向下一個 空位 所以 PU0 讀入 P Q R 三個倍準數之後, 再把他們移入 BUF1 可以寫成 : IF (MYID.EQ.0) THEN READ(7) P,Q,R IPOS = 0 ALL MPI_PAK (P, 1, MPI_REAL8, BUF1, BUFSIZE, IPOS, 1 MPI_OMM_WORLD, IERR) ALL MPI_PAK (Q, 1, MPI_REAL8, BUF1, BUFSIZE, IPOS, 1 MPI_OMM_WORLD, IERR) ALL MPI_PAK (R, 1, MPI_REAL8, BUF1, BUFSIZE, IPOS, 1 MPI_OMM_WORLD, IERR) ENDIF 然後可以使用 MPI_BAST 將字符陣列 BUF1 傳送給其他每一個 PU 如下 : IROOT=0 ALL MPI_BAST (BUF1, BUFSIZE, MPI_HARATER, IROOT, & MPI_OMM_WORLD, IERR) MPI_UNPAK 的叫用格式如下 : & ALL MPI_UNPAK (BUF1, BUFSIZE, IPOS, P, 1, MPI_REAL8, MPI_OMM_WORLD, IERR) 其中引數 BUF1 BUFSIZE IPOS 是存放集結資料的字符陣列起點 是字符陣列 BUF1 的字符長度 是要從 BUF1 取出資料的字符位置, 從零開始, 自動更新 63

64 P 是自 BUF1 取出的變數位址 1 是自 BUF1 取出的資料數量 MPI_REAL8 是自 BUF1 取出的資料類別 IPOS 在取出第一筆資料之前要設為零, 每次取出一筆資料之後, 其值會自動更新, 指向下一 筆資料的位置 所以各個 PU 在執行過 bcast BUF1 之後, 要自 BUF1 取出 P Q R 三個倍 準數可以寫成 : IF (MYID.NE.0) THEN IPOS=0 ALL MPI_UNPAK (BUF1, BUFSIZE, IPOS, P, 1, MPI_REAL8, & MPI_OMM_WORLD, IERR) ALL MPI_UNPAK (BUF1, BUFSIZE, IPOS, Q, 1, MPI_REAL8, & MPI_OMM_WORLD, IERR) ALL MPI_UNPAK (BUF1, BUFSIZE, IPOS, R, 1, MPI_REAL8, & MPI_OMM_WORLD, IERR) ENDIF 如果要 PAK 在一起的資料是同一種資料類別時, 可以設定一個陣列, 將這些資料逐一移入 該陣列, 然後將該陣列傳送給對方 對方收到該陣列後, 再依序自該陣列移出各個資料 REAL*8 BUF1(3) IF (MYID.EQ.0) THEN READ(7) P,Q,R BUF1(1)=P BUF1(2)=Q BUF1(3)=R ENDIF IROOT=0 ALL MPI_BAST (BUF1, 3, MPI_REAL8, IROOT, & MPI_OMM_WORLD, IERR) IF(MYID.NE.0) THEN P=BUF1(1) Q=BUF1(2) R=BUF1(3) ENDIF 64