程 , ,

(上海工程技術(shù)大學(xué) 電子電氣工程學(xué)院,上海 201620)

0 概述

求解大規(guī)模稀疏線性系統(tǒng)是數(shù)值計算中的一個重要問題,廣泛應(yīng)用于力學(xué)、大氣建模、地球物理、生物學(xué)、電路仿真以及計算科學(xué)與工程的其他領(lǐng)域。就大規(guī)模數(shù)值模擬而言,求解大規(guī)模稀疏線性系統(tǒng)Ax=f占據(jù)了大量的計算時間和資源,即需要計算稀疏矩陣A與向量x的乘積,這一類計算問題統(tǒng)稱為稀疏矩陣與向量乘(Sparse Matrix-Vector Multiplication,SpMV)。電磁問題復(fù)雜程度的提高使得計算規(guī)模愈加擴大,如何加快大規(guī)模稀疏線性方程組的求解速度顯得愈加重要,目前,解決該問題的重要方法之一就是運用高性能計算技術(shù)。

目前,常用的求解稀疏線性系統(tǒng)的方法有直接法與迭代法。直接法是通過對方程組的系數(shù)矩陣進(jìn)行變換,將原方程組轉(zhuǎn)化為形如三角矩陣等形式,繼而使用回代或追趕的方法得到方程組的解。迭代法是從解的某個近似值出發(fā),構(gòu)造一個無窮序列去逼近精確解的過程。在實踐中,迭代法受迭代矩陣的影響,其所適用于解的線性方程組類型不盡相同。當(dāng)稀疏矩陣為病態(tài)矩陣時,舍入誤差大大降低了收斂速度,因此,需用各種預(yù)處理技術(shù)加以改進(jìn)以達(dá)到良好的計算效果。本文結(jié)合直接法的可靠性和迭代法的易并行性,采用不完全LU分解預(yù)條件共軛梯度(Incomplete LU factorization preconditioned Conjugate Gradient,ILUCG)法進(jìn)行研究。

文獻(xiàn)[1]提出IC/ILU預(yù)條件CG/BiCGSTAB算法,并采用CUBLAS[2]、CUSPARSE[3]庫對其進(jìn)行實現(xiàn)。文獻(xiàn)[4]提出基于GPU的稀疏線性系統(tǒng)的預(yù)條件共軛梯度法,并采用對角預(yù)處理法實現(xiàn)算法。文獻(xiàn)[5]在文獻(xiàn)[1]和文獻(xiàn)[4]方法的基礎(chǔ)上,提出基于GPU加速的ICCG法。文獻(xiàn)[6]研究基于GPU的矩陣乘法優(yōu)化,提出一種基于HYB(Hybrid)格式的新存儲格式。文獻(xiàn)[7]研究了基于HYB格式稀疏矩陣與向量乘在CPU+GPU[8]異構(gòu)系統(tǒng)中的實現(xiàn)與優(yōu)化。

本文在文獻(xiàn)[6-7]研究的基礎(chǔ)上,提出一種稀疏矩陣的混合存儲格式HEC(Hybrid ELL and CSR),并在CUDA[9]平臺上實現(xiàn)基于HEC格式的ILUCG方法,以測試其加速效果。

1 存儲格式

SpMV的求解是科學(xué)計算中的一個經(jīng)典問題,為提升其運算性能,對稀疏矩陣存儲格式的優(yōu)化顯得非常重要。文獻(xiàn)[10]在CUDA平臺上實現(xiàn)了利用不同的壓縮存儲格式的SpMV運算,并對性能進(jìn)行了比較分析。文獻(xiàn)[11]提出一種新的壓縮矩陣存儲格式CDIA(Compress format DIA gonal),并成功地在CUDA平臺上對該方法進(jìn)行實現(xiàn)。文獻(xiàn)[12]針對CSR存儲格式提出在GPU上的實現(xiàn)方法和優(yōu)化策略。目前,這些對單一存儲格式優(yōu)化的研究取得了較大的進(jìn)展,NVIDIA公司的CUSPARSE線性代數(shù)庫也對這些單一的存儲格式提供了強有力的支持。

