MATLAB 到高性能C 的代碼轉換系統設計與實現*

余 澤 霖，徐 云

(1.中國科學技術大學 計算機科學與技術學院，安徽 合肥230026；2.安徽省高性能計算重點實驗室，安徽 合肥230026)

0 引言

MATLAB 由于其強大的功能、友好的開發界面和簡單易用的編程語言形式而被開發人員廣泛使用，許多企業和科研單位都存在著大量的MATLAB遺留代碼(Legacy Code)。但由于MATLAB 程序運行依賴于MATLAB 軟件系統，部分工作平臺無法滿足要求，并且它的運行速度相較于C 等更底層的語言慢，無法滿足高性能計算的需求，通常需要將這些遺留MATLAB 代碼轉換為C 代碼。而人工轉換代碼需要付出高昂的人力成本，因此，理想的解決方案是將MATLAB 代碼自動轉換成C 代碼。

程序轉換可以節省軟件開發的成本，并在軟件的移植、重用、更新和編譯等方面有著重要的應用[1]，因此在國內外有大量的對編程語言轉換的研究。本文對已有的編程語言轉換方法進行了借鑒和改進，提出了一個MATLAB 到C 的自動轉換方法，著重解決了以下三個關鍵問題：一是引入了一種基于抽象語法樹(Abstract Syntax Tree，AST)的中間表示(Intermediate Representation，IR)來解決兩語言間的差異問題，為MATLAB 和C 轉換過程中提供統一的接口，避免為兩種語言的語法結構映射制定繁瑣的規則。二是由于C 中變量需要聲明后使用，而MATLAB 不需要，因此在轉換時需要對MATLAB 中的矩陣變量的類型信息(矩陣的大小和矩陣元素的類型)進行識別。已有許多對此類高級語言的變量類型信息識別的研究成果[2-6]，本文借鑒并實現了MATLAB 到C 的轉換過程中的類型信息推導。三是MATLAB 中大部分表達式是矢量間的運算，直接轉換到C 語句會導致性能下降，并且很多內置的科學計算函數沒有直接對應的C 代碼。本文設計了轉換算法，使用C 數學核心庫(Math Kernel Library，MKL)中的函數為MATLAB 中的矢量運算和庫函數生成高性能C 代碼段，解決轉換困難和性能問題。本文實現的系統能幫助開發人員將MATLAB 代碼自動轉換生成C代碼，提升代碼性能，減少人工成本。

1 相關研究

1.1 編程語言轉換研究現狀

目前，編程語言轉換方法大致可以分為兩類。

一類是將原語言的語法結構一一轉換到目標語言的語法結構。李友仁[1]等人提出了一種語義解釋的轉換方法，為Prolog 的不同組成部分制定相應的轉換模板，并應用該方法開發了Prolog 到C 的自動轉換系統；Boyle[7]等人通過將Lisp 中的多種語言特性逐步轉換為對應的Fortran 語句實現Lisp 到Fortran的轉換；Moynihan[8]等人通過定義固定的轉換規則把Pascal-SC 中的語句轉換到對應的Ada 語句，輔助將Pascal-SC 代碼轉換到Ada 代碼；Demaine[9]提出了將C 中的指針等價地轉換成Java 中的引用的方法，并基于此實現C 到Java 代碼的自動轉換。這類方法需要為原語言的每種語言結構制定相應的轉換方案，對語言依賴性高，當出現新的需要轉換的結構時往往需要對原方案重新設計，缺乏通用性。

另一類是通過引入一層IR 來橋接語言差異。De Rose[10]和Reis[11]等人通過基于SSA(Static Single Assignment)的IR，分別實現了MATLAB 到Fortran 90 和C 代碼的轉換和優化；石學林[12]等人引入了兩層IR來實現Cobol 到Java 代碼的自動轉換器；Zhao[13]、Paulsen[14]等人分別通過基于AST 的IR 實現了PHP和MATLAB 到C++的自動轉換。這類方法使用IR作為原語言和目標語言之間的接口，存儲語言中的數據結構和語法結構信息，從而避免依賴特定語言，當某一語言的轉換需求發生變更時只需要改動其與中間層的轉換關系即可，避免對整個轉換框架做較大改動。而且在IR 中可以存儲代碼信息，方便后續的分析和優化。但是多引入一層IR 會帶來額外工作量，需要選擇合適的IR 以及生成方法。

