基于CUDA的大型實對稱矩陣并行求逆算法

霍迎秋，王武星，彭楚風，方 勇

（西北農林科技大學 信息工程學院，陜西 楊凌712100）

0 引 言

矩陣求逆被大量應用于衛星導航定位［1］、雷達脈沖壓縮［2］、圖像處理［3］及預測控制［4］等許多工程技術領域。目前主要的矩陣求逆算法如QR分解法［5，6］、Jacobi數值法［7］、LU 分解法［8，9］、極小多項式求逆［10］等算法對大型實對稱矩陣求逆時，運行速度慢，不能滿足實際工程中對實時性要求高的需求。因此，針對大型矩陣設計并行求逆算法，在保證精度滿足要求的前提下，提高算法的運行速度，在實際工程具有重要的實際應用意義。

1 分塊迭代求逆算法

對大型實對稱矩陣求逆，可用分塊迭代法實現，其原理如下：

將一個n階實對稱矩陣A ∈Rnxn分塊，設分塊形式為

式中：Wt－1——矩陣A 的前t－1階方陣，——A 矩陣的第t行的前t－1個元素，rt——的轉置，ρ——A 矩陣第n行第n 列元素。

當t＝n時，W－1t即為矩陣A 的逆矩陣。

基于CPU 實現串行分塊迭代求逆算法，流程如圖1所示。

圖1 串行分塊迭代求逆算法流程

2 CUDA架構簡介

CUDA 是由NVIDIA 在2007 年2月推出的一種新的GPU 架構，使人們可以更加方便的使用GPU 對其所進行的工作加速，其開發語言與C、C＋＋非常類似，且不用借助圖形API等額外操作。CUDA 優異的性能，使得它被推出之后廣泛應用于模式識別［11］、視頻檢測［12］、壓縮感知［13］、信號處理［14］、流體力學［15］等領域。

支持CUDA 的GPU （graphic processing unit）與傳統的GPU 最大區別在于它擁有片內共享存儲空間用以減小片外訪問延遲［16］。CUDA 源程序由在主機端 （CPU）上運行的控制程序和運行于設備端 （GPU）上的核心計算程序（kernel）兩部分組成，kernel在一個線程格 （Grid）上執行，每個線程格由大小相同的一組線程塊 （Block）組成，同一線程塊內的所有線程可以共享存儲空間用來共同協作完成計算任務，線程塊之間不能共享存儲空間。每個線程塊在一個SM （streaming multiprocessor）上運行。為了管理同時運行的數百個線程，SM 利用了一種稱為SIMT（single instruction multiple thread）的新架構［16，17］，各線程被SM 映射到一個TP （thread processor）核心上獨立執行。SIMT 以32個線程為一組進行調度，這種線程組被稱為warp塊。

CUDA 硬件模型有多個SM，每個SM 又有若干個SP（stream processor）和若干個32－bit的寄存器。此外，模型中還包括所有線程可訪問的常量內存和紋理內存，以加速對數據的讀取。CUDA 硬件模型如圖2所示［17］。

圖2 CUDA 硬件模型

3 基于CUDA的并行分塊迭代求逆算法

3.1 分塊迭代求逆算法并行性分析

分塊迭代求逆算法是一個循環迭代計算的算法，各循環之間有著強烈的依賴關系，因此算法不能完全展開使之在GPU 上運行。但經仔細分析發現，每個循環中包含大量的向量與向量相乘、矩陣與向量相乘、向量數乘和矩陣相加等運算，而迭代算法主要的計算量集中在這些運算之上。并且矩陣運算中各元素獨立計算，互相之間沒有依賴關系，非常適合基于GPU 對算法進行并行優化設計。

3.2 迭代求逆算法的并行設計

針對迭代算法每個循環中，行向量與列向量相乘、列向量與行向量相乘、矩陣與向量相乘、矩陣相加和矩陣與向量相乘等部分進行并行優化設計。求逆算法整體流程不變，基于CUDA 的并行算法流程如圖3所示。

