§5.3 矩陣的壓縮存儲
矩陣是很多科學與工程計算問題中研究的數學對象，在此，我們討論如何存儲矩陣的元，從而使矩陣的各種運算能有效第進行。對于一個矩陣結構顯然用一個二維數組來表示是非常恰當的，但在有些情況下，比如常見的一些特殊矩陣，如三角矩陣、對稱矩陣、帶狀矩陣、稀疏矩陣等，從節約存儲空間的角度考慮，對這類矩陣進行壓縮存儲，即：為多個值相同的元素或只分配一個存儲單元，對零元素不分配空間。下面從這一角度來考慮這些特殊矩陣的存儲方法。
值相同的元素或者零元素在矩陣中的分布有一定的規律，則稱此類矩陣為特殊矩陣；反之為稀疏矩陣。

5.3.1 特殊矩陣
對稱矩陣:
對稱矩陣的特點是：在一個n階方陣中，有aij=aji ，其中1≤i , j≤n，如圖5.5所示是一個5階對稱矩陣。對稱矩陣關于主對角線對稱，因此只需存儲上三角或下三角部分即可，比如，我們只存儲下三角中的元素aij，其特點是中j≤i且1≤i≤n，對于上三角中的元素aij ，它和對應的aji相等，因此當訪問的元素在上三角時，直接去訪問和它對應的下三角元素即可，這樣，原來需要n*n個存儲單元，現在只需要n(n+1)/2個存儲單元了，節約了n(n-1)/2個存儲單元，當n較大時，這是可觀的一部分存儲資源。


如何只存儲下三角部分呢？對下三角部分以行為主序順序存儲到一個向量中去，在下三角中共有n*(n+1)/2個元素，因此，不失一般性，設存儲到向量SA[n(n+1)/2]中，存儲順序可用圖5.6示意，這樣，原矩陣下三角中的某一個元素aij則具體對應一個sak，下面的問題是要找到k與i、j之間的關系。

對于下三角中的元素aij，其特點是：i≥j且1≤i≤n，存儲到SA中后，根據存儲原則，它前面有i-1行，共有1+2+…+i-1=i*(i-1)/2個元素，而aij又是它所在的行中的第j個，所以在上面的排列順序中，aij是第i*(i-1)/2+j個元素，因此它在SA中的下標k與i、j的關系為： k=i*(i-1)/2+j-1 (0≤k<n*(n+1)/2 )

若i<j，則aij是上三角中的元素，因為aij=aji ，這樣，訪問上三角中的元素aij時則去訪問和它對應的下三角中的aji即可，因此將上式中的行列下標交換就是上三角中的元素在SA中的對應關系：
k=j*(j-1)/2+i-1 (0≤k<n*(n+1)/2 )

綜上所述，對于對稱矩陣中的任意元素aij，若令I=max(i,j)，J=min(i,j)，則將上面兩個式子綜合起來得到：　k=I*(I-1)/2+J-1。稱SA[n(n+1)/2]為n階對稱矩陣A的壓縮存儲。

三角矩陣
形如圖5.7的矩陣稱為三角矩陣，其中c為某個常數。其中5.7(a)為下三角矩陣：主隊角線以上均為同一個常數；(b)為上三角矩陣，主隊角線以下均為同一個常數；下面討論它們的壓縮存儲方法。

下三角矩陣
與對稱矩陣類似，不同之處在于存完下三角中的元素之后，緊接著存儲對角線上方的常量，因為是同一個常數，所以存一個即可，這樣一共存儲了n*(n+1)+1個元素，設存入向量：SA[n*(n+1)+1]中，這種的存儲方式可節約n*(n-1)-1個存儲單元，sak 與aji 的對應關系為：


2. 上三角矩陣
對于上三角矩陣，存儲思想與下三角類似，以行為主序順序存儲上三角部分，最后存儲對角線下方的常量。對于第1行，存儲n個元素，第2行存儲n-1個元素，…，第p行存儲(n-p+1)個元素，aij的前面有i-1行，共存儲：
n+(n-1) +…+(n-i+1)= =(i-1)*(2n-i+2)/2個元素，而aij 是它所在的行中要存儲的第（j-i+1）個；所以，它是上三角存儲順序中的第 (i-1)*(2n-i+2)/2+(j-i+1)個，因此它在SA中的下標為：k=(i-1)*(2n-i+2)/2+j-i。
綜上， sak 與aji 的對應關系為：

