5.3 稀疏矩阵
稀疏矩阵的三元组表存储
设m*n 矩阵中有t 个非零元素且t<<m*n，这样的矩阵称为稀疏矩阵。很多科学管理及工程计算中，常会遇到阶数很高的大型稀疏矩阵。如果按常规分配方法，顺序分配在计算机内，那将是相当浪费内存的。为此提出另外一种存储方法，仅仅存放非零元素。但对于这类矩阵，通常零元素分布没有规律，为了能找到相应的元素，所以仅存储非零元素的值是不够的，还要记下它所在的行和列。
于是采取如下方法：将非零元素所在的行、列以及它的值构成一个三元组（i,j,v），然后再按某种规律存储这些三元组，这种方法可以节约存储空间。下面讨论稀疏矩阵的压缩存储方法。
将三元组按行优先的顺序，同一行中列号从小到大的规律排列成一个线性表，称为三元组表，采用顺序存储方法存储该表。如图5.11 稀疏矩阵对应的三元组表为图5.12。
显然，要唯一的表示一个稀疏矩阵，还需要存储三元组表的同时存储该矩阵的行、列，为了运算方便，矩阵的非零元素的个数也同时存储。这种存储的思想实现如下：
1 define SMAX 1024 /*一个足够大的数*/ 2 typedef struct 3 { 4 int i,j; /*非零元素的行、列*/ 5 datatype v; /*非零元素值*/ 6 }SPNode; /*三元组类型*/ 7 typedef struct 8 { 9 int mu,nu,tu; /*矩阵的行、列及非零元素的个数*/ 10 SPNode data[SMAX]; /*三元组表*/ 11 } SPMatrix; /*三元组表的存储类型*/
这样的存储方法确实节约了存储空间，但矩阵的运算从算法上可能变的复杂些。
下面我们讨论这种存储方式下的稀疏矩阵的两种运算：转置和相乘。
1．稀疏矩阵的转置
设SPMatrix A; 表示一m*n 的稀疏矩阵，其转置B 则是一个n*m 的稀疏矩阵，因此也有SPMatrix B; 由A 求B 需要：
A 的行、列转化成B 的列、行；
将A.data 中每一三元组的行列交换后转化到B.data 中；
看上去以上两点完成之后，似乎完成了B，没有。因为我们前面规定三元组的是按一行一行且每行中的元素是按列号从小到大的规律顺序存放的，因此B 也必须按此规律实现，A 的转置B 如图5.13 所示，图5.14 是它对应的三元组存储，就是说，在A 的三元组存储基础上得到B 的三元组表存储（为了运算方便，矩阵的行列都从1 算起，三元组表data 也从1 单元用起）。
算法思路：
①A 的行、列转化成B 的列、行；
②在A.data 中依次找第一列的、第二列的、直到最后一列，并将找到的每个三元组的行、列交换后顺序存储到B.data 中即可。
算法如下：
1 void TransM1 (SPMatrix *A) 2 { 3 SPMatrix *B; 4 int p,q,col; 5 B=malloc(sizeof(SPMatrix)); /*申请存储空间*/ 6 B->mu=A->nu; 7 B->nu=A->mu; 8 B->tu=A->tu; 9 /*稀疏矩阵的行、列、元素个数*/ 10 if (B->tu>0) /*有非零元素则转换*/ 11 { 12 q=0; 13 for (col=1; col<=(A->nu); col++) /*按A 的列序转换*/ 14 for (p=1; p<= (A->tu); p++) /*扫描整个三元组表*/ 15 if (A->data[p].j==col ) 16 { 17 B->data[q].i= A->data[p].j ; 18 B->data[q].j= A->data[p].i ; 19 B->data[q].v= A->data[p].v; 20 q++; 21 }/*if*/ 22 } /*if(B->tu>0)*/ 23 return B; /*返回的是转置矩阵的指针*/ 24 } /*TransM1*/
算法5.1 稀疏矩阵转置
分析该算法，其时间主要耗费在col 和p 的二重循环上，所以时间复杂性为O(n*t)，(设m、n 是原矩阵的行、列，t 是稀疏矩阵的非零元素个数)，显然当非零元素的个数t 和m*n 同数量级时，算法的时间复杂度为O(m*n2)，和通常存储方式下矩阵转置算法相比，可能节约了一定量的存储空间，但算法的时间性能更差一些。
算法5.1 的效率低的原因是算法要从A 的三元组表中寻找第一列、第二列、…，要反复搜索A 表，若能直接确定A 中每一三元组在B 中的位置，则对A 的三元组表扫描一次即可。这是可以做到的，因为A 中第一列的第一个非零元素一定存储在B.data[1]，如果还知道第一列的非零元素的个数，那么第二列的第一个非零元素在B.data 中的位置便等于第一列的第一个非零元素在B.data 中的位置加上第一列的非零元素的个数，如此类推，因为A 中三元组的存放顺序是先行后列，对同一行来说，必定先遇到列号小的元素，这样只需扫描一遍A.data 即可。
根据这个想法，需引入两个向量来实现：num[n+1]和cpot[n+1]，num[col]表示矩阵A中第col 列的非零元素的个数（为了方便均从1 单元用起），cpot［col］初始值表示矩阵A中的第col 列的第一个非零元素在B.data 中的位置。于是cpot 的初始值为：
1 cpot[1]=1; 2 cpot[col]=cpot[col-1]+num[col-1]; 2≤col≤n
例如对于矩阵图5.11 矩阵A 的num 和cpot 的值如下：
依次扫描A.data，当扫描到一个col 列元素时，直接将其存放在B.data 的cpot[col]位置上，cpot[col]加１，cpot[col]中始终是下一个col 列元素在B.data 中的位置。下面按以上思路改进转置算法如下：
1 SPMatrix * TransM2 (SPMatrix *A) 2 { 3 SPMatrix *B; 4 int i,j,k; 5 int num[n+1],cpot[n+1]; 6 B=malloc(sizeof(SPMatrix)); /*申请存储空间*/ 7 B->mu=A->nu; 8 B->nu=A->mu; 9 B->tu=A->tu; 10 /*稀疏矩阵的行、列、元素个数*/ 11 if (B->tu>0) /*有非零元素则转换*/ 12 { 13 for (i=1;i<=A->nu;i++) 14 num[i]=0; 15 for (i=1;i<=A->tu;i++) /*求矩阵A 中每一列非零元素的个数*/ 16 { 17 j= A->data[i].j; 18 num[j]++; 19 } 20 cpot[1]=1; /*求矩阵A 中每一列第一个非零元素在B.data 中的位置*/ 21 for (i=2;i<=A->nu;i++) 22 cpot[i]= cpot[i-1]+num[i-1]; 23 for (i=1; i<= (A->tu); i++) /*扫描三元组表*/ 24 { 25 j=A->data[i].j; /*当前三元组的列号*/ 26 k=cpot[j]; /*当前三元组在B.data 中的位置*/ 27 B->data[k].i= A->data[i].j ; 28 B->data[k].j= A->data[i].i ; 29 B->data[k].v= A->data[i].v; 30 cpot[j]++; 31 } /*for i */ 32 } /*if (B->tu>0)*/ 33 return B; /*返回的是转置矩阵的指针*/ 34 } /*TransM2*/
算法5.2 稀疏矩阵转置的改进算法
分析这个算法的时间复杂度：这个算法中有四个循环，分别执行n，t，n-1，t 次，在每个循环中，每次迭代的时间是一常量，因此总的计算量是O（n+t）。当然它所需要的存储空间比前一个算法多了两个向量。
2．稀疏矩阵的乘积
已知稀疏矩阵A(m1× n1)和B(m2× n2)，求乘积C(m1× n2)。稀疏矩阵A、B、C 及它们对应的三元组表A.data、B.data、C.data 如图5.16 所示。
由矩阵乘法规则知：
这就是说只有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(列号与行号相等)且均不为零（三元组存在）时才相乘，并且累加到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 所示。
根据以上分析，稀疏矩阵的乘法运算的粗略步骤如下：
⑴初始化。清理一些单元，准备按行顺序存放乘积矩阵;
⑵求B 的num，rpot;
⑶做矩阵乘法。将A.data 中三元组的列值与B.data 中三元组的行值相等的非零元素相乘，并将具有相同下标的乘积元素相加。
算法如下：
1 SPMatrix *MulSMatrix (SPMatrix *A, SPMatrix *B) 2 /*稀疏矩阵A(m1× n1)和B(m2× n2) 用三元组表存储，求A×B */ 3 { 4 SPMatrix *C; /* 乘积矩阵的指针*/ 5 int p,q,i,j,k,r; 6 datatype temp[n+1]; 7 int num[B->nu+1],rpot[B->nu+1]; 8 if (A->nu!=B->mu) /*A 的列与B 的行不相等*/ 9 return NULL; 10 C=malloc(sizeof(SPMatrix)); /*申请C 矩阵的存储空间*/ 11 C->mu=A->mu; C->nu=B->nu; 12 if (A->tu*B->tu==0) 13 { 14 C->tu=0; 15 return C; 16 } 17 for (i=1;i<= B->mu;i++) /*求矩阵B 中每一行非零元素的个数*/ 18 num[i]=0; 19 for (k=1;k<=B->tu;k++) 20 { 21 i= B->data[k].i; 22 num[i]++; 23 } 24 rpot[1]=1; /*求矩阵B 中每一行第一个非零元素在B.data 中的位置*/ 25 for (i=2;i<=B->mu;i++) 26 rpot[i]= rpot[i-1]+num[i-1]; 27 r=0; /*当前C 中非零元素的个数*/ 28 p=1; /*指示A.data 中当前非零元素的位置*/ 29 for ( i= 1;i<=A->mu; i++) 30 { 31 for (j=1;j<=B->nu;j++) 32 temp[j]=0; /*cij 的累加器初始化*/ 33 while (A->data[p].i==i ) ./*求第i 行的*/ 34 { 35 k=A->data[p].j; /*A 中当前非零元的列号*/ 36 if (k<B->mu) 37 t=rpot[k+1]; 38 else 39 t=B->tu+1; /*确定B 中第k 行的非零元素在B.data 中的下限位置*/ 40 for (q=rpot[k]; q<t; q++;) /* B 中第k 行的每一个非零元素*/ 41 { 42 j=B->data[q].j; 43 temp[j]+=A->data[p].v * B->data[q].v 44 } 45 p++; 46 } /* while */ 47 for (j=1;j<=B->nu;j++) 48 if (temp[j] ) 49 { 50 r++; 51 C->data[r]={i,j,temp[j] }; 52 } 53 } /*for i*/ 54 C->tu=r; 55 return C; 56 }/* MulSMatrix */
算法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 是一个稀疏矩阵的十字链表。
用十字链表表示稀疏矩阵的基本思想是：对每个非零元素存储为一个结点，结点由5个域组成，其结构如图5.19 表示，其中：row 域存储非零元素的行号，col 域存储非零元素的列号，v 域存储本元素的值，right，down 是两个指针域。
稀疏矩阵中每一行的非零元素结点按其列号从小到大顺序由right 域链成一个带表头结点的循环行链表，同样每一列中的非零元素按其行号从小到大顺序由down 域也链成一个带表头结点的循环列链表。
即每个非零元素aij 既是第i 行循环链表中的一个结点，又是第j 列循环链表中的一个结点。行链表、列链表的头结点的row 域和col 域置0。每一列链表的表头结点的down 域指向该列链表的第一个元素结点，每一行链表的表头结点的right域指向该行表的第一个元素结点。
由于各行、列链表头结点的row 域、col 域和v 域均为零，行链表头结点只用right 指针域，列链表头结点只用right 指针域，故这两组表头结点可以合用，也就是说对于第i 行的链表和第i 列的链表可以共用同一个头结点。
为了方便地找到每一行或每一列，将每行（列）的这些头结点们链接起来，因为头结点的值域空闲，所以用头结点的值域作为连接各头结点的链域，即第i 行（列）的头结点的值域指向第i+1行（列）的头结点，… ，形成一个循环表。
这个循环表又有一个头结点，这就是最后的总头结点，指针HA 指向它。总头结点的row 和col 域存储原矩阵的行数和列数。
因为非零元素结点的值域是datatype 类型，在表头结点中需要一个指针类型，为了使整个结构的结点一致，我们规定表头结点和其它结点有同样的结构，因此该域用一个联合来表示；改进后的结点结构如图5.20 所示。
综上，结点的结构定义如下:
1 typedef struct node 2 { 3 int row, col; 4 struct node *down , *right; 5 union v_next 6 { 7 datatype v; 8 struct node *next; 9 } 10 } MNode，*MLink;
让我们看基于这种存储结构的稀疏矩阵的运算。这里将介绍两个算法，创建一个稀疏矩阵的十字链表和用十字链表表示的两个稀疏矩阵的相加。
1．建立稀疏矩阵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 列)链表的头结点。
这样做可以在建立链表时随机的访问任何一行（列），为建表带来方便。
算法如下：
1 MLink CreatMLink( ) /* 返回十字链表的头指针*/ 2 { 3 MLink H; 4 Mnode *p,*q,*hd[s+1]; 5 int i,j,m,n,t; 6 datatype v; 7 scanf(“%d,%,%d”,&m,&n,&t); 8 H=malloc(sizeof(MNode)); /*申请总头结点*/ 9 H->row=m; H->col=n; 10 hd[0]=H; 11 for(i=1; i<S; i++) 12 { 13 p=malloc(sizeof(MNode)); /*申请第i 个头结点*/ 14 p->row=0; 15 p->col=0; 16 p->right=p; 17 p->down=p; 18 hd[i]=p; 19 hd[i-1]->v_next.next=p; 20 } 21 hd[S]->v_next.next=H; /*将头结点们形成循环链表*/ 22 for (k=1;k<=t;k++) 23 { 24 scanf (“%d,%d,%d”,&i,&j,&v); /*输入一个三元组，设值为int*/ 25 p=malloc(sizeof(MNode)); 26 p->row=i ; 27 p->col=j; 28 p->v_next.v=v 29 /*以下是将*p 插入到第i 行链表中去，且按列号有序*/ 30 q=hd[i]; 31 while ( q->right!=hd[i] && (q->right->col)<j ) /*按列号找位置*/ 32 q=q->right; 33 p->right=q->right; /*插入*/ 34 q->right=p; 35 /*以下是将*p 插入到第j 行链表中去，且按行号有序*/ 36 q=hd[i]; 37 while ( q->down!=hd[j] && (q->down->row)<i ) /*按行号找位置*/ 38 q=q->down; 39 p-> down =q-> down; /*插入*/ 40 q-> down =p; 41 } /*for k*/ 42 return H; 43 } /* CreatMLink */
算法5.4 建立稀疏矩阵的十字链表
上述算法中，建立头结点循环链表时间复杂度为O(S)，插入每个结点到相应的行表和列表的时间复杂度是O(t*S)，这是因为每个结点插入时都要在链表中寻找插入位置，所以总的时间复杂度为O(t*S)。
该算法对三元组的输入顺序没有要求。如果我们输入三元组时是按以行为主序（或列）输入的，则每次将新结点插入到链表的尾部的，改进算法后，时间复杂度为O(S+t)。
２．两个十字链表表示的稀疏矩阵的加法
已知两个稀疏矩阵A 和B，分别采用十字链表存储，计算C=A+B，C 也采用十字链表方式存储，并且在A 的基础上形成C。
由矩阵的加法规则知，只有A 和B 行列对应相等，二者才能相加。
C 中的非零元素cij 只可能有３种情况：或者是aij+bij，或者是aij (bij=0)，或者是bij (aij=0)，因此当B 加到A 上时，对A 十字链表的当前结点来说，对应下列四种情况：或者改变结点的值（aij+bij≠０），或者不变（bij＝０），或者插入一个新结点（aij＝０），还可能是删除一个结点（aij+bij＝０）。
整个运算从矩阵的第一行起逐行进行。对每一行都从行表的头结点出发，分别找到A 和B 在该行中的第一个非零元素结点后开始比较，然后按４种不同情况分别处理。设pa和pb 分别指向A 和B 的十字链表中行号相同的两个结点，４种情况如下：
(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 所指结点。
由前面建立十字链表算法知，总表头结点的行列域存放的是矩阵的行和列，而各行（列）链表的头结点其行列域值为零，当然各非零元素结点的行列域其值不会为零，下面分析的4 种情况利用了这些信息来判断是否为表头结点。
综上所述，算法如下：
1 MLink AddMat (Ha,Hb) 2 { 3 MLink Ha,Hb; 4 Mnode *p,*q,*pa,*pb,*ca,*cb,*qa; 5 if (Ha->row!=Hb->row || Ha->col!=Hb->col) 6 return NULL; 7 ca=Ha->v_next.next; /*ca 初始指向A 矩阵中第一行表头结点*/ 8 cb=Hb->v_next.next; /*cb 初始指向B 矩阵中第一行表头结点*/ 9 do 10 { 11 pa=ca->right; /*pa 指向A 矩阵当前行中第一个结点*/ 12 qa=ca; /*qa 是pa 的前驱*/ 13 pb=cb->right; /*pb 指向B 矩阵当前行中第一个结点*/ 14 while (pb->col!=0) /*当前行没有处理完*/ 15 { 16 if (pa->col < pb->col && pa->col !=0 ) /*第三种情况*/ 17 { 18 qa=pa; 19 pa=pa->right; 20 } 21 else if (pa->col > pb->col || pa->col ==0 ) /*第四种情况*/ 22 { 23 p=malloc(sizeof(MNode)); 24 p->row=pb->row; 25 p->col=pb->col; 26 p->v=pb->v; 27 p->right=pa;/* 新结点插入*pa 的前面*/ 28 qa->right=p; 29 pa=p; 30 /*新结点还要插到列链表的合适位置，先找位置，再插入*/ 31 q=Find_JH(Ha,p->col); /*从列链表的头结点找起*/ 32 while(q->down->row!=0 && q->down->row<p->row) 33 q=q->down; 34 p->down=q->down; /*插在*q 的后面*/ 35 q->down=p; 36 pb=pb->right; 37 } /* if */ 38 else /*第一、二种情况*/ 39 { 40 x= pa->v_next.v+ pb->v_next.v; 41 if (x==0) /*第二种情况*/ 42 { 43 qa->right=pa->right; ./*从行链中删除*/ 44 /*还要从列链中删除，找*pa 的列前驱结点*/ 45 q= Find_JH (Ha,pa->col); /*从列链表的头结点找起*/ 46 while ( q->down->row < pa->row ) 47 q=q->down; 48 q->down=pa->down; 49 free (pa); 50 pa=qa; 51 } /*if (x==0)*/ 52 else /*第一种情况*/ 53 { 54 pa->v_next.v=x; 55 qa=pa; 56 } 57 pa=pa->right; 58 pb=pb->right; 59 } 60 } /*while*/ 61 ca=ca->v_next.next; /*ca 指向A 中下一行的表头结点*/ 62 cb=cb->v_next.next; /*cb 指向B 中下一行的表头结点*/ 63 } while (ca->row==0) /*当还有未处理完的行则继续*/ 64 return Ha; 65 }
算法5.5 十字链表表示的稀疏矩阵相加
为了保持算法的层次，在上面的算法，用到了一个函数findjH。
函数Mlink Find_JH(MLink H, int j)的功能是：返回十字链表H 中第j 列链表的头结点指针，很简单，读者可自行写出。
posted on 2015-05-21 11:35 chunlanse2014 阅读(1318) 评论(0) 编辑 收藏 举报