然而,由于稀疏矩陣中非零元素分布的不規(guī)則性,導(dǎo)致采用單一的壓縮存儲格式進(jìn)行運算時,取得的效果并不理想,因此多種對矩陣壓縮存儲格式進(jìn)行優(yōu)化的方法被相繼提出:一種是對行進(jìn)行合并分塊來緩解行與行間非零元素差異過大的情況,如BCSR[13]格式、BELLPACK[14-15]格式等,此外,文獻(xiàn)[16]研究了基于GPU的SpMV運算的優(yōu)化,對分段行合并及按行分塊存儲策略進(jìn)行了深刻剖析;另一種是采取混合格式進(jìn)行存儲,以適應(yīng)不同矩陣不規(guī)則的稀疏特征。如針對準(zhǔn)對角矩陣采用的DIA和CSR混合的HDC(Hybrid Dia and CSR)[17]存儲格式以及混合ELL(ELLPACK)和COO(Coordinate)格式的HYB[18]存儲格式等。

目前,HYB存儲格式應(yīng)用尤為廣泛,現(xiàn)有的HYB Kernel函數(shù)大都是在GPU上運行ELL部分,在CPU上運行COO部分,進(jìn)行異構(gòu)計算,考慮到CSR格式存儲的數(shù)據(jù)量較COO格式更少,且GPU具有天然的并行性,因此,本文提出一種矩陣的混合存儲格式HEC,即將矩陣存儲為ELL格式和CSR格式,并在GPU上運行該存儲格式的Kernel函數(shù),以獲得一定的加速效果。

1.1 ELL存儲格式

如圖1所示,ELL格式用2個二維數(shù)組Values[ ]和Col[ ]來存儲一個n×k的矩陣(k為包含非零元素最多行的非零元數(shù)目,這里為n為4,k為2),將一個稀疏矩陣采用類稠密的方式存儲,更適合于并行計算。ELL格式對非零元較為集中于部分行的稀疏矩陣具有良好的存儲效果,二維數(shù)組Values大小為n×k,它存儲原稀疏矩陣中每行非零元的值,若該矩陣的所有行中最大非零元數(shù)為s,每行的非零元總數(shù)為c,則每行(s-c)的部分由0來填充;二維數(shù)組Col大小也為n×k,它用來存儲相應(yīng)非零元素在原系數(shù)矩陣中的列索引,它的相應(yīng)位置也可以用0來填充。

圖1 ELL存儲格式

ELL格式的優(yōu)點是只用2個數(shù)組就可以存儲一個稀疏矩陣,在一定程度上節(jié)省了存儲空間,降低了訪存開銷,且其SpMV算法在加速器內(nèi)易于并行實現(xiàn)。然而,ELL格式不適用于行之間非零元素相差較大的稀疏矩陣,若有一行每一個元素都是非零元,則ELL數(shù)組Values[ ]的列數(shù)和原稀疏矩陣的列數(shù)相等,就完全起不到壓縮的效果。

1.2 CSR存儲格式

如圖2所示,CSR格式采用3個數(shù)組Values[ ]、Rowptr[ ]、Colind[ ]來存儲一個稀疏矩陣,若矩陣的大小為n×n,且其總非零元數(shù)為m,則數(shù)組Values[ ]與Colind[ ]的長度為m,數(shù)組Values[ ]保存原稀疏矩陣中的所有非零元的值,數(shù)組Colind[ ]存儲了各非零元在原稀疏矩陣中的列索引,數(shù)組Rowptr[ ]的長度為n+1,它存儲了每一行元素在壓縮存儲格式中的起始位置。在圖2中,數(shù)組Rowptr[ ]中的元素是每一行起始位置的指針,Rowptr[0]表示第0行,即從數(shù)組Values[ ]的0號位置1開始,Rowptr[1]表示第1行,即從數(shù)組Values[ ]的2號位置3開始,以此類推。

圖2 CSR存儲格式

需要注意的是,若原稀疏矩陣有一行元素全為0(以第1行元素全為0舉例,該矩陣的行編號分別為0、1、2、3行),則Rowptr[1]與Rowptr[2]的值都是2,表示該存儲格式中沒有保存第1行的元素。CSR存儲格式可有效地表示任何形式的稀疏矩陣,因此,其在目前大規(guī)模稀疏線性系統(tǒng)的科學(xué)計算中應(yīng)用尤為廣泛。

1.3 COO存儲格式

如圖3所示,COO格式采用3個數(shù)組Values[ ]、Row[ ]、Colind[ ]來存儲一個稀疏矩陣,它與CSR格式的區(qū)別在于數(shù)組Row[ ]存儲了原稀疏矩陣中所有非零元的行索引。

