面向FT-M6678 的對稱矩陣特征值求解算法實現(xiàn)與優(yōu)化

于立，韓林，羅有才，商建東

（1.鄭州大學計算機與人工智能學院，河南 鄭州 450001；2.國家超級計算鄭州中心，河南 鄭州 450001）

0 引言

隨著計算機技術的不斷發(fā)展，其在工程領域方面的應用越來越多，其中離不開數(shù)字信號處理器（DSP）的支持［1］。DSP 由于其超低功耗，傳統(tǒng)意義上一直是許多嵌入式系統(tǒng)的基石［2-4］。在國產(chǎn)DSP 芯片中，由國防科技大學研制的飛騰系列DSP 芯片采用了更為先進的制造工藝，進一步提高了芯片主頻，其通用性能較強，并且在專用性能方面采用流處理器方式，在提高專用計算性能的同時，也能顯著降低功耗［5］。其中，F(xiàn)T-M6678 不但具有強大的浮點運算性能支持，而且具有易擴展、軟件兼容度高的特點，因此，F(xiàn)T-M6678 適合作為研究算法和優(yōu)化的平臺，同時，構(gòu)建國產(chǎn)芯片的生態(tài)體系，對提升國防安全、民生經(jīng)濟、國家綜合科技實力都有重大意義。

對稱矩陣特征值求解（SYEV）是線性代數(shù)數(shù)學庫中的一個重要組成部分，其主要負責求解給定對稱矩陣的特征值和特征向量［6-7］，解決該問題在許多科研工程中都有至關重要的作用，如量子力學、航空動力學、計算流體力學等應用領域［8-9］?？v觀近些年來，國內(nèi)學者對線性代數(shù)運算庫愈加重視，嘗試在不同國產(chǎn)平臺適配并改進其所需的線性代數(shù)相關應用的算法：文獻［10］基于飛騰2000+平臺，使用OpenMP并行技術對BLAS 庫3 級函數(shù)進行優(yōu)化；文獻［11］在神威眾核平臺SW26010 上實現(xiàn)了數(shù)學庫xMath2.0，并對其進行針對性優(yōu)化，該庫包含LAPACK 等多種數(shù)學庫，完善了神威平臺的數(shù)學庫運算需求；文獻［12］面向鯤鵬處理器平臺，通過向量化與多核并行實現(xiàn)了LAPACK 庫中對稱矩陣方程求解例程的性能優(yōu)化，完善了鯤鵬平臺線性代數(shù)運算基礎；文獻［13］基于龍芯3A 平臺對LAPACK 線性方程組求解函數(shù)并行化，并且在龍芯平臺上使得該算法在原LAPACK 基礎架構(gòu)上加速2 倍。目前未見有基于FT-M6678 平臺實現(xiàn)與優(yōu)化對稱矩陣特征值求解的相關研究，且FT-M6678 現(xiàn)有的數(shù)學庫并不能很好地滿足線性代數(shù)相關的問題需求，因此，將本文研究內(nèi)容在FT-M6678 上實現(xiàn)對滿足飛騰平臺線性代數(shù)計算需求具有重大意義。

本文研究內(nèi)容主要包含如下方面：首先介紹FT-M6678 的體系結(jié)構(gòu)，分析SYEV 的算法原理，在FT-M6678 平臺上實現(xiàn)該算法；然后對整個算法程序的運行熱點進行測試分析，針對熱點情況對SYEV進行優(yōu)化，采用基于FT-M6678 的編譯優(yōu)化、訪存優(yōu)化以及向量并行化方法，重點優(yōu)化耗時高的部分；最后驗證正確性并分析優(yōu)化性能，正確性驗證使用線性代數(shù)庫LAPACK 官方提供的驗證集，性能分析使用不同規(guī)模的輸入矩陣，矩陣由LAPACK 官方提供的生成矩陣函數(shù)生成，分析各規(guī)模數(shù)據(jù)下的優(yōu)化性能，對比優(yōu)化前后的性能差異，并與TMS320C6678平臺上的性能進行對比。

1 FT-M6678 體系結(jié)構(gòu)與算法原理

1.1 FT-M6678 體系結(jié)構(gòu)