3.2.1 矩陣與向量相乘

矩陣與相量相乘結果為向量。結果向量的每個元素由矩陣的某行與向量對應項相乘結果之和。設計kernel函數：每個線程塊計算矩陣的某行與向量對應項相乘，然后進行“塊內歸約求和”，其結果即為結果向量的一個元素值。主要設計思路如圖4所示。

歸約求和：若對數組array歸約求和，則每個線程將array的兩個元素相加，再將結果保存回array。每個線程都將array中的兩個元素合并為一個元素，此步驟完成之后，結果的數量就為原來的一半。重復執行上述操作，直至將array中所有元素歸約為一個值。歸約求和方法如圖5所示。

3.2.2 向量相乘的并行設計

算法中涉及到的向量相乘分為行向量與列向量相乘、列向量與行向量相乘，分別對這兩種向量相乘進行并行設計。

（1）行向量與列向量相乘

圖3 并行算法流程

圖4 矩陣與向量相乘

計算行向量與列向量相乘。首先，計算兩向量對應元素乘積，然后對所有乘積求和。設計核函數，計算兩向量對應項乘積時為kernel函數分配32個一維線程塊，每個線程塊分配256個線程。每個線程塊計算兩向量256個對應項元素相乘。并進行 “塊內歸約求和”，將求和結果寫入顯存，同步所有線程塊，在顯存上進行 “二次歸約求和”。求和結果即為兩向量相乘的結果。設計思路如圖6所示。

“塊內歸約求和”與 “二次歸約求和”設計思路相同，只是計算操作的區域不同，“塊內歸約求和”是在線程塊的共享內存中進行，“二次歸約求和”是在顯存上進行。歸約方法如圖5所示。

（2）列向量與行向量相乘

圖5 歸約求和

圖6 行向量與列向量相乘

列向量與行向量相乘，其結果為一個矩陣，矩陣的每個元素為兩向量對應元素的乘積。設計核函數：為kernel函數分配二維線程塊，每個維度大小為 （向量長度n ＋15）／16，每個線程塊分配16×16 個線程，線程的索引即為結果矩陣對應元素的索引，使得每個線程計算兩個對應項元素相乘，同時將計算結果寫入顯存。設計思路如圖7所示。

圖7 列向量與行向量相乘

3.2.3 向量的數乘

向量的數乘即為一個數字和某向量的乘積。設計核函數：分配線程個數不少于向量長度，每個線程計算數字與一個向量元素的乘積，將結果寫入顯存。具體設計方法類似于圖6歸約求和之前部分。

3.2.4 矩陣相加

矩陣相加即為兩矩陣相同索引處兩元素之和。設計核函數：分配 （m，m）個線程塊，每個線程塊分配 （16，16）個線程 （m×16不小于兩相加矩陣的任一維度），則對應線程的索引即為兩矩陣對應元素的索引，使得每個線程計算兩矩陣對應元素之和，將求和結果寫入顯存。設計思路如圖8所示。

圖8 矩陣相加

4 實驗結果與分析

為對比分析并行分塊迭代求逆算法的性能，本文基于C語言實現串行算法，基于CUDA 架構設計并行優化算法。硬件環境：Intel（R）Core（TM）2Quad CPU Q8400，內存4GB，GPU 采用NVIDIA （英偉達）GPU 處理器GeForce GTX 780，其主要參數：12個SM，每個SM 有192個SP，共享內存48K，顯存3GB。軟件環境：WIN7旗艦版32位操作系統，Microsoft Visual Studio 2010。

設計一個隨機生成實對稱矩陣函數，生成指定階數的實對稱矩陣，分別用串行求逆算法和并行優化算法對同一個實對稱矩陣進行求逆運算，利用C 語言的clock（）函數對串行求逆算法計時 （單位：ms），利用CUDA 事件對并行優化算法計時 （單位：ms）。計算出并行優化算法相對于串行求逆算法的加速比，計算兩種算法求逆結果的均方差，以此判斷并行優化算法計算的相對誤差。實驗結果如圖9、圖10所示。

