本論文探討利用量子點的結構位能為三維拋物線位能配合
k p
單能帶和有 效質量近似法計算激子的庫侖交互作用，經過二次量子化後的 Hamiltonian 且選 擇一組無交互作用的激子系統組態為基底，可展開為庫侖、交換能矩陣。而原本由於原子排列方式為Δr=0.38 nm，因此選擇了Δr=0.38 nm 來分析，而分析後結 果發現量子點由小到大數值解與解析解誤差也在可接受範圍內，因此將來做數值
73
受範圍內，此結果可知威格納－塞茲晶胞內的積分對於整體的積分貢獻很小，因 此 a0=0.38 nm 的結果是可信的。因此我們可選擇Δr=0.38 nm 以及 a0=0.38 nm 來 計算直接庫侖矩陣與交換能矩陣，可使我們更快、更有效率的計算出結果。
74
附錄 A、CUDA compiler 方法
Fortran 程式碼 compiler 方式會根據不同開發者而有差異，表 3.3.1 是使用 CentOS 5.4 作業系統在終端機介面 compiler 的方法：
表 A.1 不同廠商在 CentOS 5.4 一般 compiler 方式。
開發廠商 一般 compiler 方法
GNU gfortran 「檔案名稱.f」
Intel ifort 「檔案名稱.f」
The Portland Group(PGI) pgf95 「檔案名稱.f」
目前 Fortran 只有 PGI 公司開發的 PGI Fortran 支援 CUDA 資料庫，因此我們將 使用 PGI Fortran 來撰寫程式碼，以下為 PGI Fortran 使用 CUDA 資料庫 compiler 程式碼的方法：
pgf95 「檔案名稱.f」 -MCUDA=4.0 或 pgf95 「檔案名稱.cuf」
加入 CUDA 資料庫的 Fortran 程式又可稱為 CUDA Fortran。
75
附錄 B-1、CUDA 基本程式編寫
由於將一般的寫法修改成 GPU 的寫法細節有點煩雜，因此先由最簡單的數 值寫入參數中來分析，圖 B-1.1 為程式流程圖，而圖 B-1.2 為將陣列中的四個元 素寫入數值的程式碼，我們將每個元素分別利用 GPU 的四個核心來分別寫入。
單核 CPU GPU
圖 B-1.1 使用 CPU 與 GPU 的程式流程圖。
76
單核 CPU GPU[Host 端]
單核 CPU GPU[Device 端(kernel)]
單核 CPU[執行結果] GPU[執行結果]
圖 B-1.2 一般 Fortran 修改成 CUDA Fortran 簡單範例：紅色框框表示需要修改或增加的部分。
由於 CUDA 是加入 Fortran 的資料庫，表示在程式中會新增幾個指令，而一 般 Fortran 指令還是可以繼續沿用。表 B.1 是針對 Fortran 使用 CUDA 資料庫需 要增加或修改程式碼：
77
表 B-1.1 Fortran 使用 CUDA 需要修改或增加的部分。
Host 端 Device 端
1. 宣告(常數、變數) 1. 使用 MODULE 指令及定義 CUDA 副程式 2. Host 端與 Device 端記憶體互傳 2. 宣告(常數、變數)
3. 呼叫 kernel(CUDA 副程式) 3. 使用新增指令將 thread 與 block 編號 4. 釋放 device 記憶體 4. 利用平行運算的概念做計算
接下來要解釋 CUDA Fortran 程式碼要如何編輯。
78
附錄 B-2、Host 端撰寫方式
在 Host 端執行程式流程大致上分成 5 個步驟：(如圖 B-1.1 所示)
1. Host 端配置 device 記憶體 2. 上傳資料到 device 記憶體 3. 呼叫 kernel 運算
4. 下載資料回 host 記憶體 5. 釋放 device 記憶體
GPU[Host 端]
圖 B-2.1 在圖 B-1.1 中 Host 端程式碼進行分析。
將圖 B-2.1 程式碼分成 4 個部分來說明，分別是宣告方式、傳送方式、呼叫 kernel 與釋放 device 記憶體。
一、宣告方式
在 CUDA Fortran 宣告有分 2 個部分，一個是由 Host 端配置 device 記憶體(如 圖 B-2.2 的 Host 端變數「numout_d(4)」與常數「conin」)，另一個是在 Device
79
端宣告接收 Host 端的參數(此部分下一個小節會介紹)，表 B-2.1 要先介紹 CUDA 在 Host 端宣告方式：
表 B-2.1 CUDA 在 Host 端宣告方式。
變數 常數
整數 INTEGER,DEVICE INTEGER
浮點數(單精準度) REAL, DEVICE REAL
浮點數(雙精準度) DOUBLE PRECISION, DEVICE DOUBLE PRECISION 複數(單精準度) COMPLEX, DEVICE COMPLEX
複數(雙精準度) DOUBLE COMPLEX, DEVICE DOUBLE COMPLEX
邏輯 LOGICAL, DEVICE LOGICAL
字串 CHARACTER*1, DEVICE CHARACTER *1, DEVICE
這裡要注意的是字串宣告，在 CUDA 裡最多只能定義 1 個文字的字串，而它的
「numout_d=numout」表示由 CPU 讀取或計算完的資料「numout」上傳到 Device 記憶體「numout_d」，「numout=numout_d」則表示由 GPU 計算完的資料「numout_d」
下載回 Host 端記憶體「numout」。
80 GPU[Host 端]
:
圖 B-2.4 在圖 B-2.1 中 Host 端與 Device 端參數傳送方式。
三、呼叫 kernel
由 2.3.1 小節可得知，CUDA 可分為 grid、block、thread，呼叫 kernel 意思 就是創造 grid 使用 GPU 做計算，以下是在程式的編寫方式：
CALL 「副程式檔名」<<<「block 個數」,「thread 個數」>>>(「輸入輸出參數」)
範例如圖 B-2.5 所示，因此我們可以設定 block 個數與 thread 個數。
GPU[Host 端]
圖 B-2.5 在圖 B-2.1 中呼叫 kernel。
由於 CUDA 是三維結構，所以呼叫三維 CUDA 可寫成以下形式：
CALL 「副程式檔名」<<<DIM3(「block_x 個數」,「block_y 個數」,「block_z 個數」),DIM3(「thread_x 個數」,「thread_y 個數」,「thread_z 個數」)>>>(「輸入輸出的參數」)
呼叫三維 CUDA 要利用 DIM3 函式，必須使用 CUDAFOR 資料庫，如圖 B-2.6。
:
圖 B-2.6 呼叫三維 CUDA 寫法：「USA CUDAFOR」表示使用 CUDAFOR 資料庫。
四、釋放 device 記憶體
釋 放 device 記 憶 體 必 須 利 用 CUDAFREE 函 數 來 釋 放 ， 因 此 要 使 用 CUDAFOR 資料庫，如圖 B-2.7 所示。「USE CUDAFOR」表示使用 CUDAFOR 資料庫，而釋放 device 記憶體撰寫方式為：
「Host 端宣告的參數」=CUDAFREE(「Device 端要釋放記憶體的參數」)
圖 B-2.7 中的「istat」為 Host 端宣告的任意參數即可，「numout_d」為需要釋放
81
的 device 記憶體參數。
GPU[Host 端]
:
圖 B-2.7 在圖 B-2.1 中宣告 device 記憶體利用 CUDAFREE 函數釋放。
82
附錄 B-3、Device 端撰寫方式
Device 端執行程式流程大致上分成 4 個步驟：(如圖 B-3.1 所示)
1. 利用 MODULE 將整個 kernel 包裝成 GPU 可運行的格式 2. Device 端宣告接收 Host 端的參數
3. 將 block 與 thread 編號 4. 平行運算
GPU[Device 端(kernel)]
圖 B-3.1 在圖 B-1.1 中 Device 端程式碼進行分析。
將圖 B-3.1 分成 4 個部分說明，分別是包裝方式、宣告方式、CUDA 編號、計算 想法。
一、包裝方式
使用 GPU 計算必須將 Device 端的程式碼利用 MODULE 包裝起來，compiler
後會產生 檔案，之後在 Host 端利用「USE test_gpu_m」來使用此 mod，
如圖 B-3.3；定義 CUDA 副程式與一般 Fortran 副程式不同的部分是需要在
「SUBROUTINE」前面加入「ATTRIBUTES(GLOBAL)」指令，才能使 CUDA 副程式供 GPU 運算，如圖 B-3.2。
83
整數 INTEGER INTEGER,VALUE
浮點數(單精準度) REAL REAL,VALUE
浮點數(雙精準度) DOUBLE PRECISION DOUBLE PRECISION,VALUE
複數(單精準度) COMPLEX COMPLEX,VALUE
複數(雙精準度) DOUBLE COMPLEX DOUBLE COMPLEX,VALUE
邏輯 LOGICAL LOGICAL,VALUE
字串 CHARACTER *1 CHARACTER *1
字串宣告在 Device 端的常數宣告方式與其他常數宣告方式不同。
84
85
圖 B-3.5 此圖說明 CUDA 編號指令範例。我們選擇 Grid 內有 5 個 block，每個 block 有 6 個 thread，
當呼叫 CUDA 副程式時，block「BLOCKIDX」會自動編號 1~5，thread「THREADIDX」會自動 編號 1~6；grid 維度「GRIDDIM」=5，block 維度「BLOCKDIM」=6。
我們可自定編號將 Grid 內所有 thread 編成不同的號碼，這樣可以更明確決 定每個 thread 要做的計算，如圖 B-3.6。
GPU[Device 端(kernel)]
圖 B-3.6 在圖 B-3.1 中自定的廣義編號：三維的廣義編號方式只是多加 2 個變數來編號，例如 多加一個維度可多加一行「j=(BLOCKIDX%Y-1)*BLOCKDIM%Y+THREADIDX%Y」。
以圖 B-3.6 範例比較 CUDA 編號指令定義的編號與我們自定的廣義編號的差別，
如表 B-3.3：
86
表 B-3.3 CUDA 編號與自定編號差異。
CUDA 編號指令定義的編號 自定的廣義編號
BLOCKIDX THREADIDX i=(BLOCKIDX-1)* BLOCKDIM+ THREADIDX
1 1 (1-1)*6+1=1
87
四、計算想法
由「三、CUDA 編號」可知「i」為我們自定的編號，圖 B-3.7 是利用自定 編號命令 4 個 thread 執行計算，其他 thread 閒置，如圖 B-3.8 示意圖。
GPU[Host 端]
圖 B-3.7 在圖 B-3.1 中 Device 端命令 4 個 thread 做運算。
圖 B-3.8 在圖 B-3.6 中自定編號對應到圖 3.4.2.7 的平行運算。
由圖 B-3.7 看出可利用 IF 指令命令哪些 thread 做哪些事，在附錄 C 要討論 如何使用 DO 指令將程式做平行運算。
88
89
GPU 計算所撰寫出的程式碼。
單核 CPU GPU[Host 端]
單核 CPU GPU[Device 端(kernel)]
單核 CPU[執行結果] GPU[執行結果]
圖 C.3 一般 Fortran 修改成 CUDA Fortran 迴圈範例：紅色框框表示相較圖 B-1.1 迴圈新增寫法。
由附錄 B 得知要如何宣告 CUDA Fortran 程式碼、Host 端與 Device 端如何 互傳、如何呼叫 kernel、釋放 device 記憶體方法以及撰寫 Device 端的基本格式 與想法，因此在此不作探討。表 C.1 為除了附錄 B 所修改的部分，對於迴圈 CUDA Fortran 寫法和想法又要修改哪些部分：
90
因此需要宣告一個專門供 thread 計算數值結果存放的陣列，此陣列大小=( thread 總個數)，如圖 C.4。 兩者比較：與一般 Fortran 比較只有在一開始宣告有差異其他都一樣，而可變矩陣釋放 device
記憶體不能使用之前 CUDAFREE 函數的方法，直接使用 DEALLOCATE 即可。
2. 迴圈分成(thread 總個數)等份
圖 C.6 表示將迴圈分成(thread 總個數)等份，迴圈採取跳躍式的疊加，其跳 躍間距等於「BLOCKDIM*GRIDDIM」也就是(thread 總個數)；計算的部份使用 在 Host 端宣告陣列大小=(thread 總個數)的陣列，利用此陣列元素儲存每個 thread 計算後的結果，這樣才不會產生資料碰撞。圖 C.7 為每個 thread 運行的狀況，使
91
用陣列「a_d」將每個 thread 計算結果儲存起來。
GPU[Device 端(kernel)]
圖 C.6 在圖 C.3 中將迴圈分成(thread 總個數)「BLOCKDIM*GRIDDIM」等份，利用每個 thread 做計算。
圖 C.7 每個 thread 運行的狀況：此程式的 grid 內 block 個數「griddx」=3，每個 block 內 thread 個數「tPBdx」=4，因此「BLOCKDIM*GRIDDIM」=12，每個 thread 迴圈跳躍間距為 12，這
樣就不會重複計算到了。
3. 所有 thread 計算數值總加起來
GPU[Host 端]
圖 C.8 在圖 C.3 中每個 thread 總加起來的結果。
在 GPU 內做平行計算，為了不發生資料碰撞，必須將每個 thread 計算出來 的結果儲存起來，最後再傳回 Host 端利用如圖 C.8 方法做總加的步驟。
實際測試結果：
92
Fortran 本身就有一個參數能設定的數值極限，為了防止出錯，因此在此程 式碼迴圈內又加入一個迴圈，如圖 C.9 程式碼所示。
單核 CPU GPU[Device 端(kernel)]
圖 C.9 在迴圈內加入一個迴圈，使它不會因為數值的極限而計算錯誤。
利用圖 C.9 修改過的程式碼，改變迴圈數比較 CPU(編輯器：Intel Fortran) 與 GPU(編輯器：PGI Fortran)的數值結果、時間和倍率(CPU 時間/GPU 時間)，
執行後的結果如圖 C.10 所示。(選取 grid 內 block 個數「griddx」=50，每個 block 內 thread 個數「tPBdx」=512。)
93
107 108 109 1010 1011
107 108 109 1010 1011 1012
value
N
CPU_value GPU_value
圖 C.10 比較單核 CPU 與 GPU 數值、時間、倍率關係：rt為增快倍率，N 為迴圈總個數。
由於結果都使用到 GPU 的優點(核心數多、Device 端內記憶體寬頻快)，因 此增快效應會非常大。假如考慮 device 記憶體很大及 GPU 記憶體不足需分批計 算(詳細 2.2.2 節)、register 不足需使用到 global memory (詳細 2.3.2 節)時，增快
94
效應就不會這麼大了，而庫侖交互作用程式需使用較多的 device 記憶體運算(每 個核心都會消耗記憶體)，必須將程式耗時較久的部分分批計算，由於耗費 device 記憶體多，register 不足要使用到 global memory，因此增快效應沒有此範例結果 大。
95
2. 擁有較複雜的計算，需要較多 register 與使用 global memory 儲存暫存檔，不能將六重迴圈全 部平行運算，因此要考慮讓 GPU 處於最大效益。
(相關內容在 2.3.2 小節與 2.4 小節有提到。)
3. 因為資料龐大，需要創造較多的 grid 來讓 GPU 計算，又為了不讓資料碰撞，每個 thread 的 資訊必須儲存，在做整合時會消耗許多時間，因此必須利用 Tree Reduction 演算法使時間變 快。
(此演算法需使用到 shared memory 與同步概念。)
接下來分成三項探討，分別是傳送方式、平行方式、Tree Reduction 演算法。
一、傳送方式
96 算是最恰當的結果，register 與 global memory 的存取及釋放暫存檔剛好能使 GPU 處於最大效益。而 2.3.2 小節所下的結論是盡量不要在 global memory 存取暫存檔，
但庫侖交互作用程式計算量龐大，必須使用到 global memory 才能做運算，減少 使用 global memory 的方法就是減少需要存取的暫存檔，例如將「A=B+C+D」在 Host 端先計算完「E=C+D」後，再利用 GPU 計算 Device 端的「A=B+E」，這樣 可以減少存取的暫存檔，但壞處就是會破壞原有程式碼的可讀性，因此目前不採 用此法。
三、Tree Reduction 演算法[12]
此部分必須先介紹 shared memory 的宣告與同步的概念。由附錄 B 可知宣告
97
分為 Host 端與 Device 端的宣告，而 shared memory 的宣告只適用於 Device 端的
分為 Host 端與 Device 端的宣告，而 shared memory 的宣告只適用於 Device 端的