稀疏矩阵的存储方法之二

转自：http://study.hhit.edu.cn/subject/CourseWare_Detail.aspx?TeachCourseWareID=257

1.基本概念

稀疏矩阵(SparseMatrix)：是矩阵中的一种特殊情况，其非零元素的个数远小于零元素的个数。
设m行n列的矩阵含t个非零元素，则称
 
以二维数组表示高阶的稀疏矩阵时，会产生零值元素占的空间很大且进行了很多和零值的运算的问题。 
特殊矩阵：值相同的元素或0元素在矩阵中的分布有一定的规律。如下三角阵、三对角阵、稀疏矩阵。
压缩存储：为多个值相同的元素只分配一个存储空间；对0元素不分配空间。目的是节省大量存储空间。
n x n的矩阵一般需要n2个存储单元，当为对称矩阵时需要n(1+n)/2个单元。

2.三元组顺序表——压缩存储稀疏矩阵方法之一（顺序存储结构）

三元组顺序表又称有序的双下标法，对矩阵中的每个非零元素用三个域分别表示其所在的行号、列号和元素值。它的特点是，非零元在表中按行序有序存储，因此便于进行依行顺序处理的矩阵运算。当矩阵中的非0元素少于1/3时即可节省存储空间。
（1） 稀疏矩阵的三元组顺序表存储表示方法
#define MAXSIZE 12500 // 假设非零元个数的最大值为12500
typedef struct {
int i, j; // 该非零元的行下标和列下标
ElemType e; //非零元素的值
} Triple; // 三元组类型
typedef union { //共用体
Triple data[MAXSIZE + 1]; // 非零元三元组表，data[0]未用
int mu, nu, tu; // 矩阵的行数、列数和非零元个数
} TSMatrix; // 稀疏矩阵类型
（2） 求转置矩阵的操作
◆用常规的二维数组表示时的算法
for (col=1; col<=nu; ++col)
for (row=1; row<=mu; ++row)
T[col][row] = M[row][col];
其时间复杂度为: O(mu×nu)
◆ 用三元组顺序表表示时的快速转置算法
Status FastTransposeSMatrix(TSMatrix M, TSMatrix &T) {
// 采用三元组顺序表存储表示，求稀疏矩阵M的转置矩阵T
T.mu = M.nu; T.nu = M.mu; T.tu = M.tu;
if (T.tu) {
for (col=1; col<=M.nu; ++col) num[col] = 0;
for (t=1; t<=M.tu; ++t) ++num[M.data[t].j];// 求 M 中每一列所含非零元的个数
cpot[1] = 1;
for (col=2; col<=M.nu; ++col) cpot[col] = cpot[col-1] + num[col-1];
// 求 M 中每一列的第一个非零元在 b.data 中的序号
for (p=1; p<=M.tu; ++p) { // 转置矩阵元素
col = M.data[p].j; q = cpot[col];
T.data[q].i =M.data[p].j; T.data[q].j =M.data[p].i;
T.data[q].e =M.data[p].e; ++cpot[col]; 
} // for
} // if
return OK;
} // FastTransposeSMatrix
其时间复杂度为: O(mu +nu) 

3.行逻辑联接的顺序表——压缩存储稀疏矩阵方法之二（链接存储结构）

行逻辑联接的顺序表：稀疏矩阵中为了随机存取任意一行的非0元素，需要知道每一行的第一个非0元素在三元组表中的位置，因此将上述快速转置算法中指示行信息的辅助数组cpot固定在稀疏矩阵的存储结构中，让每一行对应一个单链表，每个单链表都有一个表头指针，这种“带行链接信息”的三元组表即称为行逻辑联接的顺序表。
（1）行逻辑联接的顺序表的表示法
#define MAXMN 500 // 假设矩阵行数和列数的最大值为500
typedef struct {
Triple data[MAXSIZE + 1]; // 非零元三元组表，data[0]未用
int rpos[MAXMN + 1]; // 指示各行第一个非零元的位置
int mu, nu, tu; // 矩阵的行数、列数和非零元个数
} RLSMatrix; // 行逻辑链接顺序表类型
（2） 求矩阵乘法的操作
◆ 矩阵乘法的精典算法:
for (i=1; i<=m1; ++i)
for (j=1; j<=n2; ++j) {
Q[i][j] = 0;
for (k=1; k<=n1; ++k) Q[i][j] += M[i][k] * N[k][j];
}
其时间复杂度为：O(m1 n1 n2)
◆ 用行逻辑联接的顺序表表示时的矩阵乘法
Status MultSMatrix(RLSMatrix M, RLSMatrix N, RLSMatrix &Q) {
//求矩阵乘积Q=M*N，采用行逻辑链接存储表示。
if (M.nu != N.mu) return ERROR;
Q.mu = M.mu; Q.nu = N.nu; Q.tu = 0; // Q初始化
if (M.tu*N.tu != 0) { // Q是非零矩阵
for (arow=1; arow<=M.mu; ++arow) { // 处理M的每一行
ctemp[] = 0; // 当前行各元素累加器清零
Q.rpos[arow] = Q.tu+1; 
if (arow<M.mu) tp=M.rpos[arow+1];
else {tp=M.tu+1}
for (p=M.rpos[arow]; p<M.rpos[arow+1];++p) { 
//对当前行中每一个非零元找到对应元在N中的行号
brow=M.data[p].j; 
if (brow < N.nu ) t = N.rpos[brow+1];
else { t = N.tu+1 }
for (q=N.rpos[brow]; q< t; ++q) {
ccol = N.data[q].j; // 乘积元素在Q中列号
ctemp[ccol] += M.data[p].e * N.data[q].e;
} // for q
} // 求得Q中第crow( =arow)行的非零元
for (ccol=1; ccol<=Q.nu; ++ccol) // 压缩存储该行非零元
if (ctemp[ccol]) {
if (++Q.tu > MAXSIZE) return ERROR;
Q.data[Q.tu] = {arow, ccol, ctemp[ccol]};
} // if
} // for arow
} // if 
return OK;
} // MultSMatrix
上述算法的时间复杂度分析：
◆ 累加器ctemp初始化的时间复杂度为O (M.mu x N.mu)
◆求Q的所有非零元的时间复杂度为O (M.tu x N.tu/N.mu)
◆ 进行压缩存储的时间复杂度为O (M.mu x N.nu)
总的时间复杂度就是O (M.mu x N.nu + M.tu x N.tu/N.mu)。
若M是m行n列的稀疏矩阵，N是n行p列的稀疏矩阵，则M中非零元的个数 M.tu = d M x m x n，N中非零元的个数 N.tu = d N x n x p，相乘算法的时间复杂度就是 O (m x p x(1+nd Md N)) ，当d M<0.05 和d N<0.05及 n <1000时，相乘算法的时间复杂度就相当于 O (mxp)。
显然，这是一个相当理想的结果。如果事先能估算出所求乘积矩阵Q不再是稀疏矩阵，则以二维数组表示Q，相乘的算法也就更简单了。 

