5.3 矩阵的压缩存储
视频二维码（扫码观看）
在科学与工程计算问题中，矩阵是一种常用的数学对象，在高级语言编程时，通常将一个矩阵描述为一个二维数组。这样，可以对其元素进行随机存取，各种矩阵运算也非常简单。
压缩存储：
①多个相同的非零元素只分配一个存储空间；
②零元素不分配空间。
一、特殊矩阵
特殊矩阵是指值相同的元素或者0元素有规律分布的矩阵。
1对称矩阵
若一个n阶方阵A＝（aij）n×n中的元素满足性质：aij＝aji，其中1≤i，j≤n则称A为对称矩阵。
对称矩阵中的元素关于主对角线对称，因此，给每一对对称元素aij和aji（i≠j）分配一个存储空间，则n2个元素压缩存储到n（n＋1）/2个存储空间，能节约近一半的存储空间。
不失一般性，假设按“行优先顺序”存储下三角形（包括对角线）中的元素。
设用一维数组（向量）sa[n（n＋1）/2]存储n阶对称矩阵。为了便于访问，矩阵元素aij和向量sa[k]之间的对应关系为
对于矩阵中的任意元素aij，均可在一维数组sa中唯一确定其位置k；反之，对所有k＝1，2，…，n（n＋1）/2，都能确定sa[k]中的元素在矩阵中的位置（i，j）。则称sa[n（n＋1）/2]为n阶对称矩阵A的压缩存储。
图5-3 对称矩阵的压缩存储
2三角矩阵
对称矩阵的压缩存储方式同样适用于三角矩阵。
上三角矩阵的下三角（不包括主对角线）中的元素均为常数c（一般为0）。只需在对称矩阵存储基础上，增加常数c的存储空间。
3对角矩阵
矩阵中，除了主对角线和主对角线上或下方若干条对角线上的元素之外，其余元素皆为零。即所有的非零元素集中在以主对角线为了中心的带状区域中。
对角矩阵可按行优先顺序或对角线顺序，将其压缩存储到一个向量中，并且也能找到每个非零元素和向量下标的对应关系。
上述各种特殊矩阵，其非零元素的分布都是有规律的，因此总能找到一种方法将它们压缩存储到一个向量中，并且一般都能找到矩阵中的元素与该向量的对应关系，通过这个关系，仍能对矩阵的元素进行随机存取。
二、稀疏矩阵
【定义】对于稀疏矩阵，目前还没有一个确切的定义。设矩阵A是一个n×m的矩阵中有s个非零元素，设δ＝s/（n×m），称δ为稀疏因子，如果某一矩阵的稀疏因子δ满足δ≤0.05时称为稀疏矩阵。
对于稀疏矩阵，采用压缩存储方法时，只存储非0元素。必须存储非0元素的行下标值、列下标值、元素值。因此，一个三元组（i，j，aij）唯一确定稀疏矩阵的一个非零元素。
稀疏矩阵：
其三元组线性表为：（（1，2，12），（1，3，9），（3，1，－3），（3，8，4），（4，3，24），（5，2，18），（6，7，－7），（7，4，－6））。
1三元组顺序表
若以行序为主序，稀疏矩阵中所有非0元素的三元组，就可以构成该稀疏矩阵的一个三元组顺序表。相应的数据结构定义如下：
（1）三元组结点定义
#define MAX_SIZE 101
typedef struct
{
int row;/*行下标*/
int col;/*列下标*/
ElemType value;/*元素值 */
}Triple;
（2）三元组顺序表定义
typedef struct
{
int rn;/*行数*/
int cn;/*列数*/
int tn;/*非0元素个数*/
Triple data[MAX_SIZE+1];//非零元三元组表，data[0]未用
}TSMatrix;
矩阵的运算包括矩阵的转置、矩阵求逆、矩阵的加减、矩阵的乘除等。在此，先讨论在这种压缩存储结构下的求矩阵的转置的运算。
一个m×n的矩阵A，它的转置B是一个n×m的矩阵，且b[i][j]＝a[j][i]，0≤i≤n，0≤j≤m，即B的行是A的列，B的列是A的行。
设稀疏矩阵A是按行优先顺序压缩存储在三元组表a.data中，若仅仅是简单地交换a.data中i和j的内容，得到三元组表b.data，b.data将是一个按列优先顺序存储的稀疏矩阵B，要得到按行优先顺序存储的b.data，就必须重新排列三元组表b.data中元素的顺序。
图5-4 稀疏矩阵A分别按行、列优先顺序压缩
求转置矩阵的基本算法思想是：
①将矩阵的行、列下标值交换。即将三元组表中的行、列位置值i、j相互交换；
②重排三元组表中元素的顺序。即交换后仍然是按行优先顺序排序的。
方法一：
算法思想：按稀疏矩阵A的三元组表a.data中的列次序依次找到相应的三元组存入b.data中。
每找转置后矩阵的一个三元组，需从头至尾扫描整个三元组表a.data。找到之后自然就成为按行优先的转置矩阵的压缩存储表示。
算法描述：
void TransMatrix(TSMatrix a,TSMatrix b)
{
int p,q,col;
b.rn=a.cn;
b.cn=a.rn;
b.tn=a.tn;/*置三元组表b.data的行、列数和非0元素个数*/
if(b.tn==0)
printf("The Matrix A=0\n");
else
{
q=0;
for(col=1;col<=a.cn;col++)/*每循环一次找到转置后的一个三元组*/
for(p=0;p<a.tn;p++)/*循环次数是非0元素个数*/
if(a.data[p].col==col)
{
b.data[q].row=a.data[p].col;
b.data[q].col=a.data[p].row;
b.data[q].value=a.data[p].value;
q++;
}
}
}
算法分析：本算法主要的工作是在p和col的两个循环中完成的，故算法的时间复杂度为O（cn×tn），即矩阵的列数和非0元素的个数的乘积成正比。
而一般传统矩阵的转置算法为：
for(col=1;col<=n;++col)
for(row=0;row<=m;++row)
b[col][row]=a[row][col];
其时间复杂度为O（n×m）。当非零元素的个数tn和m×n同数量级时，算法TransMatrix的时间复杂度为O（m×n2）。
由此可见，方法一虽然节省了存储空间，但时间复杂度却大大增加。所以上述算法只适合于稀疏矩阵中非0元素的个数tn远远小于m×n的情况。
视频二维码（扫码观看）
方法二（快速转置的算法）
算法思想：直接按照稀疏矩阵A的三元组表a.data的次序依次顺序转换，并将转换后的三元组放置于三元组表b.data的恰当位置。
前提：若能预先确定A中每一列的（即B中每一行）第一个非0元素在b.data中应有的位置，则在作转置时就可直接放在b.data中恰当的位置。因此，应先求得A中每一列的非0元素个数。
附设两个辅助数组num[ ]和cpot[ ]。
num[col]：统计A中第col列中非0元素的个数；
cpot[col]：指示A中第col列第一个非0元素在b.data中的恰当位置。
显然有位置对应关系：
cpot[1]＝1
cpot[col]＝cpot[col－1]＋num[col－1]，其中2≤col≤a.cn
矩阵A的num[col]和cpot[col]的值，如表5-1所示。
表5-1 num[col]和cpot[col]的值表
算法描述：
void FastTransMatrix(TMatrix a,TMatrix b)
{
int p,q,col,k;
int num[MAX_SIZE],copt[MAX_SIZE];
b.rn=a.cn;
b.cn=a.rn;
b.tn=a.tn;/*置三元组表b.data的行、列数和非0元素个数*/
if(b.tn==0)
printf("The Matrix A=0\n");
else
{
for(col=1;col<=a.cn;++col)
num[col]=0;/*向量num[]初始化为0*/
for(k=1;k<=a.tn;++k)
++num[a.data[k].col];/*求原矩阵中每一列非0元素个数*/
for(cpot[0]=1,col=2;col<=a.cn;++col)
cpot[col]=cpot[col-1]+num[col-1];/*求第col列中第一个非0元在b.data中的序号*/
for(p=1;p<=a.tn;++p)
{
col=a.data[p].col;
q=cpot[col];
b.data[q].row=a.data[p].col;
b.data[q].col=a.data[p].row;
b.data[q].value=a.data[p].value;
++cpot[col];/*至关重要，本列中其他非0元素 */
}
}
}
2行逻辑链接的三元组顺序表
将上述方法二中的辅助数组cpot[ ]固定在稀疏矩阵的三元组表中，用来指示“行”的信息。得到另一种顺序存储结构：行逻辑链接的三元组顺序表。其类型描述如下：
#define MAX_ROW 100
typedef struct
{
Triple data[MAX_SIZE+1];/*非0元素的三元组表*/
int rpos[MAX_ROW+1];/*各行第一个非0位置表*/
int rn,cn,tn;/*矩阵的行、列数和非0元个数*/
}RLSMatrix;
稀疏矩阵的乘法
设有两个矩阵：A＝（aij）m×n，B＝（bij）n×p，则C＝（cij）m×p，其中cij＝∑aik×bkj，1≤k≤n，1≤i≤m，1≤j≤p
经典算法是三重循环：
for(i=1;i<=m;++i)
for(j=1;j<=p;++j)
{
c[i][j]=0;
for(k=1;k<=n;++k)
c[i][j]=c[i][j]+a[i][k]×b[k][j];
}
此算法的复杂度为O（m×n×p）。
设有两个稀疏矩阵A＝（aij）m×n，B＝（bij）n×p，其存储结构采用行逻辑链接的三元组顺序表。
算法思想：对于A中的每个元素a.data[p]（p＝1，2，…，a.tn），找到B中所有满足条件：
a.data[p].col＝b.data[q].row的元素b.data[q]，求得a.data[p].value×b.data[q].value，该乘积是cij中的一部分。求得所有这样的乘积并累加求和就能得到cij。
为得到非0的乘积，只要对a.data[1…a.tn]中每个元素（i，k，aik）（1≤i≤a.rn，1≤k≤a.cn），找到b.data中所有相应的元素（k，j，bkj）（1≤k≤b.rn，1≤j≤b.cn）相乘即可。则必须知道矩阵B中第k行的所有非0元素，而b.rpos[ ]向量中提供了相应的信息。
b.rpos[row]指示了矩阵B的第row行中第一个非0元素在b.data[ ]中的位置（序号），显然，b.rpos[row＋1]－1指示了第row行中最后一个非0元素在b.data[ ]中的位置（序号）。最后一行中最后一个非0元素在b.data[ ]中的位置显然就是b.tn。
算法描述：
void MultsMatrix(RLSMatrix a,RLSMatrix b,RLSMatrix c)
{/*求矩阵A、B的积C＝A×B，采用行逻辑链接的顺序表*/
ElemType ctemp[Max_Size];
int p,q,arow,ccol,brow,t;
if(a.cn!=b.rn)
{
printf("Error\n");
exit(0);
}
else
{
c.rn=a.rn;
c.cn=b.n;
c.tn=0;/*初始化C*/
if(a.tn*b.tn!=0)
{/*C是非零矩阵*/
for(arow=1;arow<=a.rn;++arow)
{
ctemp[arow]=0;/*当前行累加器清零*/
c.rpos[arow]=c.tn+1;
p=a.rops[arow];
for(;p<a.rpos[arow+1];++p)
{/*对第arow行的每一个非0元素*/
brow=a.data[p].col;/*找到元素在b.data[]中的行号*/
if(brow<b.cn)
t=b.rpos[brow+1];
else
t=b.tn+1;
for(q=b.rpos[brow];q<t;++q)
{
ccol=b.data[q].col;/*乘积元素在c中的列号*/
ctemp[ccol]+=a.data[p].value*b.data[q].value;
}
}/*求出c中第arow行中的非0元素*/
for(ccol=1;ccol<=c.cn;++ccol)
{
if(ctemp[ccol]!=0)
{
if(++c.tn>MAX_SIZE)
{
printf("Error\n");
exit(0);
}
else
c.data[c.tn]=(arow,ccol,ctemp[ccol]);
}
}
}
}
}
}
视频二维码（扫码观看）
3十字链表
对于稀疏矩阵，当非0元素的个数和位置在操作过程中变化较大时，采用链式存储结构表示，比三元组的线性表更方便。
稀疏矩阵的十字链表存储表示：
typedef struct Clnode
{
int row,col;/*行号和列号*/
ElemType value;/*元素值*/
struct CLNode *down,*right;//该非0元素所在的行表和列表的后继链域
}OLNode;/*非0元素结点*/
typedef struct
{
int rn;/*矩阵的行数*/
int cn;/*矩阵的列数*/
int tn;/*非0元素总数*/
OLink *chead,*rhead;//行、列表头指针向量，基址由CreateSMatrix分配
}CrossList;
由定义知，稀疏矩阵中同一行的非0元素的由right指针域链接成一个行链表，由down指针域链接成一个列链表。则每个非0元素既是某个行链表中的一个结点，同时又是某个列链表中的一个结点，所有的非0元素构成一个十字交叉的链表。称为十字链表。
此外，还可用两个一维数组分别存储行链表的头指针和列链表的头指针。
对于稀疏矩阵M：
对应的十字交叉链表如图5-5所示。
图5-5 稀疏矩阵M的十字交叉链表