在MATLAB 代碼轉換中針對如何將高性能運算轉換到C 代碼的問題也有相應的研究。Banerjee[15]等人的方法是使用目標語言開發相應的接口，該方法工作量較大，難以應對需求變化，且轉換后代碼難以保證高性能；Paulsen[14]使用Armadillo 庫中的函數來替換相對應的MATLAB 函數，能有效避免轉換困難，但是存在轉換后代碼性能問題，同時，對非函數形式的表達式的轉換也有所欠缺。對此，需要設計相應的算法解決。

另外，商用的工具如MATLAB Coder[16]還沒有公開的實現方法。

1.2 矩陣變量類型識別研究現狀

目前已有許多針對高級語言的矩陣型變量的類型識別研究。Schwartz[2]通過對高級語言SETL 做數據流分析，獲取程序上下文中向量的長度信息并生成相應的優化代碼；Chatterjee[3]在VCODE 編譯器中實現了對向量的長度推斷；Joisha[4]、Chevalier-Boisvert[5]和Bispo[6]等分別在MATLAB 的優化或轉換工作中對矩陣變量的類型信息進行了推斷。這些方法的推斷過程基本上是通過代碼中入口函數的參數對代碼語句進行數據流分析，逐條推導表達式結果變量的類型信息。推導方法可以應用于本文的研究中。

2 方法設計與實現

2.1 基于AST 的中間表示設計

本文的研究中引入了一層中間表示(Intermediate Representation，IR)來避免對語言特征的依賴，該IR以AST 為基礎構建。AST 是源代碼語法結構的一種抽象表示，以樹狀的方式呈現，其中每個節點表示源代碼中的一個結構，一個AST 的示例如圖1 所示。AST 被廣泛應用于代碼編譯等領域。

圖1 表達式y=1+2*x 的AST

在構建C 代碼時，需要對變量先聲明后使用，所以需要知道相應MATLAB 變量的類型信息，本文將這些信息存儲到AST 中代表矩陣變量的節點中，從而構成以AST 為基礎的IR。從源代碼生成AST需要對代碼進行語法分析，大部分研究中都是手動實現這一過程，是一項繁瑣的工作。本文使用ANTLR(ANother Tool for Language Recognition)來實現。ANTLR是一款強大的語法分析器生成工具，可用于讀取、處理、執行和翻譯結構化的文本或二進制文件，被廣泛應用于學術領域和工業生產實踐，是眾多語言、工具和框架的基石[17]。使用ANTLR 生成MATLAB的語法分析器后，就能進一步生成AST。

與編譯器技術中IR 會深入到底層表示不同，本文的IR 層只是概念模型，MATLAB 代碼轉換到IR 后會直接轉換到C 代碼，不涉及更底層的語言，減少了轉換的層次。

2.2 功能模塊劃分與工作流程

由于引入了IR 層，需要分別將MATLAB 代碼轉換到IR 以及將IR 轉換到C 代碼，將這兩個階段的工作劃分為前端轉換模塊和后端轉換模塊，通過中間層模塊連接起來。

在前端轉換模塊中，需要從MATLAB 代碼生成AST，并識別代碼中矩陣變量類型信息。上述工作分別劃分給AST 生成模塊、矩陣大小推導模塊和矩陣元素類型推導模塊。

后端轉換模塊中，需要對IR 中的AST 進行遍歷，轉換為相應的C 代碼。該過程中的矢量運算和函數調用需要使用MKL 庫中相同功能的C 函數進行替換生成相應代碼段。上述工作分別劃分給解析模塊和替換模塊。

各模塊及其子模塊的關系如圖2 所示。

圖2 系統功能模塊劃分

圖3 系統處理流程

