稀疏矩阵算法

1、稀疏矩阵的压缩存储
　为了节省存储单元，可只存储非零元素。由于非零元素的分布一般是没有规律的，因此在存储非零元素的同时，还必须存储非零元素所在的行号、列号，才能迅速确定一个非零元素是矩阵中的哪一个元素。稀疏矩阵的压缩存储会失去随机存取功能。
　其中每一个非零元素所在的行号、列号和值组成一个三元组(i，j，aij)，并由此三元组惟一确定。
　稀疏矩阵进行压缩存储通常有两类方法：顺序存储和链式存储。链式存储方法。

2、三元组表
　将表示稀疏矩阵的非零元素的三元组按行优先(或列优先)的顺序排列(跳过零元素)，并依次存放在向量中，这种稀疏矩阵的顺序存储结构称为三元组表。
注意：
　 以下的讨论中，均假定三元组是按行优先顺序排列的。
　　　　　　
（1）三元组表的类型说明
　 为了运算方便，将矩阵的总行数、总列数及非零元素的总数均作为三元组表的属性进行描述。其类型描述为：
#define MaxSize 10000 //由用户定义
typedef int DataType； //由用户定义
typedef struct { //三元组
int i，j；//非零元的行、列号
DataType v； //非零元的值
}TriTupleNode；

typedef struct{ //三元组表
TriTupleNode data[MaxSize]； //三元组表空间
int m,n，t； //矩阵的行数、列数及非零元个数
}TriTupleTable；

(2) 压缩存储结构上矩阵的转置运算
一个m×n的矩阵A，它的转置矩阵B是一个n×m的矩阵，且：
A[i][j]=B[j][i]，0≤i<m，0≤j<n，
即A的行是B的列，A的列是B的行。

①三元组表表示的矩阵转置的思想方法
　　第一步：根据A矩阵的行数、列数和非零元总数确定B矩阵的列数、行数和非零元总数。
　　第二步：当三元组表非空（A矩阵的非零元不为0）时，根据A矩阵三元组表的结点空间data（以下简称为三元组表），将A的三元组表a->data置换为B的三元组表b->data。

②三元组表的转置
　方法一:简单地交换a->data中i和j中的内容，得到按列优先顺序存储倒b->data;再将b->data重排成按行优先顺序的三元组表。
　方法二:由于A的列是B的行，因此，按a->data的列序转置，所得到的转置矩阵B的三元组表b->data必定是按行优先存放的。
　按这种方法设计的算法，其基本思想是：对A中的每一列col(0≤col≤a->n-1)，通过从头至尾扫描三元组表a->data，找出所有列号等于col的那些三元组，将它们的行号和列号互换后依次放人b->data中，即可得到B的按行优先的压缩存贮表示。

下面是完整的稀疏矩阵算法，包括加法，转置，乘法等

#include <stdio.h>
#include <stdlib.h>
#define A 6
#define B 7
#define A1 3
#define B1 4
#define A2 4
#define B2 2
#define OK 1
#define ERROR -1;
#define MAXSIZE 1000

typedef int Status;
typedef int ElemType;
typedef struct {
int i,j;
ElemType e;
}Triple;
typedef struct {
Triple data[MAXSIZE+1];
int rpos[MAXSIZE+1];
int mu,nu,tu;
}TSMatrix;

void CreatTripleTable (int array_a[A][B],TSMatrix *a)
/*array_aÊÇÒ»¸öÏ¡Êè¾ØÕó,aÊDzúÉúµÄÏà¶ÔÓ¦µÄÈýÔª×é´æ´¢*/
{
int i,j,k=1;
//a->data=(Triple *)malloc((MAXSIZE+1)*sizeof(Triple));
for (i=0;i<A;i++) /*°´ÐÐÓÅÏÈ˳ÐòɨÃèarray_aµÄÔªËØ,²»Îª0Õß´æÈëBÖÐ*/
for (j=0;j<B;j++)
if (array_a[i][j]!=0)
{
a->data[k].i=i+1;
a->data[k].j=j+1;
a->data[k].e = array_a[i][j];
//printf("i=%d/n",a->data[k].i);
//printf("j=%d/n",a->data[k].j);
//printf("e=%d/n",a->data[k].e);
//printf("k=%d/n",k);
//printf("/n");
k++;
}
a->mu=A;
a->nu=B;
// data[0][1]=N;
a->tu=k-1;/*´æÈë·Ç0ÔªËظöÊý*/
}/*CreatTripleTable*/

