科學計算 | Matlab 使用 GPU 並行計算
Matlab下直接使用GPU並行計算（預告）<-- 這預告也貼出來太久了，然而我的大論文還是沒有寫完，但是自己挖的坑一定要填上，我可不是寫小說的。
小引言
說它小是因為它只是博士論文的附錄一部分，但是其實我還是用了很久才學明白的
中心處理器（CentralProcessing Unit, CPU）是計算機系統的計算和控制核心，在軌道設計中使用計算機程序進行研究和仿真都依賴於CPU的強大計算能力，是研究者經常接觸並且熟悉的概念。圖形處理器（Graphics Processing Unit, GPU）是計算機設備中用來進行繪圖運算，支持顯示器設備的芯片。GPU的計算速度和指令復雜度遠不及CPU，但是由於其支持高刷新率高分辨率顯示設備的需求，GPU具有高並行數、大數據吞吐量的特征，並且在這方面的能力遠高於CPU。早在2001年左右，便有人提出了基於GPU的通用計算（General-purpose computing ongraphics processing units, GPGPU）的概念，利用GPU的強大的並行計算能力和相對的低功耗來加速科學計算。基本思路是把需要計算的數據打包成GPU可以處理的圖像信息，然后利用處理圖像信息的運算來實現科學計算。直到2008年左右，NVIDIA公司推出了CUDA（Compute Unified Device Architecture）並行計算架構，代替了GPGPU的概念，並發布了支持C，C++，Fortran的函數庫和編譯器。CUDA是NVIDIA發明的一種並行計算平台和編程模型。它通過利用圖形處理器 (GPU) 的處理能力，可大幅提升計算性能。此后經過一系列的發展和迭代，現在CUDA已經能夠支持符合IEEE標准的雙精度運算，並且支持越來越豐富的流程控制，至於是多維數組和顯存動態分配等功能。
了解GPU的物理結構和動作方式，在一定程序上有助於提高代碼的效率，也有且於開發和調試程序。下面基於本論文使用的GPU型號GTX850M，簡要介紹一下GPU的物理結構。如下圖所示，該GPU核心的實際運算單位有640個，即每個綠色的小塊代表的最基礎的計算單元流處理器（Streaming Process，SP），它們在GPU內部被組織成不同的結構，來實現常用的圖像處理功能。
上圖是計算結構的放大圖，一行中的每4個SP和一些共用的存儲資源等被組織成一個流多處理器（Stream Multiprocessor，SM）；每8個SM，即每8行，添加一些指令調度功能構成一個warp，是GPU中最小的調度單元。在直接調用CUDA進行C/C++編程時，多個計算thread被組成一個block，然后多個block被組織一個grid，最終GPU按照block來調度所有任務在每個warp上進行計算。可以想象，每個block中thread個數，每個grid中block個數的划分，都會影響到能否充分利用GPU的計算能力，而合理的選擇即需要不斷的嘗試，也需要對並行計算算法的理解和較高水下的編輯能力，是一個極其復雜的過程。
好了，看完這些，實際上並不能幫助你了解GPU大概有什么作用，怎么用，還不如直接去看Matlab幫助文檔，上述基本都是廢話，然而我相信等你熟悉了用Matlab調用GPU跑程序以后，你會發現這些廢話基本還算是精髓。
GPU怎么工具
GPU有顯存，還有緩存，也不管它對應上面那個圖中哪些塊塊兒，反正我們也涉及不到這么深入的編程，統稱為顯存好啦。GPU只能直接計算顯存中的數據，就好像CPU只能計算內存中的數據，顯存中的數據或者由CPU從內存移入到顯存，或者直接由GPU在顯存中創建。（很好理解吧，總不能拿個電池在顯存上搓一搓就寫入了數據吧，你知道我在講哪個段子對吧，很冷對吧。）
然后，GPU就可以對顯存里的數據進行各種計算，但是很顯然可以執行哪些計算決定於GPU的能力。比如我的執行 gpuDevice 這個命令后，可以查看顯卡的能力 ComputeCapability 是5.0，我只記得1.3以上就可以支持雙精度運算，即 double，因此 Matlab 也要求必須是 1.3 以上的 GPU 才可以調用。（一般 gtx 8xx 以上的顯卡都是5.0）。
平時不用顯卡做數值計算的時候，它其實就是每秒算60幾張屏幕大小的圖，然后不停地輸送給顯示器。所以可以想象，顯卡的計算能力有多大，然而遺憾的是，喂顯示器數據的計算都很簡單，簡單倒可以把一些算法硬件化來提速，而且單精度就足夠了，而且偶爾錯一兩個像素也無所謂，而且幀數算不過來就跳幾幀也看不出來……總之因為它的出身是顯卡，這些職業特性導致讓它做數值計算會有很多麻煩。當然，Nvidia也在解決，你們肯花錢就好。
好在，用Matlab調用GPU計算，一切都簡單了，反正除了修改代碼和調研，上述困難我都暫時沒有碰到~
在 Matlab 中使用 ode45
這就是本文主要實現的功能，相信有很多人想過用並行計算來加速積分器ode45，但是很遺憾，RKxx這種積分器的數學性質決定它一定是不可以並行的，可以去網上搜，一般中文的答非所問，英文的告訴你不行（有一些其它的積分器有學者在研究並行，基本理念都是用特殊函數基底去近似微分方程的解，然后基函數的系數是有辦法可以並行迭代求解的，不展開，外行）。
此處實現的GPU上並行計算ode45是這樣，每個GPU的小核計算一個ode45，所以同時可以有1000多個ode45在跑，這個意義上的並行。我的應用場景是搜索參數空間里的格點，每個格點都要積分較長時間，並且在積分中判斷每步的狀態。
（寫得好累，語言的表達總是跟不上思路……）
越寫越覺得這個事情不是一句兩句能說清楚的，但是實際真的不難……
啰嗦的詳細步驟
還是一樣的工具箱，在熟悉的parfor下面不遠處，有個GPU Computing的條目，建議把它通讀一遍，內容很少，不比一篇高水平的文章多。
具體要學**的目錄如下
一點都不多，而且只需要學**前三個就好，為了把ode45改寫到可以在GPU上運行，我選中的這第三個才要重點理解。
點進去首先看這部分，學**怎么在GPU上建立矩陣
”Transfer Arrays Between Workspace and GPU“基本就是一句語法，
agpu = gpuArray(a)
這樣就把正常的一個矩陣a變成了GPU上的矩陣agpu（不用非這么命名，也可以叫axianka，嗯，我是在吐槽）。這是之前講的CPU把數據從內存移動到顯存，可以想象，GPU自動生成數據的功能就是上面那個”Create GPU Arrays Directly“，包含一些隨機數函數啥的，而且Matlab一再強調與CPU下的不同，用的時候再看就好了，反正我是用不到。
然后就是是學**怎么用這些gpuArray類型的數據。哦對，上述這些顯存上的數據，類型是gpuArray，如下圖whos a agpu的結果
想要讓CPU用這些gpuArray的時候，用gather可以把它們取回到內存變成double（此處不要問問題，不要多想，我假裝啥也不知道，用到時候查查文檔或者google一下就可能解決的，相信我~）
繼續
然后就是是學**怎么用這些gpuArray類型的數據，學**下圖中的Run Built-in Functions on a GPU部分
這部分介紹了Matlab內嵌的，可以支持gpuArray數據類型的數據的函數。就是說，你給這些函數輸入之前的agpu，它們是會傻傻地算的，並且計算是發生在GPU上，結果也是在顯存中。但是，……算了一會兒下一小段再但是。
再來學**如何在Matlab上運行自定義函數，這里也就是我們將要修改的ode45函數
可以看到，我在這里強調了element-wise函數，在Matlab上運行的自定義函數，必須是每個輸入都是1維的，element-wise，就像 .*，./，.^ 什么的一樣的element-wise operator。這個element-wise的要求，是對自定義函數內部全部操作而言的，即沒有向量操作，沒有矩陣操作，沒有矩陣乘，沒有求逆，沒有各種……甚至不能有索引出現，當然，沒有前面那些還要索引做什么。我的GPU版本的ode45中矩陣都是折成分量，把乘法一點點寫出來的，就跟當年學C時候一樣，感覺自己像個小孩子。
哦，值得說一下，在內部不能動態定義數組，比如a = [a,1]這種操作是不可以的，因為CUDA好像不支持這功能，而且，這不就出現數組了嘛，記住原則是element-wise。
然后最令人崩潰的一個表格來了，下圖，支持的Matlab code
我當時看到這個的時候， 直接略過，這什么嘛，又重復一遍，和前面那個圖不是一樣樣的嘛。
當我在GPU什么都跑不起來的時候，自定義函數不停報錯的時候，我又看到這部分的時候，唉，讓我靜靜……再次告訴我，別以為Matlab里有一句是廢話，一定仔細看，不然全是坑……
說正經的，這些函數Matlab稱之為Supported Matlab Code的意思（我揣測）是說這些函數是所有你可以在GPU上運行的自定義代碼里使用的函數，即你寫GPU自定義code只能用這些，別的不支持。比如**慣了寫成sind(30)的必須現在寫成sin(30/180*pi)，明白了？好在pi是支持。仔細看看，其實跟C的基礎數學庫差不多，多個inf，pi，randi啥的，值得注意的是支持全部的常用流程控制for，while，continue，然而switch，try，catch這些高級貨是沒有的。
好了，有這些知識 ，基本上就可以改你手頭的代碼了，改成符合上述要求的：
-
只用了這些基本函數和流程控制
-
通篇element-wise操作
-
其它……
（寫其它是為了嚴謹，畢竟我也是新手門外漢二把刀，純粹是為了提高公眾號知名度才來寫這**的，才不是因為真心熱愛學術喜歡搞技術呢~~話都說這兒了，得幫我推廣一下吧，您呢~）
改完以后，用arrayfun這個函數調用，像這樣：
arrayfun的第一個參數是在GPU上運行的自定義函數fungpu，這里它內部調用了下面要講的ode45修改版，后續的參數依次是fungpu的參數，其中只要有一個是gpuArray類型的（這里的t10vector），其它參數全部會被傳送到顯存，然后fungpu就會在GPU上運行。
下面我就開始介紹修改ode45的過程，
（其實想想你跟我修改ode45並沒有多大用，重點還是上面講的那些多看看多琢磨，平時科研碰到用得上的問題留個心，有功夫了就嘗試一下，有問題多思考或者去google或者來交流討論，這樣才有效果）
考慮到這么深切的原因而不是因為我懶，我就精簡地講一下。
先看下ode45的源碼，我知道其實很多人都已經讀過了。
沒讀過也不要緊，相信我，如果不是有需求被逼，你情願自己寫一個rk45或者rk78也不願意讀這功能強大的代碼（不得不說官方的代碼寫得很好，好用又好讀，還好改）
上面紅杠的部分，基本是要刪掉的，因為在GPU上不支持，這些都是一些高級功能，比如Event啊啥的，還有一些通用性的功能。
我們需要的功能主要在紅框中，主要是ode45的rk45算法部分。
（注意上面的代碼我折疊了很多，折疊了很多，折疊了很多，而且不是嚴格按照我划的改，隨手一划，你懂得，隨手）
修改這部分代碼
比如我是積分三維空間中的軌道，所以輸入的變量是6維，時間是1維，還有些其它質量比例什么的參數。
改后的結果如下：可以看到此時精度也作為獨標量直接輸入，而不能再用Name-Value pair的形式輸入，其它的功能都被精簡掉了，因為GPU不支持。
可以看到原先的矩陣操作都被寫成了分量形式，其實也不難，對吧。
下面看到的一些中文注釋，很多是在讀ode45代碼時根據當時理解添加的。注釋總不嫌多，對吧。
看，RK45每步求微分的過程，也要用分量形式寫出來了，其中的OdeRtbpElementWise函數應該能猜到是什么吧，就是要積分的動力學方程，也被寫成了分量形式的輸入輸出，1個時間，6個狀態。
誤差分析也要一點一點拆，就當復**了一遍數值分析，對吧。你看這大段大段的英文原版注釋，我都沒有動過的，說明修改的工作量其實也不大，對吧。
這就結束了，有沒有覺得很神奇？這貨就可以在GPU上積分啦！
等等，並沒有，還沒有講輸出。因為GPU不支持自定義函數的索引功能，不支持動態數組，我們不能像正常ode45一樣輸出軌道的時刻序列，只能輸出某一個狀態，這里我們選擇末狀態[t,y1,y2,y3,y4,y5,y6]。
這樣便存在一個問題，如何驗證整個過程中軌道是正確的？很簡單。
用ode45從t0到t1積分一遍，取ii=1:1000個點畫圖。
用ode45gpu積分一遍，分別從t0到t0+ii/1000積分ii=1:1000條軌道，注意這1000條軌道是同時被計算的，所以和從t0到t1直接積分的耗時只差一點點一點點。
然后gather回來重疊畫圖比較
這個圖就是這么畫的，所以你看這么多小點點~~怎么樣，不難吧？
還有一個調試的問題沒有講到，是程序就免不了調試。其實也不用講，都在Matlab下了，還要怎么講，直接在CPU下串行調試就好啦，保證element-wise的原則，正常調試嘍。
應該就講完了，再講就是各種細節了，我覺得沒必要。
最后上個圖說明一下有效性：
上圖是在星歷模型下搜索得到的某種結果圖，下圖是在某種近似模型下搜索得到的同種結果圖，關鍵點在於分布規律是一樣的！說明GPU的積分精度和直接在CPU上的RK78積分精度是可比擬的，是可用的。
上圖大約14小時，分12x16線程計算，2.3GHz的CPU，用的UPC的服務器，還是FORTRAN77的計算效率；
下圖大約13小時，顯卡GTX850M超頻后1.2GHz，Matlab自動划分各種Warp，Block啥的。
可見，加速能力還是相當強悍的，才是個不到1000塊錢的筆記本顯卡。
當然這里有一個問題，GPU如何使用星歷模型，竊以為這個問題並不存在一個簡單的可以僅依靠Matlab解決的途徑，所以，管它呢……
寫得着實太粗糙，如果有問題或者有錯誤，還請右下角隨意“留言”，我會認真回復的。
如果有想以我這個為基礎研究一下的，我可以考慮分享我的代碼給你，暫時還沒有想好形式，而且應該要在6月份以后了。
-
歡迎留言，歡迎討論。
-
歡迎推薦給其它人。
-
點擊“閱讀原文”查看更多。