圖9 加速比

通過對圖9的實驗結果分析發現，并行優化算法的運行速度明顯要比串行迭代求逆算法的運行速度快，在矩陣階數為8000×8000 時，串行分塊迭代求逆算法需要運行88.50分鐘，而并行優算法只需要19.03s，與串行求逆算法相比，其加速比可以達到279，且隨著矩陣階數的增加，其加速比也越來越大，可以很好的滿足實際工程中對實時性的要求。

圖10 均方差

通過對圖10實驗結果分析發現：本文所設計的并行優化算法所用浮點型數據采用單精度浮點型表示時，其與串行求逆算法兩者計算結果的均方差在10－8級別，均方差的大小也并未隨著矩陣階數的增加而增加。當矩陣階數更大時，本文所設計并行求逆算法與串行求逆算法計算結果的均方差大概穩定在10－8級別，此精度完全可以滿足工程的實際需要。

綜上所述，本文所設計的并行優化算法與串行求逆算法相比，雖然有些許誤差，但是可以滿足實際工程對精度的要求，而且大大提高了求逆算法的運行速度，隨著矩陣階數的增加，加速比也越來越大。因此，此并行優化算法具有很高的實際應用價值。

5 結束語

本文簡要說明了分塊迭代求逆算法的原理，給出了串行求逆算法和并行優化算法的流程圖，同時對并行優化算法做了詳細的介紹，并基于C 語言實現了串行分塊迭代求逆算法，基于CUDA 架構設計了并行優化算法。對10個隨機產生的大型實對稱矩陣進行求逆運算，統計了兩種算法求逆結果的均方差和并行優化算法相對于串行求逆算法的加速比。實驗結果表明，并行優化算法表現出了優異的性能，特別算法是運算時間得到了極大的提高，與串行算法相比，在矩陣階數為8000×8000 時，其加速比高達279，并且隨著矩陣階數的增大，加速比越來越大；而兩種算法求逆結果的均方差表明，本文所設計的并行分優化算法可以滿足工程實際中對精度的要求。因此，本文設計的并行優化算法在工程實際中具有重要的應用價值。

［1］ZENG Degui.Matrix inversion and its application in beidou double－star positioning system ［J］.Information and Computer，2010 （9）：31－32 （in Chinese）. ［曾德貴.矩陣求逆及其在北斗雙星定位系統上的應用 ［J］.信息與電腦，2010 （9）：31－32.］

［2］SONG Yixin，YAO Zhendong.RMMSE radar pulse compression fast algorithm of matrix inversion of the FPGA implementation ［J］.Journal of Chengdu Information Engineering College，2009，24 （5）：435－439 （in Chinese）. ［宋一鑫，姚振東.RMMSE雷達脈沖壓縮快速算法中矩陣求逆的FPGA 實現［J］.成都信息工程學院學報，2009，24 （5）：435－439.］

［3］ZHENG Zuoyong，ZHANG Ruixia.A quick method for inversing a circulant matrix on GPU ［J］.Engineering and Computer Science，2012，34 （7）：84－88 （in Chinese）. ［鄭作勇，張瑞霞.GPU 上循環的快速求逆算法 ［J］.計算機工程與科學，2012，34 （7）：84－88.］

［4］SONG Jian.The implementation of the predictive control algorithm based on FPGA ［D］.Dalian：Dalian Maritime University，2010 （in Chinese）. ［宋見.基于FPGA 的預測控制算法的實現 ［D］.大連：大連海事大學，2010.］