FT-M6678 是由國防科技大學自主研發(fā)的數(shù)字信號處理器，單核理論計算峰值為16GFlops，具有高性能、易擴展的特點，可用于機載信號處理［14］、圖像處理［15］、電子對抗等多個領域，滿足嵌入式高性能計算需求，且兼容德州儀器（TI）的DSP 芯片TMS320C 6678 處理器 的軟件 開發(fā)環(huán) 境［16-17］。FT-M6678 采 用新型Key Stone 多核架構(gòu)和增強型內(nèi)核M66x，為8 核DSP，每個核的處理頻率為1 GHz。FT-M6678 的架構(gòu)如圖1 所示，其從功能上主要分為CorePac 內(nèi)核、核外存儲系統(tǒng)、互連網(wǎng)絡、高速接口、低速接口、集成外設、全局控制寄存器、自舉復位［18-19］等8 個部分。每個核具有1 級和2 級緩存，核外存儲有共享內(nèi)存通信（SMC）存儲以及內(nèi)存DDR3 存儲，每個SMC 共享存儲空間為4 MB，每個DDR3 最高可用2 GB。

由于FT-M6678 的硬件特性，其強大的浮點運算能力適用于線性代數(shù)中各問題求解的矩陣運算，這能很好地發(fā)揮飛騰平臺的優(yōu)勢。

1.2 SYEV 算法原理

SYEV 算法求解特征值與特征向量使用QL、LQ或QR 分解算法。QR 算法的具體形式如式（1）所示，其中，A為對稱矩陣，Q為正交矩陣，R為上三角矩陣。QL 或LQ 算法與此原理相同，不同之處為L代表使用的是下三角矩陣，且LQ 算法為右乘操作。

QR 分解算法是一個不斷迭代的過程［20］，迭代過程如式（2）所示，其中，上標T 表示該矩陣的轉(zhuǎn)置（下文相同），變量k為迭代次數(shù)。

使用正交矩陣Q不斷左乘A，將A的某些元素消零，其中正交矩陣Q的生成方法是通過HouseHolder反射變換。HouseHolder 反射變換是通過初等反射矩陣的連乘將非奇異矩陣轉(zhuǎn)化為上三角矩陣［21］，其只需求解n-1 個反射矩陣及其乘積，復雜度低于Givens 分解［22-23］。HouseHolder 反射變換原理如下：

假設非0 向量α∈Rn，I為單位矩陣，則滿足：

如 式（3）所 示，形式如n階矩陣P即 是HouseHolder 矩陣，該矩陣滿足對稱且正交的性質(zhì)。將式（3）修改為式（4）所示的形式，H即為算法每一次反射的反射因子，其中，I為單位矩陣，v為單位正交矩陣，t為變換輸入的標量參數(shù)。

因此，迭代過程可轉(zhuǎn)化為式（5）到式（6）所列出的公式，其中，H（x）為第x代的HouseHolder 反射因子，R為上三角矩陣，每一次迭代將一列的對角線以下元素清零，使得矩陣無限趨近于一個上三角矩陣。

經(jīng)過上述不斷迭代，最終可收斂簡化成式（7）所示的形式：

其中：Λ為對角矩陣，其對角線元素即為所求對稱矩陣的特征值，而XT中包含與之所對應的特征向量。

2 面向FT-M6678 的SYEV 實現(xiàn)與優(yōu)化

2.1 面向FT-M6678 的SYEV 算法實現(xiàn)

SYEV 算法實現(xiàn)過程如圖2 所示。算法先將對稱矩陣轉(zhuǎn)換為三對角形式，轉(zhuǎn)換方式為正交相似變換，這一步是為了簡化數(shù)據(jù)量，減少后續(xù)計算的運算開銷。

圖2 SYEV 算法實現(xiàn)過程Fig.2 Implementation process of SYEV algorithm

程序在只求解特征值時，直接使用QL 或QR 算法計算一個對稱矩陣的所有特征值；在既求解特征值又求解特征向量時，SYEV 會先生成正交矩陣，生成的方式是通過采用HouseHolder 反射變換，與將對稱矩陣轉(zhuǎn)換為三對角步驟中使用的正交矩陣生成方式相同，然后計算所有特征值和對應的特征向量，使用的算法是隱式QL 或QR 算法。