系統將MATLAB 代碼轉換為C 代碼的工作流程如圖3 所示。首先輸入待轉換的MATLAB 代碼文件，進入到前端轉換階段，該階段第一步是將MATLAB代碼經過詞法和語法解析后生成AST，然后添加MATLAB 代碼的輸入參數矩陣信息，通過AST 結合輸入信息推導出其他矩陣變量的大小和元素類型信息，加入到AST 中得到IR。然后進入到后端轉換階段，對IR 進行解析，將每個AST 節點轉換為相應的C 代碼結構，并輸出為相應的C 代碼。同時，對代表矢量運算和函數調用的節點，將其交由替換模塊處理，輸出相應的MKL 庫函數。最終得到高性能的C 代碼文件。

2.3 功能模塊實現

2.3.1 前端轉換模塊實現

下面說明前端轉換模塊中三個子模塊的設計與實現。

(1)AST 生成模塊實現。本文使用ANTLR 來實現MATLAB 的詞法和語法分析器。使用ANTLR 時需要提供待處理語言的語法文件，通過閱讀MATLAB文檔，編寫出MATLAB 的語法文件，圖4 是該語法文件的部分示例。

圖4 MATLAB 語法文件部分示例

(2)矩陣大小推導模塊實現。矩陣大小的推導方法是根據矩陣運算以及MATLAB 函數的規則來制定的。首先是矩陣運算，當知道參與運算的矩陣大小以及運算符，就可以推導出結果矩陣的大小，比如兩矩陣A 和B 相乘，結果矩陣的行數為A 的行數，列數為B 的列數。這些矩陣運算的規則有限，可以實現到代碼流程中。

對于函數調用，新產生的矩陣的大小與參數矩陣的大小和函數有關，可以通過這樣一個三元組來定義：＜函數名，參數矩陣大小，輸出矩陣大小＞，即函數名和參數矩陣大小可以確定輸出矩陣的大小。為便于管理和擴展，將這些規則按照三元組形式寫入到文件中構成規則文件，推導模塊在遇到函數調用時通過規則文件查詢相應的規則。當待處理的MATLAB 代碼中出現新的函數調用時，在規則文件中加入新的規則即可。

有了上述推導方法，在給出MATLAB 代碼輸入參數矩陣的大小之后，對2.1 節中生成的AST 進行遍歷，對其中遇到的計算節點和函數調用節點推導結果矩陣的大小，遍歷完成后就能得到所有矩陣變量的大小并存儲到相應的矩陣變量節點中。

MATLAB 代碼入口輸入參數矩陣的大小信息通過一個三元組來定義：＜矩陣名，行數，列數＞，將所有輸入參數矩陣的大小按照三元組形式寫入到文件中，提供給推導模塊使用。

(3)矩陣元素類型推導模塊實現。在MATLAB中，矩陣的元素可以有多種數據類型，具體如表1所示。

表1 MATLAB 數據類型

MATLAB 規定了矩陣運算時，結果矩陣的元素數據類型與運算數的元素數據類型的關系，規則如下：①如果運算數中含有整數矩陣，其他運算數也必須是同一類型的整數矩陣，或者是雙精度的浮點數標量，而結果仍然是原整數型；②浮點數運算時，如果運算數中有單精度矩陣，那么結果矩陣也是單精度的。這兩個有限的規則實現到代碼流程中。對于函數調用的輸出矩陣元素類型規則，同前文一樣，定義為一個三元組：＜函數名，參數矩陣元素類型，輸出矩陣元素類型＞，通過文件的方式提供給推導模塊。

矩陣元素類型推導過程與前文中矩陣大小推導過程相似，也是對AST 進行遍歷，同時進行推導。該過程需要輸入參數矩陣的元素類型信息，通過一個二元組來定義：＜矩陣名，元素類型＞。將所有輸入參數矩陣的元素類型按照二元組形式寫入到文件中，提供給推導模塊使用。

2.3.2 后端轉換模塊實現

下面說明后端轉換模塊中兩個子模塊的設計與實現。

(1)解析模塊實現。解析模塊的工作是對IR 中的AST 進行遍歷并生成相應的C 代碼。MATLAB 語句生成對應的C 語句時，由于變量需要聲明，因此先定義MATLAB 變量類型在C 中的映射規則。