［5］NI Tao，DING Haifeng，RUAN Liting，et al.Implementation of arbitrary order complex matrix inversion based on the QR decomposition algorithm in DSP ［J］.Electronic Science and Technology，2010，23 （4）：99－101 （in Chinese）.［倪濤，丁海峰，阮黎婷，等.基于QR 分解算法的任意階復矩陣求逆的DSP實現 ［J］.電子科技，2010，23 （4）：99－101.］

［6］HAN Liangliang，YE Ping.A Jacobi matrix inversion method for redundant robotic manipulator base on QR decomposition［J］.Software，2013，34 （11）：64－66 （in Chinese）. ［韓亮亮，葉平.基于QR 分解的冗余度機械臂雅可比矩陣求逆方法［J］.軟件，2013，34 （11）：64－66.］

［7］WANG Jiexian，FENG Baoxin.Application of Jacobi numerical method in inverse problem of real symmetry matrix ［J］.Geodesy and Geodynamics，2010，30 （1）：74－76 （in Chinese）.［王解先，馮寶新.Jacobi數值方法在實矩陣求逆中的應用 ［J］.大地測量與地球動力學，2010，30 （1）：74－76.］

［8］LI Tao，ZHANG Zhongpei.Matrix inversion by FPGA ［J］.Communication Technology，2010，43 （11）：147－149 （in Chinese）.［李濤，張忠培.矩陣求逆的FPGA 實現 ［J］.通信技術，2010，43 （11）：147－149.］

［9］ZHAO Liqun.The inverse and determinant computation of some sparse matrix ［D］.Zhangzhou：Zhangzhou Normal College，2011 （in Chinese）.［趙立群.一些稀疏矩陣的逆和行列式的計算 ［D］.漳州：漳州師范學院，2011.］

［10］YIN Yunxing，ZHAO Qing.Minimal polynomial in the application of matrix inversion ［J］.Science，Technology and Engineering，2009，9 （5）：1217－1218 （in Chinese）.［殷云星，趙清.極小多項式在矩陣求逆中的應用 ［J］.科學技術與工程，2009，9 （5）：1217－1218.］

［11］ZHANG Shu.Pattern recognition parallel algorithm and the GPU implementation study at a high speed ［D］.Chengdu：University of Electronic Science and Technology，2009 （in Chinese）.［張舒.模式識別并行算法與GPU 高速實現研究［D］.成都：電子科技大學，2009.］

［12］FU Cheng.Research on video detection algorithm for vehicles red－light running based on GPU ［D］.Chengdu：Xihua University，2012 （in Chinese）. ［付誠.基于GPU 的闖紅燈車輛視頻檢測算法研究 ［D］.成都：西華大學，2012.］

［13］Yong Fang，Liang Chen.GPU implementation of orthogonal matching pursuit for compressive sensing ［C］／／Proc IEEE ICPADS，2011：1044－1047.

［14］ZHAO Fangbin.Research on passive radar signal processing platform based on FPGA and GPU ［D］.Chengdu：University of Electronic Science and Technology，2013 （in Chinese）.［趙芳斌.基于FPGA 和GPU 的外輻射源雷達信號處理平臺的研究 ［D］.成都：電子科技大學，2013.］

［15］DONG Yanxing，LI Xinliang，LI Sen，et al.Acceleration of computational fluid dynamics codes on GPU ［J］.The Computer System Application，2011，20 （1）：104－109 （in Chinese）.［董延星，李新亮，李森，等.GPU 上計算流體力學的加速 ［J］.計算機系統應用，2011，20 （1）：104－109.］

［16］GAN Xinbiao，SHEN Li，WANG Zhiying.Parallelizing full search algorithm for motion estimation using CUDA ［J］.Journal of Computer Graphics，2010，22 （3）：457－460 （in Chinese）.［甘新標，沈立，王志英.基于CUDA 的并行全搜索并行估計算法 ［J］.計算機圖形學學報，2010，22 （3）：457－460.］

［17］NVIDIA.CUDA programming guide 2.0 ［M］.Santa Clara：NVIDIA Corporation，2008.