本文選用單精度浮點型數(shù)據(jù)類型進行算法實現(xiàn)，算法部分計算需要借助底層一些基本矩陣操作（如對稱矩陣與向量乘等操作）來完成，該操作部分使用BLAS 二級函數(shù)［24-25］輔助完成。

2.2 熱點分析

分別生成100×100、500×500、1 000×1 000 規(guī)模的對稱矩陣，測試SYEV 的性能并且分析函數(shù)熱點，單個函數(shù)運行耗時占比如圖3 所示。

圖3 運行熱點分析圖Fig.3 Analysis diagram of operational hotspots

圖3 僅顯示出耗時排在前9 位的操作運算，再后面的操作運算耗時幾乎可以忽略不計。由該圖可知，SYEV 的運行熱點分布比較規(guī)律，注意到子程序矩陣平面旋轉(zhuǎn)操作是耗時最多的，單個函數(shù)耗時占比超過總程序運行時間的50%，并且隨著計算規(guī)模的增大，其占比也會變大，而剩余耗時的部分是一些BLAS 二級函數(shù)，如計算矩陣向量乘以及秩更新等，但都與矩陣平面旋轉(zhuǎn)操作無法相比，因此，本文重點優(yōu)化該操作部分。矩陣平面旋轉(zhuǎn)功能上負責實現(xiàn)將一個矩陣從左側(cè)或右側(cè)進行平面旋轉(zhuǎn)，在計算HouseHolder 反射因子部分會大量運行該操作。

2.3 編譯優(yōu)化

對于編譯階段的優(yōu)化主要是根據(jù)不同的編譯選項指導編譯器來進行不同程度的編譯優(yōu)化［26］。FT-M6678 編譯器提供了部分直接影響或控制程序優(yōu)化的編譯選項，通過合理地配置這些優(yōu)化選項，可以達到提升程序運行效率的目的。編譯優(yōu)化是對程序最基礎的優(yōu)化，后續(xù)的優(yōu)化均建立在編譯優(yōu)化之上，其選項的合理配置需要針對程序特點以及實驗嘗試進行驗證。本文基于SYEV 算法的運行特點，經(jīng)過測試，選擇表1 列舉的編譯優(yōu)化選項，以取得較好的加速效果，其中，-o3 會指導編譯器根據(jù)程序語句、關鍵字等有效信息來源了解到循環(huán)執(zhí)行次數(shù)和數(shù)據(jù)存放方式等信息，在編譯階段自動進行循環(huán)展開以及指令流水線等操作，但是效果上因程序而異，在本文研究內(nèi)容中，經(jīng)測試，須配合手動展開才能取得更好的優(yōu)化效果。

表1 編譯優(yōu)化選項 Table 1 Compilation optimization options

2.4 訪存優(yōu)化

2.4.1 緩存優(yōu)化

DSP 的運行耗時主要分布在兩部分：一是計算所需的時間開銷；二是數(shù)據(jù)存取的時間開銷［27］。緩存優(yōu)化是通過對Cache 加速，提高運算數(shù)據(jù)的存取效率。

FT-M6678 每個核內(nèi)存儲都采用了兩級Cache 結(jié)構(gòu)，其中L2 為512 KB，L1 分為一級程序緩存（L1P）和一級數(shù)據(jù)緩存（L1D），容量均為32 KB。此外，還有4 MB 的共享內(nèi)存SMC 以及2 GB 的外部空間DDR3。FT-M6678 的存儲結(jié)構(gòu)與各級存儲器的工作頻率如圖4 所示。

圖4 FT-M6678 多級內(nèi)核存儲結(jié)構(gòu)Fig.4 FT-M6678 multilevel kernel storage structure

FT-M6678 內(nèi)核從L1P 中讀取程序指令，從L1D中加載數(shù)據(jù)，L1D 和L1P 與L2 可直接通信，L2 可 直接與共享內(nèi)存通信，而共享內(nèi)存直接與外部DDR3通信。整體的核內(nèi)通信分級分層劃分明確，且離DSP 內(nèi)核越遠，存儲器工作時鐘越低，即內(nèi)存訪問的速度越慢。

