董延亮

【摘要】稀疏方程組一般都是大型系數(shù)矩陣方程組，如果按照直接法求解，往往需要耗費(fèi)大量的內(nèi)存資源和計(jì)算時(shí)間.本文則給出了一種使用變帶寬壓縮存儲(chǔ)的方式求解對(duì)稱正定稀疏方程組的高效解法.

【關(guān)鍵詞】變寬帶存儲(chǔ)；稀疏方程組；解法

稀疏方程組一般都是大型系數(shù)矩陣方程組，如果按照直接解法的一般原理，往往需要耗費(fèi)大量的內(nèi)存資源和計(jì)算時(shí)間.但是對(duì)于一些特殊形式的系數(shù)矩陣，我們可以根據(jù)其自身特點(diǎn)尋找相應(yīng)的矩陣存儲(chǔ)形式，進(jìn)而得到一種更為節(jié)省內(nèi)存、提高運(yùn)算速度的新解法.這里我們以對(duì)稱正定稀疏方程組為例，應(yīng)用變帶寬存儲(chǔ)的方法壓縮儲(chǔ)存系數(shù)矩陣，給出相應(yīng)方程組的高效算法.

下用變帶寬壓縮存儲(chǔ)法求解對(duì)稱正定稀疏線性方程組

AX=B，（1-1）

其中A∈Rn×n是對(duì)稱正定稀疏陣，X，B∈Rn×m（m≥1）.

先對(duì)A，再對(duì)B分別作如下的分解：

A=LDLT，（1-2）

A=LDB～.（1-3）

其中L為單位下三角矩陣（lii=1，i=1，2，…，n，lij=0，i

設(shè)矩陣A各行的第一個(gè)非零元素是aimi，i=1，2，…，n.記mij=max（mi，mj），則由式（1-2）和（1-3），根據(jù)矩陣乘法規(guī)則，推出矩陣L，D，B～的元素的下列計(jì)算公式.

lij=（aij-∑j-1k=mijlikdkljk）/dj，j=mi，mi+1，…，i-1，i=2，3，…，n，

di=aii-∑i-1k=mil2ikdk，i=1，2，…，n，

b～ik=（bik-∑i-1j=milijdjb～jk）/di，i=1，2，…，n，k=1，2，…，m.

將式（1-2）和（1-3）代入式（1-1），得到

LTX=B～.（1-4）

這樣，方程組（1-1）的求解就歸結(jié)為方程組（1-4）的求解.

以下給出對(duì)稱正定稀疏方程組（1-1）的算法，其中A∈Rn×n是對(duì)稱正定稀疏陣，X，B∈Rn×m（m≥1）.

1.采用變帶寬壓縮存儲(chǔ)方式

用一維數(shù)組a（1：p）和一維整型數(shù)組d（0：n）存放對(duì)稱矩陣A.

這里a（1：p）=（a11，a2m2，a22，a3m3，…，a33，…，aimi，…，aii，…，anmn，…，ann），其中p=∑ni=1（i-mi+1）是數(shù)組a中元素的個(gè)數(shù).另外用d（0）=0，d（i），i=1，2，…，n來(lái)標(biāo)記對(duì)角元素aii在數(shù)組a中的序號(hào).矩陣A的元素aij處于數(shù)組a中的第d（i）-i+j個(gè)位置上，即

aij=a（d（i）-i+j），

aimi中的列下標(biāo)mi的計(jì)算公式為

mi=i-（d（i）-d（i-1））+1，i=1，2，…，n.

矩陣B存放在二維數(shù)組b（1：n，1：m）的位置.最終解X也存放在b（1：n，1：m）的位置.

2.計(jì)算矩陣L，D和B～

L的元素存放在數(shù)組a（1：p）中A元素的相應(yīng)位置，L的對(duì)角線元素不存.D的元素di存放在a（1：p）中aii相應(yīng)的位置.B～存放在b（1：n，1：m）的位置.

對(duì)于i=1，2，…，n

（1）計(jì)算I=d（i）-i；MI←mi=i-（d（i）-d（i-1））+1.

（2）對(duì)于j=MI，MI+1，…，i

J=d（j）-j；MJ←mj=j-（d（j）-d（j-1））+1；

M←mij=max（MI，MJ）=max（mi，mj）

a（I+j）：=a（I+j）-∑j-1k=Ma（I+k）·a（d（k））·a（J+k）

如果j=i，則轉(zhuǎn)（3），否則a（I+j）：=a（I+j）/a（d（j））；

b（i，k）：=b（i，k）-a（I+j）·a（d（j））·b（j，k），k=1，2，…，m.

（3）如果a（I+i）=0，則計(jì)算停止（這時(shí)A非正定）

否則b（i，k）：=b（i，k）/a（d（i）），k=1，2，…，m.

3.求解方程組LTX=B～

對(duì)于i=n，n-1，…，1

（1）計(jì)算I=d（i）-i；

MI←mi=i-（d（i）-d（i-1））+1.

（2）對(duì)于j=MI，MI+1，…，i-1

b（j，k）：=b（j，k）-a（I+j）·b（i，k），k=1，2，…，m.