圖3 COO存儲格式

1.4 HYB存儲格式

如圖4所示,HYB存儲格式由ELL與COO存儲格式混合而成,用2個二維數(shù)組Values和Col來存儲一個n×k的矩陣(k為包含非零元素最多行的非零元數(shù)目),二維數(shù)組Values大小為n×k,它存儲稀疏矩陣中每行非零元的值;二維數(shù)組Col大小也為n×k,它存儲稀疏矩陣中非零元的列索引。

圖4 HYB存儲格式

HYB格式依據(jù)閾值K將矩陣劃分為ELL格式與COO格式存儲,即將數(shù)組Values第K(閾值)列左邊的部分(包括K列)存儲為ELL格式,右邊的部分存儲為COO格式。對于此,產(chǎn)生了一個重要的問題是閾值K的劃分,即如何劃分可以獲取最優(yōu)的性能,這里采取界限法,劃分規(guī)則如下:

先求得ELL格式存儲的矩陣,在這之中,若第K列的元素加到劃分的ELL部分后,若ELL部分的總非零元數(shù)超過其總元素的二分之一,那么第K列將會被劃到ELL部分。

閾值K劃分算法描述如下:

求出第1列的非零元數(shù)

for j=2:1:sf

//j為矩陣Values的列,sf為矩陣Values的總列數(shù)

求出前j列的總非零元數(shù)s

求出前j列的元素的總數(shù)量c

IF s>(c/2),則 K=j;ELSE break

求出閾值K

例如,在圖4中,在將稀疏矩陣的第1列劃分為ELL格式后(共有4列,編號為1、2、3、4),若將第2列劃分為ELL格式,則ELL部分元素的總數(shù)量為10,總非零元數(shù)為6>5,因此,第2列應(yīng)該劃分為ELL格式;若將第3列劃也為ELL部分,則總元素為15,總非零元數(shù)為7<7.5。因此,第3列不應(yīng)該被劃分為ELL部分,即此時的閾值K應(yīng)取2。

1.5 HEC存儲格式

如圖5所示,HEC格式由ELL與CSR存儲格式混合而成,根據(jù)閾值K將矩陣劃分為ELL格式與CSR格式兩部分存儲,閾值K的劃分規(guī)則與HYB格式相同,即K仍取2。ELL部分的內(nèi)容仍與HYB格式相應(yīng)的部分相同,遍歷剩余部分所有非零元的所在行,將其劃分到CSR部分,由于剩余部分的3個數(shù)2、3、1位于同一行,因此CSR部分的數(shù)組Rowptr只需存儲2個數(shù)0、3,即為起始位置索引0和劃分到CSR部分元素的個數(shù)3。

圖5 HEC存儲格式

若存儲為HYB格式,則COO部分的數(shù)組Row需存儲 3個數(shù)0、0、0,即為原稀疏矩陣中相應(yīng)部分各非零元的行索引。在本例中,稀疏矩陣的行列數(shù)及非零元數(shù)較少,相較于HYB格式只少存儲了1個數(shù),但在大規(guī)模稀疏線性系統(tǒng)的求解中,隨著稀疏矩陣行列數(shù)及非零元數(shù)的增多,HEC存儲格式的效率將會顯著提升。

HEC格式保留了ELL格式適合于并行算法的優(yōu)勢,并且它可克服ELL格式只適合于矩陣各行間非零元數(shù)相差不大情況的缺點,從而減少填充元素的數(shù)量,因而其具有良好的性能。

此外,HEC存儲格式比CSR格式多存儲了2組數(shù)據(jù)Val_ELL和Col_ELL(分別為以ELL格式存儲的值和列),會對讀取帶來一定的開銷,但在真實應(yīng)用中,SpMV操作會在下一個迭代求解器中對矩陣A反復(fù)進(jìn)行。在每一次迭代求解時,向量x和向量y(y=Ax)會發(fā)生改變,但矩陣A不變,因此,生成HEC存儲格式的工作可攤銷在多次的迭代求解中。

2 不完全LU分解的預(yù)條件共軛梯度法

2.1 預(yù)條件共軛梯度法

本文采取的預(yù)條件共軛梯度法是通過采取一系列的預(yù)處理技術(shù)減少共軛梯度法中的迭代次數(shù),從而加速其收斂的一種方法。