本文在優(yōu)化SYEV 算法中開啟了L1、L2 緩存，設置寄存器與Cache 之間的通信，并將外部存儲器設置為可緩存屬性，在數(shù)據(jù)運算過程中將矩陣數(shù)據(jù)提前加載到緩存中，可在緩存層面優(yōu)化數(shù)據(jù)讀取，減少時間開銷。

2.4.2 訪存優(yōu)化的資源分配

根據(jù)FT-M6678 的硬件特點，對硬件資源（.MEMORY）和段（.SECTIONS）進行分配，可根據(jù)程序段與數(shù)據(jù)段的大小將其分配到合適的位置，提高程序運行時的訪存效率。

首先在硬件資源分配上分別把L1P 和L1D 的部分或全部設置為Cache，容量取最大為32 KB；L2 存儲控制器是L1 存儲器和更高級存儲器的接口，把L2設為Cache，容量也取最大為512 KB；對于共享內(nèi)存SMC 取最大為4 MB；對于外部存儲DDR3，考慮到矩陣大小，取1 GB 空間。

其次是段分配。在text 段中存放程序代碼，將該段放入DDR3 中；對于存放全局變量和常量的段，如cinit 段或const 段等，考慮到此類段數(shù)據(jù)量很小，故放入共享內(nèi)存SMC 中；將sysmem 段用于程序中的函數(shù)動態(tài)分配存儲空間，如malloc 空間等，該段也是堆空間存放的段?？筛鶕?jù)輸入的對稱矩陣的大小將數(shù)據(jù)放入不同的段中：當矩陣規(guī)模不大時，可將sysmem 段放入共享內(nèi)存SMC 中，數(shù)據(jù)無須從外部存儲器再讀入；當矩陣規(guī)模較大時，共享內(nèi)存空間不足以支撐矩陣容量，將sysmem 段放入DDR3 內(nèi)。

2.5 向量并行化優(yōu)化

向量并行化優(yōu)化主要針對耗時占比最高的矩陣平面旋轉(zhuǎn)進行優(yōu)化。首先分析該部分的代碼特點，考慮循環(huán)展開，探究最佳循環(huán)展開次數(shù)，并使用單指令多數(shù)據(jù)流（SIMD）并行操作指令對計算部分以及讀取部分進行向量化操作。

2.5.1 循環(huán)展開

矩陣平面旋轉(zhuǎn)核心操作部分會大量使用雙重for 循環(huán)，每一次循環(huán)實現(xiàn)2 個數(shù)組元素的互換，并且伴隨乘法、加減計算。分別對該核心運算部分進行4、8、12 次循環(huán)展開，選取合適的展開因子，在保證性能的同時也要防止代碼體過大。循環(huán)展開將迭代間并行轉(zhuǎn)為了迭代內(nèi)并行，可以在展開后的循環(huán)體內(nèi)發(fā)掘數(shù)據(jù)級并行，生成向量訪存和運算指令以提高性能。對比這3 次不同的展開運行效率發(fā)現(xiàn)，在1 000 階矩陣的輸入測試下，展開4 次與8 次效果相近，8 次稍優(yōu)，展開12 次性能稍有下降，但整體上相差不大。展開4 次的偽代碼示例如下：

2.5.2 向量化

FT-M6678 系 列支持32 位數(shù) 據(jù)的SIMD 指令，并且允許其在128 位向量上進行操作，如圖5 所示，通過向量化操作將4 個32 位數(shù)據(jù)組成一個128 位的向量，運算時通過SIMD 指令實現(xiàn)128 位的向量乘，從而達到4 對32 位數(shù)據(jù)并行乘操作的效果，這樣做的好處是可以充分利用FT-M6678 的寄存器增大處理器指令調(diào)度的空間，從而提高運算效率。

圖5 FT-M6678 128 位并行乘示意圖Fig.5 Schematic diagram of FT-M6678 128 bit parallel multiplication