4. 十字链表——压缩存储稀疏矩阵方法之三（链接存储结构）

（1） 基本概念
十字链表：是既带行指针又带列指针的链接存储方式，每个三元组结点处于所在行单链表与列单链表的交点处，当矩阵的非零元个数和位置在操作过程中变化较大时，用这种存储结构更为恰当。
在十字链表中，每个非零元可用一个含五个域的结点表示，其中 i, j 和e 三个域分别表示该非零元所在的行、列和非零元的值，向右域 right 用以链接同一行中下一个非零元，向下域down 用以链接同一列中下一个非零元。同一行的非零元通过 right 域链接成一个线性链表，同一列的非零元通过 down 域链接成一个线性链表，每个非零元既是某个行链表中的一个结点，又是某个列链表中的一个结点，整个矩阵构成了一个十字交叉的链表，故称这样的存储结构为十字链表，可用两个分别存储行链表的头指针和列链表的头指针的一维数组表示之。例如：矩阵M的十字链表如下图所示。
 
　　假设非空指针 pa和 pb分别指向矩阵A和B中行值相同的两个结点，pa == NULL 表明矩阵A在该行中没有非零元，则上述四种情况的处理过程为：
(1) 若pa==NULL或pa->j 〉pb->j，则需要在A矩阵的链表中插入一个值为bi,j的结点。此时，需　　改变同一行中前一结点的right域值，以及同一列中前一结点的down域值。
(2) 若pa->j〈 pb->j，则只要将pa指针往右推进一步。
(3) 若pa->j == pb->j且pa->e＋pb->e !=0，则只要将ai,j＋bi,j 的值送到pa所指结点的e域即　　可，其它所有域的值都不变。
(4) 若pa->j == pb->j且pa->e＋pb->e == 0，则需要在A矩阵的链表中删除pa所指的结点。此时　　，需改变同一行中前一结点的right域值，以及同一列中前一结点的down域值。
　　为了便于插入和删除结点，还需要设立一些辅助指针。其一是，在A的行链表上设pre指针，指　　示pa所指结点的前驱结点；其二是，在A的每一列的链表上设一个指针hl[j]，它的初值和列链　　表的头指针相同，即hl[j]＝chead[j]。
（2） 稀疏矩阵的十字链表表示与建立的算法 （P：104）
（3） 两个矩阵相加的算法描述
(1) 初始令pa和pb分别指向A和B的第一行的第一个非零元素的结点，即
　　pa＝A.rhead[1]; pb＝B.rhead[1]; pre = NULL;
　　且令hl初始化 for (j=1; j<=A.nu; ++j) hl[j]＝A.chead[j];
(2) 重复本步骤，依次处理本行结点，直到B的本行中无非零元素的结点，即pb==NULL为止：
① 若pa==NULL或pa->j〉pb->j（即A的这一行中非零元素已处理完），则需在A中插入一个pb所指　　结点的复制结点。假设新结点的地址为p，则A的行表中的指针作如下变化：
if pre == NULL rhead[p->i]=p;
else { pre->right＝p; }
p->right＝pa; pre = p;

A的列链表中的指针也要作相应的改变。首先需从hl[p->j]开始找到新结点在同一列中的前驱结点，并让hl[p->j]指向它，然后在列链表中插入新结点：
if chead[p->j] == NULL 
{ chead[p->j] = p; p->down = NULL; }
else { 
p->down＝hl[p->j]->down; hl[p->j]->down＝p; }
hl[p->j] = p;
② 若pa->j〈pb->j且pa->j!=0，则令pa指向本行下一个非零元结点，即 pre＝pa; pa＝pa->right;
③ 若pa->j == pb->j，则将B中当前结点的值加到A中当前结点上，即
pa->e＋＝pb->e;

此时若pa->e!=0，则指针不变，否则删除A中该结点，即行表中指针变为：
if pre == NULL rhead[pa->i] = pa->right;
else { pre->right＝pa->right; }
p=pa; pa=pa->right; 

同时，为了改变列表中的指针，需要先找到同一列中的前驱结点，且让hl[pa->j]指向该结点，然后如下修改相应指针：
if chead[p->j] == p
chead[p->j] = hl[p->j] = p->down; 
else { hl[p->j]->down＝p->down; }
free (p);
(3) 若本行不是最后一行，则令pa和pb指向下一行的第一个非零元结点，转(2)；否则结束。
此算法时间复杂度：O(ta+tb)