對角矩陣
在這種矩陣中，所有的非零元都集中在以主對角線為中心的帶狀區域中，即除了主對角線上和直接在對角線上方、下方若干條對角線上的元外，所有的其它元皆為零。對于這種矩陣，我們也可按某個原則（以行為主，或以對角線為主的順序）將其壓縮存儲道一維數組中。見P96 圖5.4

5.3.2 稀疏矩陣
一、定義：
設m*n矩陣中有t個非零元素且t<<m*n，這樣的矩陣稱為稀疏矩陣。令,稱為矩陣的稀疏因子，通常認為≤0.05時稱為稀疏矩陣。很多科學管理及工程計算中，常會遇到階數很高的大型稀疏矩陣。

二、抽象數據類型稀疏矩陣的定義
ADT SparseMatrix{
數據對象：
數據關系：
基本操作:
CreateSMatrix(&M)
創建稀疏矩陣
DestroySMatrix(&M)
銷毀稀疏矩陣
PrintSMatrix(&M)
輸出稀疏矩陣
CopySMatrix(M，&T)
復制稀疏矩陣
AddSMatrix(M，N，&Q)
求稀疏矩陣的和
SubSMatrix(&M)
求稀疏矩陣的差
MultSMatrix(&M)
求稀疏矩陣的積
TransposeSMatrix(&M)
求稀疏矩陣的轉置矩陣
} ADT SparseMatrix

三、稀疏矩陣的三元組表存儲
如果按常規分配方法，順序分配在計算機內，那將是相當浪費內存的。為此提出另外一種存儲方法，僅僅存放非零元素。但對于這類矩陣，通常零元素分布沒有規律，為了能找到相應的元素，所以僅存儲非零元素的值是不夠的，還要記下它所在的行和列。于是采取如下方法：將非零元素所在的行、列以及它的值構成一個三元組（i,j,v），然后再按某種規律存儲這些三元組，這種方法可以節約存儲空間。下面討論稀疏矩陣的壓縮存儲方法。
將三元組按行優先的順序，同一行中列號從小到大的規律排列成一個線性表，稱為三元組表（有序的雙下標法），采用順序存儲方法存儲該表。如圖5.11稀疏矩陣對應的三元組表為圖 5.12。

顯然，要唯一的表示一個稀疏矩陣，還需要存儲三元組表的同時存儲該矩陣的行、列，為了運算方便，矩陣的非零元素的個數也同時存儲。這種存儲的思想實現如下
define MAXSIZE 12500 /*一個足夠大的數*/
typedef struct
{ int i,j; /*非零元素的行、列*/
elemtype e; /*非零元素值*/
}Triple; /*三元組類型*/
typedef struct
{ int mu,nu,tu; /*矩陣的行、列及非零元素的個數*/
Triple data[MAXSIZE]; /*三元組表*/
} TSMatrix; /*三元組表的存儲類型*/
這樣的存儲方法確實節約了存儲空間，但矩陣的運算從算法上可能變的復雜些。下面我們討論這種存儲方式下的稀疏矩陣的兩種運算：轉置和相乘。
稀疏矩陣的轉置
設TSMatrix A; 表示一m*n的稀疏矩陣，其轉置B則是一個n*m的稀疏矩陣，因此也有 TSMatrix B; 由A求B需要：
A的行、列轉化成B的列、行；
將A.data中每一三元組的行列交換后轉化到B.data中；
看上去以上兩點完成之后，似乎完成了B，沒有。因為我們前面規定三元組的是按一行一行且每行中的元素是按列號從小到大的規律順序存放的，因此B也必須按此規律實現，A的轉置B如圖5.13所示，圖5.14是它對應的三元組存儲，就是說，在A的三元組存儲基礎上得到B的三元組表存儲（為了運算方便，矩陣的行列都從1算起，三元組表data也從1單元用起）。


轉置實現算法思路：
①將矩陣的行列值相互交換；
②將每個三元組中的I和j相互調換；
③重排三元組之間的次序便可實現矩陣的轉置。
第三步的實現：
A、按照b.data中三元組的次序依次在a.data中找到相應的三元組進行轉置。具體算法描述見P98 算法5.1。其僅適于tu《mu×nu的情況。
B、按照a.data中的三元組的次序進行轉置，并將轉置后的三元組置入b中恰當的位置。為求出每列中b.data的恰當位置，在轉置前，應先求得M的每一列中非零元的個數。
需要附設num和cpot兩個向量。Num[col]表示矩陣M中第col列中非零元的個數，cpot[col]指示M中第col列的第一個非零元在b.data中的恰當位置。顯然有：copt[1]=1
copt[col]= copt[col-1]+ Num[col-1] 2≤col≤a.nu
算法見P100 算法5.2
將這種轉置稱為快速轉置。