MATLAB 中的矩陣使用C 中的數組進行模擬。一個大小為m×n 的矩陣在C 中定義為一個長度為m×n 的數組。一行一列的矩陣是特例，在C 中當做一個標量來定義。在本文使用的硬件平臺是32 位的情況下，MATLAB 的數據類型對應的C 的數據類型總結于表2 中。

表2 MATLAB 數據類型對應的C 數據類型

當解析過程中發現新的變量時，它的大小和數據類型已經在前端轉換過程中被推導出來并存儲在IR 中的AST 節點中，將其聲明補充上即可。

在遍歷AST 的時候，當前AST 節點決定應該進行的相應操作。MATLAB 整體代碼生成的AST 的根節點應該是語句塊，解析從該語句塊開始。解析過程的算法如下：

算法1：AST 解析算法

輸入：IR 中AST 根節點Root；

輸出：C 代碼。

其中，ASSIGNMENT 代表賦值語句類型，VEC_OP和FUNC_CALL 代表MATLAB 的矢量運算和庫函數調用，IF 代表if 語句類型，LOOP 代表循環語句類型(for/while 循環)。

(2)替換模塊實現。替換模塊的工作是將MATLAB中的矢量運算和庫函數調用替換為MKL 庫函數形成的C 代碼段。首先是庫函數調用，可以在MKL 庫中尋找相應功能的函數來進行替換，這部分的替換較為簡單，替換規則與下文矢量運算的替換規則一起考慮。

對于矢量運算，首先考慮表達式不存在嵌套(即一個表達式只對應一個MKL 函數)的情況。替換的難點主要在于三方面：①確定矢量運算所對應的MKL 函數；②同一個矢量運算表達式由于矩陣類型不同會對應不同的MKL 函數(一對多)；③找到MKL函數后變量如何填入。

本文設計了一個匹配算法解決上述問題。首先定義MATLAB 語句替換為MKL 函數的替換規則，該規則通過一個三元組來定義：＜MATLAB 語句模板，參數類型，MKL 庫函數模板＞，并通過文件來管理和擴展。對上述問題解決方法為：①MATLAB語句模板和參數類型共同匹配一個對應的MKL 庫函數，相應的匹配算法解決了函數的確定問題；②參數類型也參與了函數的確定，為不同的參數類型添加新的替換規則，來解決一對多的問題；③規則中使用的是模板語句，參數使用占位符，當輸出具體的C 語句時才根據實際代碼中的變量名來替換。

匹配算法根據實際語句的AST(下稱“SAST”)與MATLAB 語句模板的AST(下稱“TAST”)進 行。替換模塊初始時會讀取替換規則文件并將MATLAB 語句模板轉換為TAST。判斷AST 是否相同即是判斷樹結構是否相同，當前的節點類型相同并且子樹也相同時說明AST 相同，可以通過對AST 做遍歷完成。在該過程中，變量獲取以及類型信息判斷較為簡單，在算法說明中略去。

如果考慮表達式存在嵌套，則需要在AST 匹配失敗時對嵌套的表達式先輸出對應MKL 函數，然后相應修改AST，之后繼續匹配，直到輸出最終結果。匹配算法如下：

算法2：匹配算法

輸入：SAST 和替換規則文件；

輸出：MKL 庫函數調用語句。

該算法是遞歸的，這是由于表達式的嵌套（如加法和乘法混合等）。嵌套部分將被提取出來輸出對應的MKL 庫函數調用語句，然后修改AST，再重新進行匹配，直到完成匹配過程，算法在第4 行結束。嵌套的情況舉例如圖5 所示，表達式y=a+b-c中(a、b、c 均為長為n 的向量)，a+b 的結果先存放到臨時變量t 中并輸出相應的MKL 函數語句，然后AST 中用新節點t 替換子樹，重新匹配y=t-c 得到最終結果。

圖5 嵌套運算匹配示例

3 實驗結果與分析

3.1 實驗配置

本文測試轉換系統對MATLAB 代碼轉換生成C代碼的正確性以及性能。由于MATLAB 在科學計算和數據挖掘等方面有著廣泛應用，本文選取迭代法解線性方程組和K-means 聚類算法的MATLAB 程序來測試系統。測試1(迭代法求方程組解)分別測試解100、200、400 階線性方程組的時間，測試2(K-means 聚類算法)分別測試將100、200、400 個四維點劃分為4 個聚類的時間。C 程序的數據精度均為double。