FT-M6678 編譯器提供同時支持128 位的并行操作數(shù)據(jù)類型_x128_t，指令示例為FT_fto128；支持的4 個float 并行乘的SIMD 指令示例為FT_qmpysp。FT_fto128 指令將取4 個float 類型的數(shù)據(jù)并定義為一個_x128_t 的向量數(shù)據(jù)放入寄存器中，F(xiàn)T_qmpysp指令則是實現(xiàn)2 個_x128_t 數(shù)據(jù)類型的相乘，即4 對float 數(shù)據(jù)的并行乘運算操作，運算的結(jié)果依然返回一個_x128_t。在本文算法的矩陣平面旋轉(zhuǎn)核心操作上運用該向量化，將耗時的乘法操作并行執(zhí)行。以4 次循環(huán)展開為例，4 次展開的代碼中有16 次指令乘法操作，而向量化將這16 次指令乘法轉(zhuǎn)為4 次指令執(zhí)行，從而提升了運算的效率，以下為4 次展開的向量化偽代碼：

3 實驗結(jié)果與分析

本節(jié)從正確性驗證、優(yōu)化性能分析2 個方面來分析實驗結(jié)果。正確性驗證采用權威線性代數(shù)運算庫LAPACK 官方提供的正確性驗證集，將該驗證程序移植到FT-M6678 平臺上編譯并驗證結(jié)果；優(yōu)化性能分析則生成不同規(guī)模的對稱矩陣作為輸入，測試運行時間，矩陣的生成選用LAPACK 官方提供的生成矩陣的方式，對比方式分為縱向?qū)Ρ群蜋M向?qū)Ρ?，即FT-M6678 優(yōu)化前后的運行速度對比、FT-M6678優(yōu)化后與TMS320C6678 平臺對比，選用的對比平臺TMS320C6678 是DSP 主流的芯片之一，與該平臺對比能凸顯FT-M6678 的可替代性與國產(chǎn)自主性。

3.1 實驗環(huán)境

實驗環(huán)境分別采用目標平臺FT-M6678 與對比平臺TMS320C6678，且設置C6678 的主頻為1.25 GHz，2 個平臺開發(fā)環(huán)境保持一致以便于控制對比變量，實驗在單核運行環(huán)境下進行，具體參數(shù)如表2 所示。

表2 實驗環(huán)境參數(shù) Table 2 Experimental environment parameters

3.2 正確性驗證

驗證集位于LAPACK 文件TESTING 下的EIG文件部分，驗證過程是由MATGEN 功能部分生成測試矩陣，EIG 調(diào)用測試程序進行運算，判斷運算結(jié)果是否正確。由于EIG 部分測試范圍廣泛，測試的是所有有關矩陣特征問題的例程，而對稱矩陣特征求解算法部分的驗證由EIG 下的子程序SEP 測試程序負責，因此只需要執(zhí)行該部分子程序即可。

在FT-M6678 上運行移植后的測試程序，輸入SEP 默認測試參數(shù)，驗證實現(xiàn)且優(yōu)化后的算法，結(jié)果如圖6 所示，在默認給定的不同測試規(guī)模下，正確性測試用例對算法部分進行4 662 次測試，對驅(qū)動例程部分進行14 256 次測試，F(xiàn)T-M6678 平臺均通過測試，驗證了算法的正確性。

圖6 FT-M6678 正確性測試結(jié)果Fig.6 FT-M6678 correctness test results

3.3 優(yōu)化性能分析

分別生 成1 000×1 000、2 000×2 000、3 000×3 000、4 000×4 000、5 000×5 000 規(guī)模的對稱矩陣作為輸入，測試矩陣生成使用LAPACK 官方提供的MATGEN 功能部分，選用生成矩陣函數(shù)LATMS，該函數(shù)可根據(jù)用戶指定的特征值與特征向量生成對應的對稱矩陣，并且該函數(shù)提供種子參數(shù)，輸入LAPACK 提供的不同種子生成的測試矩陣數(shù)值也不同。分別在FT-M6678 以及TMS320C6678 平臺上運行測試矩陣，相同階數(shù)運行3 次以上取結(jié)果平均值，數(shù)據(jù)精確到小數(shù)點后3 位。