總結：非零元在表中按行序有序存儲，因此便于進行以行順序處理的矩陣運算，然而若序按行號存取某一行的非零元，則需從頭開始進行查找。

2．稀疏矩陣的乘積
為了實現矩陣相乘，應便于隨機存取任意一行的非零元，則需知道每一行的第一個非零元在三元組表中的位置，為此，將快速轉置矩陣算法中創建的，指示“行”信息的輔助數組copt固定在稀疏矩陣的存儲結構中，稱這種“帶行鏈接信息”的三元組表稱為邏輯連接的順序表，其類型描述如下：
Typedef struct{
Triple data[MAXSIZE+1];
Int rpos[MAXRC+1];
Int mu,nu,tu;
}RLSMatrix;
下面討論兩個矩陣相乘已知稀疏矩陣A(m1× n1)和B(m2× n2)，求乘積C(m1× n2)。
稀疏矩陣A、B、C及它們對應的三元組表A.data、B.data、C.data如圖5.16所示。

由矩陣乘法規則知：
C（i，j）=A(i,1)×B(1,j)+A(i,2)×B(2,j)+…+A(i,n)×B(n,j)
=
這就是說只有A(i,k)與B(k,p)（即A元素的列與B元素的行相等的兩項）才有相乘的機會，且當兩項都不為零時，乘積中的這一項才不為零。
矩陣用二維數組表示時，傳統的矩陣乘法是：A的第一行與B的第一列對應相乘累加后得到c11，A的第一行再與B的第二列對應相乘累加后得到c12，…，因為現在按三元組表存儲，三元組表是按行為主序存儲的，在B.data中，同一行的非零元素其三元組是相鄰存放的，同一列的非零元素其三元組并未相鄰存放，因此在B.data中反復搜索某一列的元素是很費時的，因此改變一下求值的順序，以求c11和c12為例，因為

即a11只有可能和B中第1行的非零元素相乘，a12只有可能和B中第2行的非零元素相乘，…，而同一行的非零元是相鄰存放的，所以求c11和c12同時進行：求a11*b11累加到c11，求a11*b12累加到c12，再求a12*b21累加到c11，再求a12*b22累加到c22.，…，當然只有aik和bkj(列號與行號相等)且均不為零（三元組c11= c12= 解釋
a11*b11+ a11*b12+ a11只與B中1行元素相乘
a12*b21+ a12*b22+ a12只與B中2行元素相乘
a13*b31+ a13*b32+ a13只與B中3行元素相乘
a14*b41 a14*b42 a14只與B中4行元素相乘
存在）時才相乘，并且累加到cij當中去。
為了運算方便，設一個累加器：datatype temp[n+1]；用來存放當前行中cij的值，當前行中所有元素全部算出之后，再存放到C.data中去。
為了便于B.data中尋找B中的第k行第一個非零元素，與前面類似，在此需引入num和rpot兩個向量。num[k]表示矩陣B中第k行的非零元素的個數；rpot[k]表示第k行的第一個非零元素在B.data中的位置。
于是有
rpot[1]=1
rpot[k]=rpot[k-1]+num[k-1] 2≤k≤n
例如 ，對于矩陣B的num和rpot如圖5.17所示。

col 1 2 3 4
num[col] 2 0 2 0
cpot[col] 1 3 3 5
圖5.17 矩陣B的num與rpot值
根據以上分析，稀疏矩陣的乘法運算的粗略步驟如下：
⑴初始化。清理一些單元，準備按行順序存放乘積矩陣;
⑵求B的num，rpot;
⑶做矩陣乘法。將A.data中三元組的列值與B.data中三元組的行值相等的非零元素相乘，并將具有相同下標的乘積元素相加。見P102
算法5.3
分析上述算法的時間性能如下：（1）求num的時間復雜度為O(B->nu+B->tu)；（2）求rpot 時間復雜度為O(B->mu)；（3）求temp時間復雜度為O(A->mu*B->nu)；（4）求C的所有非零元素的時間復雜度為O(A->tu*B->tu/B->mu)；（5）壓縮存儲時間復雜度為O(A->mu*B->nu)；所以總的時間復雜度為O(A->mu*B->nu+(A->tu*B->tu)/B->nu)。
四、稀疏矩陣的十字鏈表存儲
三元組表可以看作稀疏矩陣順序存儲，但是在做一些操作（如加法、乘法）時，非零項數目及非零元素的位置會發生變化，這時這種表示就十分不便。在這節中，我們介紹稀疏矩陣的一種鏈式存儲結構——十字鏈表，它同樣具備鏈式存儲的特點，因此，在某些情況下，采用十字鏈表表示稀疏矩陣是很方便的。
圖5.18是一個稀疏矩陣的十字鏈表。