共軛梯度法[19]是迭代法的一種,迭代法采用逐次逼近并用迭代公式得到一系列近似解,從而逐漸地逼近于真實解,最終可得滿足一定誤差條件的近似解。共軛梯度法是依賴于向量正交性的,對舍入誤差的影響十分敏感,其收斂速度和系數(shù)矩陣的條件數(shù)緊密相關(guān)。理論上,它用小于n(矩陣的大小為n×n)的迭代次數(shù)即可得到求解結(jié)果。但當(dāng)稀疏矩陣呈現(xiàn)病態(tài)時,在舍入誤差的影響下,迭代次數(shù)大于n,收斂速度較慢,會出現(xiàn)數(shù)值解與理論解嚴(yán)重偏離的現(xiàn)象。

文獻(xiàn)[20]對共軛梯度法進(jìn)行分塊處理,最終得到解決最小二乘問題的算法。

2.2 不完全LU分解

不完全LU分解是一種相當(dāng)高效的預(yù)處理技術(shù)。LU分解是矩陣分解的一種,可將一個矩陣分解為一個下三角矩陣和上三角矩陣的乘積,或為上、下三角矩陣和一個置換矩陣的乘積。對式(1),若n階稀疏矩陣A對稱且正定,則如式(2)所示,可確定一個正定矩陣M,在求解向量x時可比式(1)獲得更好的性能,矩陣M稱為預(yù)處理矩陣,構(gòu)造預(yù)處理矩陣的方法[21]有多種,如對角預(yù)處理構(gòu)造法、多項式預(yù)處理構(gòu)造法和不完全LU分解預(yù)處理法,本文正是采用了不完全LU分解預(yù)處理法。

在式(1)中,對n階稀疏矩陣A進(jìn)行分裂式近似分解,可得式(3),其中,L為下三角矩陣,P為剩余矩陣。式(2)為預(yù)處理后的方程組,矩陣M與矩陣A具有相似的稀疏性,因此,為提高算法的求解速率,用矩陣M代替矩陣A進(jìn)行迭代運算。若將矩陣M進(jìn)行LU分解,可得式(4),其中,LM為下三角矩陣,UM為上三角矩陣。

Ax=f

(1)

M-1Ax=M-1f

(2)

A=LLT+P=M+P

(3)

A≈M=LMUM

(4)

2.3 ILUCG算法

基于LU分解的ILUCG算法描述如下:

設(shè)定初值x0

r=f-Ax0,p=r,rw=r

for i=0,1,2,3,…,maxit

ρi=rwT×r

if ρi==0.0 then?

method failed

end if

求解 M×ph=p,得到ph

計算q,q=A×ph

s=r-α×q

x=x+α×ph

if ‖s‖2≤tol then

method converged

endif

求解 M×sh=s,得到sh

計算 t=A×sh

求得 x=x+ω×sh

求得 r=s-ω×t

end for

3 算法實現(xiàn)

3.1 基于CUSPARSE和CUBLAS庫的ILUCG法

利用CUDA中的CUSPARSE和CUBLAS庫可以對ILUCG算法進(jìn)行實現(xiàn),在運算過程中,GPU負(fù)責(zé)迭代前和每次迭代過程中的矩陣與向量,向量與向量間的并行計算,CPU負(fù)責(zé)迭代循環(huán)和收斂條件的控制以及標(biāo)量相除等操作。算法描述如下:

初始化CUSPARSE和CUBLAS庫

對矩陣A進(jìn)行不完全LU分解,得到一個下三角矩陣L和上三角矩陣U,分別讀取矩陣A、L、U的值、行索引和列,分別為(Data,Row_ptr,Col_index)、(data_L,row_ptrL,col_indexL)、(data_U,row_ptrU,col_indexU)

創(chuàng)建descr_L、descr_U、descr_A、info_L、info_U并初始化

用cudaMalloc函數(shù)給各變量分配內(nèi)存

用cudaMemcpy進(jìn)行數(shù)據(jù)傳輸

用cusparseDcsrsv_analysis函數(shù)生成info_L、info_U

用cusparseDcsrmv函數(shù)計算r=A×x0

用cusparseDcsr2hyb函數(shù)將矩陣A以HYB格式保存在hybA中

用cusparseDhybmv函數(shù)計算r=A×x0//以HYB格式計算

用cublasDscal函數(shù)計算r=-r

用cublasDaxpy函數(shù)計算r=b+r=-A×x0+b

用cublasDcopy函數(shù)賦值p=r,rw=r