本文算法在FT-M6678 上基于不同對稱矩陣規(guī)模的優(yōu)化性能對比結(jié)果如表3 所示。由表3 可知，矩陣規(guī)模在3 000×3 000 后加速比趨于穩(wěn)定，故選取該階數(shù)下的測試數(shù)據(jù)分析單步優(yōu)化性能。

表3 FT-M6678 優(yōu)化前后性能對比 Table 3 Performance comparison of FT-M6678 before and after optimization

表4 列出了單步優(yōu)化分析數(shù)據(jù)，其中數(shù)據(jù)為多次運行取的平均值，誤差范圍保持在毫秒級。表4中的優(yōu)化是一個疊加的過程，即先進行編譯優(yōu)化，在編譯優(yōu)化的基礎上依次進行訪存優(yōu)化和向量并行化。在實驗中嘗試去除單個優(yōu)化后均不能達到最佳效果。以編譯優(yōu)化為例，若單獨去除編譯優(yōu)化，除了編譯優(yōu)化這部分的性能加速會受到影響外，向量并行化優(yōu)化的加速效果也會下降，原因在于編譯優(yōu)化會對算法次熱點的部分進行總體性優(yōu)化，提升了主熱點子程序的耗時占比，若去除編譯優(yōu)化，主熱點子程序的耗時占比會下降，對該部分的向量化優(yōu)化性能也會隨之下降。

表4 單步優(yōu)化數(shù)據(jù) Table 4 Single step optimization data

在FT-M6678 平臺上對1 000×1 000 以上規(guī)模對稱矩陣運行測試后，編譯優(yōu)化加速比約為3.697～3.794，訪存優(yōu)化加速比可達到7.104～7.351，其中主要的性能優(yōu)化來源于Cache 的數(shù)據(jù)預取，開啟緩存優(yōu)化可以很大程度上降低矩陣數(shù)據(jù)的訪存開銷，而向量并行化優(yōu)化可再將速度提升約1.930～2.149 倍。

表5 為算法在FT-M6678 與TMS32C6678 上的對比結(jié)果。結(jié)合表3 分析可知，在優(yōu)化前FT-M6678平臺與TMS320C6678 運行性能差距較大，但在優(yōu)化后FT-M6678 運行性能比TMS320C6678 性能更好。

表5 FT-M6678 與TMS32C6678 性能對比 Table 5 Performance comparison between FT-M6678 and TMS32C6678

圖7為縱向與橫向加速比對比，由該圖可知，本文算法在FT-M6678 上的整體加速效果較好，優(yōu)化前后的加速比可達51.399～58.346，優(yōu)化后較TMS32C6678平臺運行加速比可達1.923～2.053。

圖7 縱向與橫向加速比對比Fig.7 Comparison of longitudinal and lateral acceleration ratios

4 結(jié)束語

本文面向FT-M6678 平臺實現(xiàn)對稱矩陣特征值求解算法并對其進行優(yōu)化。通過分析FT-M6678 的體系結(jié)構(gòu)以及求解算法的原理與過程，在FT-M6678平臺實現(xiàn)該算法，并對算法進行針對性優(yōu)化，充分發(fā)揮FT-M6678 平臺的運算優(yōu)勢，分別使用基于FT-M6678 平臺的編譯優(yōu)化、緩存層面優(yōu)化、基于FT-M6678 和程序運行的硬件資源和段分配優(yōu)化、針對算法的循環(huán)展開以及向量并行化優(yōu)化。對實現(xiàn)且優(yōu)化后的算法進行正確性驗證與優(yōu)化性能分析，經(jīng)測試，正確性驗證全部通過，性能提升明顯，F(xiàn)T-M6678 最終優(yōu)化效果較優(yōu)化前速度提升51.399～58.346 倍，對 比TMS32C6678 平臺速 度提升1.923～2.053 倍，實驗結(jié)果證明FT-M6678 具有可行性，其在性能上可以替代TMS32C6678 平臺，且具有國產(chǎn)自主可控性。后續(xù)可考慮使用OpenMP 多核并行以及對算法底層所調(diào)用的BLAS 函數(shù)進行優(yōu)化，進一步提升優(yōu)化性能。