1、▲：用十字鏈表表示稀疏矩陣的基本思想是：對每個非零元素存儲為一個結點，結點由5個域組成，其結構如圖5.19表示，其中：row域存儲非零元素的行號，col域存儲非零元素的列號，
域存儲本元素的非零元值，right，down是兩個指針域。
圖5.19 十字鏈表的結點結構

圖5.19 十字鏈表的結點結構
稀疏矩陣中每一行的非零元素結點按其列號從小到大順序由right域鏈成一個帶表頭結點的循環行鏈表，同樣每一列中的非零元素按其行號從小到大順序由down域也鏈成一個帶表頭結點的循環列鏈表。即每個非零元素aij既是第i行循環鏈表中的一個結點，又是第j列循環鏈表中的一個結點。整個矩陣構成一個十字交叉的鏈表，顧稱這樣的存儲結構為十字鏈表。可用兩個存儲行鏈表的頭指針和列鏈表的頭指針的一維數組之。見P104 圖5.6

2、稀疏矩陣的十字鏈表存儲表示：見P104 算法5.4

3、矩陣運算：
①．建立稀疏矩陣A的十字鏈表
首先輸入的信息是：m（A的行數），n（A的列數），r（非零項的數目），緊跟著輸入的是r個形如（i,j,aij）的三元組。
算法的設計思想是：首先建立每行（每列）只有頭結點的空鏈表，并建立起這些頭結點拉成的循環鏈表；然后每輸入一個三元組（i，j，aij），則將其結點按其列號的大小插入到第i個行鏈表中去，同時也按其行號的大小將該結點插入到第j個列鏈表中去。在算法中將利用一個輔助數組MNode *hd[s+1]; 其中 s=max(m , n) , hd [i]指向第i行(第i列)鏈表的頭結點。這樣做可以在建立鏈表時隨機的訪問任何一行（列），為建表帶來方便。
見P 104 算法5.4 建立稀疏矩陣的十字鏈表
上述算法中，建立頭結點循環鏈表時間復雜度為O(S)，插入每個結點到相應的行表和列表的時間復雜度是O(t*S)，這是因為每個結點插入時都要在鏈表中尋找插入位置，所以總的時間復雜度為O(t*S)。該算法對三元組的輸入順序沒有要求。如果我們輸入三元組時是按以行為主序（或列）輸入的，則每次將新結點插入到鏈表的尾部的，改進算法后，時間復雜度為O(S+t)。
②．兩個十字鏈表表示的稀疏矩陣的加法 v已知兩個稀疏矩陣A和B，分別采用十字鏈表存儲，計算C=A+B，C也采用十字鏈表方式存儲，并且在A的基礎上形成C。
由矩陣的加法規則知，只有A和B行列對應相等，二者才能相加。C中的非零元素cij只可能有3種情況：或者是aij+bij，或者是aij　(bij=0)，或者是bij　(aij=0)，因此當B加到A上時，對A十字鏈表的當前結點來說，對應下列四種情況：或者改變結點的值（aij+bij≠0），或者不變（bij＝0），或者插入一個新結點（aij＝0），還可能是刪除一個結點（aij+bij＝0）。整個運算從矩陣的第一行起逐行進行。對每一行都從行表的頭結點出發，分別找到A和B在該行中的第一個非零元素結點后開始比較，然后按4種不同情況分別處理。

設pa和pb分別指向A和B的十字鏈表中行號相同的兩個結點，4種情況如下：
(1) 若pa->col=pb->col且pa->v+pb->v≠0，則只要用aij+bij的值改寫pa所指結點的值域即可。
(2) 若pa->col=pb->col且pa->v+pb->v=0，則需要在矩陣A的十字鏈表中刪除pa所指結點，此時需改變該行鏈表中前趨結點的right域，以及該列鏈表中前趨結點的down域。
(3) 若pa->col < pb->col且pa->col≠0（即不是表頭結點），則只需要將pa指針向右推進一步，并繼續進行比較。
(4) 若pa->col > pb->col或pa->col＝0（即是表頭結點），則需要在矩陣A的十字鏈表中插入一個pb所指結點。

為了便于運算，在A的行鏈表上設pre指針，指示pa所指結點的前驅結點；在A的每一列的鏈表上設一個指針hl[j]，它的初值和列鏈表的頭指針相同，即hl[j]=chead[j]。矩陣相加操作的概要描述，見P105。