面向FT-M6678 的對稱矩陣特征值求解算法實現與優化

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

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

0 引言

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

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

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

1 FT-M6678 體系結構與算法原理

1.1 FT-M6678 體系結構

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

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

1.2 SYEV 算法原理

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

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

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

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

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

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

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

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

2 面向FT-M6678 的SYEV 實現與優化

2.1 面向FT-M6678 的SYEV 算法實現

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

圖2 SYEV 算法實現過程Fig.2 Implementation process of SYEV algorithm

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

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

2.2 熱點分析

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

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

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

2.3 編譯優化

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

表1 編譯優化選項 Table 1 Compilation optimization options

2.4 訪存優化

2.4.1 緩存優化

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

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

圖4 FT-M6678 多級內核存儲結構Fig.4 FT-M6678 multilevel kernel storage structure

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

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

2.4.2 訪存優化的資源分配

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

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

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

2.5 向量并行化優化

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

2.5.1 循環展開

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

2.5.2 向量化

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

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

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

3 實驗結果與分析

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

3.1 實驗環境

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

表2 實驗環境參數 Table 2 Experimental environment parameters

3.2 正確性驗證

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

在FT-M6678 上運行移植后的測試程序，輸入SEP 默認測試參數，驗證實現且優化后的算法，結果如圖6 所示，在默認給定的不同測試規模下，正確性測試用例對算法部分進行4 662 次測試，對驅動例程部分進行14 256 次測試，FT-M6678 平臺均通過測試，驗證了算法的正確性。

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

3.3 優化性能分析

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

本文算法在FT-M6678 上基于不同對稱矩陣規模的優化性能對比結果如表3 所示。由表3 可知，矩陣規模在3 000×3 000 后加速比趨于穩定，故選取該階數下的測試數據分析單步優化性能。

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

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

表4 單步優化數據 Table 4 Single step optimization data

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

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

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

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

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

4 結束語

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