用cublasDnrm2計算nrmr0=(r12+r22+…+rN2)0.5

for(i=0;i

rhop=rho,進(jìn)行rho值的更新,使得rhop[i]=rho[i-1]

用cublasDdot函數(shù)計算內(nèi)積(r,rw)

if (i>0) then

β=(rho/rhop)×(α/ω)

用cublasDaxpy、cublasDscal函數(shù)計算p=r+β×(p-ω×q)

end if

用cusparseDcsrsv_solve函數(shù)計算ph,ph=M-1×p

A.用cusparseDcsrmv函數(shù)計算q=A×ph

B.用cusparseDhybmv函數(shù)計算q=A×ph

用cublasDdot函數(shù)計算內(nèi)積(r,q),即temp=rT×q

計算a=rho/temp

用cublasDaxpy函數(shù)計算s=r-a×q,x=x+a×ph

用cublasDnrm2計算nrmr=(r12+r22+...+rN2)0.5

若(nrmr/nrmr0)

用cusparseDcsrsv_solve函數(shù)計算sh,sh=M-1×s

用cusparseDcsrmv函數(shù)計算t=A×sh

用cusparseDhybmv函數(shù)計算t=A×sh

用cublasDaxpy函數(shù)計算x,x=x+w×sh,r=s-w×t

用cublasDnrm2計算nrmr=(r12+r22+...+rN2)0.5

If (nrmr/nrmr0)

//tol為自設(shè)允許的誤差范圍,這里取10-6

end for

3.2 基于CUDA的SpMV算法

3.2.1 基于CSR存儲格式的SpMV算法

在不完全Cholesky分解預(yù)條件共軛梯度算法中SpMV運算的操作可由自作的GPU Kernel函數(shù)完成,使用cudaStream_t(CUDA流)以提高程序的執(zhí)行效率。若將矩陣A以CSR格式存儲,則有Kernel函數(shù)如下:

__global__void SpMV_CSR(double *data ,int *col_index,int *row_ptr,double *x,double *y)

{

int elem;

int row =blockIdx.x*blockDim.x + threadIdx.x;

if(row< N) //N為矩陣A的行數(shù)

{

float dot=0;

int row_start=row_ptr[row];

int row_end=row_ptr[row+1];

for(elem=row_start;elem

{

dot+=data[elem]* x[col_index[elem]];

}

y[row]=dot;

}

}

3.2.2 基于HEC存儲格式的SpMV算法

若將稀疏矩陣A以HEC格式存儲,可調(diào)用SpMV_ELL和SpMV_CSR 函數(shù)分別計算劃分到ELL和CSR部分的值并保存在數(shù)組y[row]和y2[row]中,最終,將y[row]與y2[row]的相應(yīng)值相加,可得一次SpMV運算的結(jié)果。

基于HEC存儲格式的SpMV算法如下:

__global__voidSpMV_ELL(double *data ,int *col_index,double *x,double *y)//GPU ELL Kernel函數(shù)

{

int i;

int row =blockIdx.x*blockDim.x + threadIdx.x;

if(row< P)//P為所劃分出的ELL格式的總行數(shù)

{

float dot=0;

for(i=0;i

{

dot+=data[row+i*P]* x[col_index[row+i*P]];

}

y[row]=dot;

}

}

__global__ voidSpMV_CSR(double *data2 ,int *col_index2,int *row_ptr,double *x,double *y2) //GPU CSR Kernel函數(shù)

{

int elem;

int row2 =blockIdx.x*blockDim.x + threadIdx.x;

if(row2< P)

{

float dot=0;

int row_start=row_ptr[row2];

int row_end=row_ptr[row2+1];

for(elem=row_start;elem

{

dot+=data2[elem]* x[col_index2[elem]];

}

y2[row2]=dot;

}

}

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

本文實驗的計算平臺,其CPU為Intel Core(TM) i5-6500 @3.20 GHz;GPU為NVIDIA GeForce GTX1050;運行內(nèi)存(RAM)的大小為8 GB;系統(tǒng)為Windows7 64位;實驗環(huán)境為Visual Studio 2012,CUDA8.0。測試所用稀疏矩陣全部來自UF Sparse Matrix Collection[22],即佛羅里達(dá)大學(xué)的稀疏矩陣集,如表1所示。

表1 測試所用稀疏矩陣

4.1 總時間對比

分別運用調(diào)用CUSPARSE庫函數(shù)cusparse Dcsrmv()、cusparseDhybmv()和調(diào)用GPU Kernel SpMV_CSR、SpMV_HEC的方式進(jìn)行SpMV運算,并調(diào)用計時函數(shù)clock_t()記錄下基于LU分解的ILUCG法的總運算時間,運行50次,取平均值。其中:

1)在調(diào)用cusparseDhybmv()函數(shù)進(jìn)行SpMV運算前,需先調(diào)用cusparseDcsr2hyb()函數(shù)將CSR格式存儲的稀疏矩陣信息轉(zhuǎn)換為HYB格式進(jìn)行存儲并保存在hybA中。

2)總時間=CPU向GPU數(shù)據(jù)傳輸時間+GPU計算時間。

上述指標(biāo)對比如圖6所示。其中,LU_CSR、LU_HYB分別代表以CSR、HYB為稀疏矩陣的存儲格式并調(diào)用CUSPARSE庫函數(shù)進(jìn)行運算的實現(xiàn)方式。LU_K_CSR、LU_K_HEC分別代表以CSR、HEC為存儲格式并調(diào)用GPU Kernel進(jìn)行運算的實現(xiàn)方式。

圖6 4種實現(xiàn)方式的總時間對比

由圖6可知,LU_K_HEC的實現(xiàn)方式對本次實驗所測試的矩陣具有良好的性能,當(dāng)稀疏矩陣為HB/bcsstk29 時,以該實現(xiàn)方式在4種方式中耗時最短,可獲得最優(yōu)的加速效果,為10.4%。一方面,這是由于矩陣HB/bcsstk29具有對角特征,行中的非零元素較為均衡,分布在對角線之外的元素稀少,分割出來的ELL矩陣進(jìn)行數(shù)據(jù)(0)的填充較少,因而具有較好的加速效果;另一方面,相較于HYB的存儲方式,該方式存儲的非零元更少,更利于其發(fā)揮在內(nèi)存讀取與數(shù)據(jù)傳輸方面的優(yōu)勢;相較于CSR的存儲方式,矩陣 HB/bcsstk29有較多的非零元,可抵消函數(shù)調(diào)用帶來的損耗,更利于GPU的并行化處理,從而縮短GPU的計算時間。

4.2 GPU計算時間對比

圖7展示了GPU計算時間的對比情況,GPU計算時間即運用GPU處理整個預(yù)條件共軛梯度法的運算時間,計時方法采用調(diào)用CUDA函數(shù)cudaEvent_t()進(jìn)行計時。

圖7 4種實現(xiàn)方式的GPU計算時間對比

從圖7可以看出,對于多數(shù)矩陣,LU_K_CSR的實現(xiàn)方式耗時最短,這是由于這些矩陣的非零元數(shù)較少,以CSR格式進(jìn)行存儲相較于HEC格式進(jìn)行存儲,只需調(diào)用一個函數(shù),減少了計算損耗。對于矩陣HB/bcsstk33,它比HB/bcsstk28有更多的非零元,花費的計算時間卻相近,原因是該矩陣經(jīng)LU分解后所得到的矩陣L和U有更少的非零元數(shù),減少了GPU的計算時間。而對于矩陣HB/bcsstk29,LU_K_HEC的實現(xiàn)方式耗時最短,表明HEC存儲格式對大型稀疏線性矩陣有良好的存儲效率。

5 結(jié)束語

本文提出一種稀疏矩陣的存儲格式HEC,并將其以ILUCG法中調(diào)用GPU Kernel的方式應(yīng)用于大型稀疏線性系統(tǒng)的求解中,與利用CUSPARSE庫函數(shù)cusparseDcsrmv()、cusparseDhybmv()的實現(xiàn)方式相比,最高可得到10.4%的加速效果。

在應(yīng)用領(lǐng)域,矢量有限元法是求解電磁問題的重要方法之一,該方法將一個離散后的線性系統(tǒng)歸結(jié)為方程Ax=b的求解,而電磁問題的復(fù)雜化使得系數(shù)矩陣A愈加不規(guī)則,因而可以考慮通過多次劃分ELL格式進(jìn)行存儲,以用于更加不規(guī)則的線性方程組計算。同時,建立一個稀疏矩陣存儲格式的優(yōu)選模型,根據(jù)不同矩陣的特點選取最優(yōu)的存儲格式以提升計算效率,也將是下一步的研究方向。