void CreatTripleTable1 (int array_a[A1][B1],TSMatrix *a)
/*array_aÊÇÒ»¸öÏ¡Êè¾ØÕó,aÊDzúÉúµÄÏà¶ÔÓ¦µÄÈýÔª×é´æ´¢*/
{
int i,j,k=1;
//a->data=(Triple *)malloc((MAXSIZE+1)*sizeof(Triple));
for (i=0;i<A1;i++) /*°´ÐÐÓÅÏÈ˳ÐòɨÃèarray_aµÄÔªËØ,²»Îª0Õß´æÈëBÖÐ*/
for (j=0;j<B1;j++)
if (array_a[i][j]!=0)
{
a->data[k].i=i+1;
a->data[k].j=j+1;
a->data[k].e = array_a[i][j];
k++;
}
a->mu=A1;
a->nu=B1;
a->tu=k-1;/*´æÈë·Ç0ÔªËظöÊý*/
}

void CreatTripleTable2 (int array_a[A2][B2],TSMatrix *a)
/*array_aÊÇÒ»¸öÏ¡Êè¾ØÕó,aÊDzúÉúµÄÏà¶ÔÓ¦µÄÈýÔª×é´æ´¢*/
{
int i,j,k=1;
//a->data=(Triple *)malloc((MAXSIZE+1)*sizeof(Triple));
for (i=0;i<A2;i++) /*°´ÐÐÓÅÏÈ˳ÐòɨÃèarray_aµÄÔªËØ,²»Îª0Õß´æÈëBÖÐ*/
for (j=0;j<B2;j++)
if (array_a[i][j]!=0)
{
a->data[k].i=i+1;
a->data[k].j=j+1;
a->data[k].e = array_a[i][j];
k++;
}
a->mu=A2;
a->nu=B2;
a->tu=k-1;/*´æÈë·Ç0ÔªËظöÊý*/
}

Status TransposeSMatrix(TSMatrix *M,TSMatrix *T)
{
int p,q,col;
T->mu=M->nu; T->nu=M->mu; T->tu=M->tu;
if (T->tu) {
q=1;
for (col=1;col<=M->nu;++col)
for (p=1;p<=M->tu;++p)
if (M->data[p].j==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; ++q;
}
}
return OK;
}

Status FastTransposeSMatrix(TSMatrix *M,TSMatrix *T)
{
int t,p,q,col,num[A],cpot[M->nu];
T->mu=M->mu; 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];
cpot[1]=1;
for (col=2;col<=M->nu;++col) cpot[col]=cpot[col-1]+num[col-1];
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];
}
}
return OK;
}

Status TSMatrix_Add(TSMatrix M,TSMatrix N,TSMatrix *C)//ÈýÔª×é±íʾµÄÏ¡Êè¾ØÕó¼Ó·¨
{
int i,j,x,ce,pa,pb,pc;
if (M.mu!=N.mu || M.nu!=N.nu) return ERROR;
C->mu=M.mu;C->nu=M.nu;C->tu=0;
pa=1;pb=1;pc=1;
for(x=1;x<=M.mu;x++) //¶Ô¾ØÕóµÄÿһÐнøÐмӷ¨
{
while(M.data[pa].i<x) pa++;
while(M.data[pb].i<x) pb++;
while(M.data[pa].i==x&&N.data[pb].i==x)//ÐÐÁÐÖµ¶¼ÏàµÈµÄÔªËØ
{
if(M.data[pa].j==N.data[pb].j)
{
ce=M.data[pa].e+N.data[pb].e;
if(ce) //ºÍ²»Îª0
{
C->data[pc].i=x;
C->data[pc].j=M.data[pa].j;
C->data[pc].e=ce;
pa++;pb++;pc++;
}
}//if
else if(M.data[pa].j>N.data[pb].j)
{
C->data[pc].i=x;
C->data[pc].j=N.data[pb].j;
C->data[pc].e=N.data[pb].e;
pb++;pc++;
}
else
{
C->data[pc].i=x;
C->data[pc].j=M.data[pa].j;
C->data[pc].e=M.data[pa].e;
pa++;pc++;
}
}//while
while(M.data[pa].i==x) //²åÈëAÖÐÊ£ÓàµÄÔªËØ(µÚxÐÐ)
{
C->data[pc].i=x;
C->data[pc].j=M.data[pa].j;
C->data[pc].e=M.data[pa].e;
pa++;pc++;
}
while(N.data[pb].i==x) //²åÈëBÖÐÊ£ÓàµÄÔªËØ(µÚxÐÐ)
{
C->data[pc].i=x;
C->data[pc].j=N.data[pb].j;
C->data[pc].e=N.data[pb].e;
pb++;pc++;
}
}//for
C->tu=pc;
return OK;
}//TSMatrix_Add