生成的C 程序運行的硬件環境為：Intel Core i7-8550U，1.80 GHz，4 核8 處理器CPU；8.0 GB 內存。

3.2 與人工編寫代碼對比

為了測試本文系統轉換生成的C 代碼的正確性以及與人工編寫代碼的性能對比，使用專家級水平的人工代碼進行對比測試。本文系統轉換生成的C 代碼與人工編寫的C 代碼均使用MKL 函數庫。對兩個測試用例，運行系統轉換生成的C 代碼和人工編寫的C 代碼，使用相同的輸入，兩者的運行結果一致，說明了本文系統能生成正確的C 代碼。

兩個C 程序在測試用例上的運行耗時如表3和表4 所示。

表3 與人工編寫代碼對比測試1 結果

表4 與人工編寫代碼對比測試2 結果

由結果可以看出，系統轉換生成的代碼運行時間和專家級水平的人工代碼運行時間相差不多，排除運行環境浮動干擾帶來的誤差，兩者的性能相當。本文的系統可以幫助開發人員將MATLAB 代碼轉換到性能更好的C 代碼上，并提供相當于專家級水平的人編寫的代碼性能，降低開發成本。

3.3 與已有方法對比

選擇兩個已有方法將MATLAB 轉換生成C++代碼，這兩個方法生成的C++代碼實際上并未使用較多C++特有的頭文件，因此將其與本文的C 代碼進行比較。第一個是MATLAB Coder，使用時需要指定矩陣變量的元素類型和大小，針對不同的測試需要生成多份C 代碼。第二個是文獻[14]中的matlab2cpp，使用過程中需要根據它的自定義規則將部分不能自動轉換的代碼段替換為正確的C 代碼段。生成的代碼運行耗時測試結果如表5 和表6 所示，其中本文方法的結果沿用3.2 節中的數據。

表5 與已有方法對比測試1 結果

表6 與已有方法對比測試2 結果

結果表明，MATLAB Coder 轉換生成的C 代碼的性能最差，以它的運行時間為基準，matlab2cpp 和本文方法對它的加速比如表7 所示。

表7 測試結果相比于MATLAB Coder 的加速比

結果顯示本文方法生成的C 代碼的性能略優于matlab2cpp，兩者的效果均優于MATLAB Coder。對MATLAB Coder 生成的C 代碼分析可以得知，MATLAB Coder 基本上是將MATLAB 代碼中的矢量運算以及庫函數展開為循環，這樣C 代碼的性能對比使用了性能庫的C 代碼性能顯著較差。而matlab2cpp使用了armadillo 庫來提供大部分MATLAB 矢量運算和庫函數的轉換，該庫的語法與MATLAB 相似，在轉換上能提供部分便利性，且該庫底層也集成了高性能的線性代數庫，但是由于該庫對底層的封裝，在性能上會帶來部分損失。同時，matlab2cpp 也不能保證轉換的完全自動化，因為armadillo 庫對于MATLAB中部分矢量運算和庫函數有所缺失。而本文通過對矢量運算和庫函數設計的替換算法能夠保證轉換的自動化以及擴展性，同時轉換后的C 代碼性能在MKL 的提升下也能優于已有的方法。

4 結論

本文設計并實現了一個MATLAB 代碼自動轉換到高性能C 代碼的系統。針對大部分自動轉換方法對特定語言依賴較大的問題，本文引入IR 層，并將系統的整體框架分為前端轉換模塊、中間層模塊和后端轉換模塊三大模塊進行實現，對傳統方法中存在的一些問題作出了改進。實驗結果表明，該系統自動生成的高性能C 代碼性能與專家級水平的人工實現的相當，且優于已有的方法。

目前，該系統只支持MATLAB 代碼的矩陣運算代碼轉換，對于含更高維度的矢量代碼轉換還未實現，并且轉換后的C 代碼僅支持使用MKL 庫提升性能，局限性較大，這是本文進一步的研究方向。