void PrintSMatrix(TSMatrix M)
{
int i,j,pa;
int a[M.mu][M.nu];
pa=1;
for (i=0;i<M.mu;i++)
for (j=0;j<M.nu;j++)
a[i][j]=0;
for (pa=1;pa<=M.tu;pa++) {
//printf("%d/n",M.data[pa].i);
//printf("%d/n",M.data[pa].j);
//printf("%d/n",M.data[pa].e);
a[M.data[pa].i-1][M.data[pa].j-1]=M.data[pa].e;
}
for (i=0;i<M.mu;i++) {
for (j=0;j<M.nu;j++)
printf("%6d",a[i][j]);
printf("/n");
}
printf("/n");
}

void GetRpos(TSMatrix *M) {
int row,t,j,num[M->mu],cpot[M->nu];
for (row=1;row<=M->mu;++row) num[row]=0;
for (t=1;t<=M->tu;++t) ++num[M->data[t].i];
M->rpos[1]=1;
for (row=2;row<=M->mu;++row) M->rpos[row]=M->rpos[row-1]+num[row-1];
for (t=1;t<=M->mu;++t) printf("%d/n",M->rpos[t]);
}

Status MultSMatrix(TSMatrix M,TSMatrix N,TSMatrix *Q) {
int arrow,brow,tp,p,q,t,ccol,ctemp[Q->nu];
if (M.nu!=N.mu) return ERROR;
Q->mu=M.mu; Q->nu=N.nu; Q->tu=0;
if (M.tu*N.tu!=0) {
for (arrow=1;arrow<=M.mu;++arrow) {
for (p=1;p<=N.nu;p++) ctemp[p]=0;
Q->rpos[arrow]=Q->tu+1;
if (arrow<M.mu) tp=M.rpos[arrow+1];
else tp=M.tu+1;
for (p=M.rpos[arrow];p<tp;++p) {
brow=M.data[p].j;
if (brow<N.mu) t=N.rpos[brow+1];
else {t=N.tu+1;}
for (q=N.rpos[brow]; q<t; ++q) {
ccol=N.data[q].j;
//printf("ctemp[ccol]=%d/n",ctemp[ccol]);
ctemp[ccol]+=M.data[p].e * N.data[q].e;
//printf("ccol=%d/n",ccol);
//printf("M.data[p].e * N.data[q].e=%d * %d/n",M.data[p].e,N.data[q].e);
//printf("ctemp[ccol]=%d/n",ctemp[ccol]);
//printf("")
}
}
for (ccol=1; ccol<=Q->nu;++ccol)
if (ctemp[ccol]){
if (++Q->tu>MAXSIZE) return ERROR;
Q->data[Q->tu].i=arrow;
Q->data[Q->tu].j=ccol;
Q->data[Q->tu].e=ctemp[ccol];
}
}
}
}

void DestroySMatrix(TSMatrix M)
{
M.mu=0; M.nu=0; M.tu=0;
}

int main(int argc, char *argv[])
{
TSMatrix a,b,c,d,e; //Èç¹ûдΪ*a£¬Ôò±ØÐë¼ÓÉÏa=(TSMatrix *)malloc(sizeof(TSMatrix)),»ò¼ÓÒ»TSMatrix b,a=&b;
TSMatrix aa,bb,cc;
int array[A][B]={{0,12,9,0,0,0,0},{0,0,0,0,0,0,0},{-3,0,0,0,0,14,0},{0,0,24,0,0,0,0},{0,18,0,0,0,0,0},{15,0,0,-7,0,0,0}};
int array1[A][B]={{0,12,3,4,0,0,0},{0,0,0,3,0,0,0},{-3,0,0,0,0,14,0},{0,0,24,0,0,0,0},{0,18,0,0,0,0,0},{15,0,0,-7,0,0,0}};

int array2[A1][B1]={{3,0,0,5},{0,-1,0,0},{2,0,0,0}};
int array3[A2][B2]={{0,2},{1,0},{-2,4},{0,0}};

//a=(TSMatrix *)malloc(sizeof(TSMatrix));
CreatTripleTable(array,&a);
PrintSMatrix(a);
printf("Transpose/n");
TransposeSMatrix(&a,&c);
PrintSMatrix(c);
printf("Fast Transpose/n");
FastTransposeSMatrix(&a,&d);
PrintSMatrix(d);

printf("Add Transpose/n");
CreatTripleTable(array1,&b);
TSMatrix_Add(a,b,&e);
PrintSMatrix(e);

CreatTripleTable1(array2,&aa);
PrintSMatrix(aa);
CreatTripleTable2(array3,&bb);
PrintSMatrix(bb);
GetRpos(&aa);
GetRpos(&bb);
MultSMatrix(aa,bb,&cc);
PrintSMatrix(cc);

DestroySMatrix(a);
DestroySMatrix(b);
DestroySMatrix(c);
DestroySMatrix(aa);
DestroySMatrix(bb);
DestroySMatrix(cc);
system("PAUSE");
return 0;
}