稀疏矩阵的转置操作及其乘法操作

#include <stdio.h>
#include <stdlib.h>

#define MAX_COL 20
#define MAX_TERMS 50
#define EXIT_FALLURE -1

struct term{                 //稀疏矩阵的单个元素
int row;				//该元素所在行
int col;                //该元素所在列
int value;				//该元素的值
};

typedef struct term term;
term a[4]={{3,3,3},{0,0,1},{1,0,1},{2,0,1}};            //矩阵的初始化
term b[4]={{3,3,3},{0,0,1},{0,1,1},{0,2,2}};
term d[MAX_TERMS];
term newB[MAX_TERMS];

void output(term c[])
{
int i,j,n=1;
for(i=0; i<c[0].row; i++){
for(j=0; j<c[0].col; j++)
{
if((i==c
.row)&&(j==c
.col)){
printf("%d  ",c[n++].value);
}
else printf("0  ");
}
printf("\n");
}
}

int compare(int a,int b){
if(a<b) return -1;
else if(a>b) return 1;
else return 0;
}

void transpose(term a[],term b[])    //稀疏矩阵的转置
{//b is set to the transpose of a
int n,i,j,currentb;
n=a[0].value; //total number of elements
b[0].row = a[0].col; //rows in b = cols in a
b[0].col = a[0].row; //cols in b = rows in a
b[0].col = n;
if(n>0) { //non zero matrix
currentb = 1;
for (i=0; i<a[0].col; i++)
//transpose by the cols in a
for (j=1;j<=n;j++)
//transpose by the cols in a
if (a[j].col==i) {
//element is in current col,add it to b
b[currentb].row = a[j].col;
b[currentb].col = a[j].row;
b[currentb].value = a[j].value;
currentb++;
}
}
}

/*void fastTranspose(term a[],term b[])    //稀疏矩阵的快速转置
{//the transpose of a is placed in b
int rowTerms[MAX_COL],startingPos{MAX_COL];
int i,j,numCols = a[0].col, numTerms = a[0].value;
b[0].row = numCols;
b[0].col = a[0].row;
b[0].value = numTerms;
if (numTerms>0) {//nonzero matrix
for (i=0; i<numCols; i++)
rowTerms[i] = 0;
for (i=1; i<=numTerms; i++)
rowTerms[a[i].col]++;
startingPos[0] =1;
for (i=1; i<=numCols; i++)
startingPos[i] = startingPos[i-1] + rowTerms[i-1];
for (i=1;i<=numTerms; i++){
j = startingPos[a[i].col]++;
b[j].row = a[i].col;
b[j].col = a[i].row;
b[j].value = a[i].value;
}
}
}*/

void storeSum(term d[], int *totalD, int row, int col, int *sum)
{	/*if *sum!=0,then it along with its row and col.position is stored as the *totalD+1 entry in d */
if (*sum)
if (*totalD<MAX_TERMS) {
d[++*totalD].row = row;
d[*totalD].col = col;
d[*totalD].value = *sum;
*sum = 0;
}else {
//fprintf(stderr,"Numbers of terms in productexceeds %d\n",MAX_TERMS);
exit(EXIT_FAILURE);
}
}

void mmult(term a[], term b[], term d[])               //矩阵的乘法
{/*multiply two sparse matrices*/
int i, j, col, totalB=b[0].value, totalD=0;
int rowsA=a[0].row, colsA=a[0].col;
int colsB=b[0].col, totalA=a[0].value;
int rowBegin=1, row=a[1].row, sum=0;

/*if (colsA!=b[0].row)
{
fprintf(stderr,"Incompatible matrices\n");
exit(EXIT_FALLURE);
}*/

transpose(b, newB);
/*set boundary condition*/
a[totalA+1].row = rowsA;
newB[totalB+1].row = colsB; newB[totalB+1].col = 0;
for (i=1; i<=totalA; ) {
col = newB[1].row;
for (j=1; j<=totalB+1; ) {
/*multiply row of a by col of b */
if (a[i].row!=row) {//如果a的这一行取完了，就直接开始从newB的下一行的第一个开始取，并把i回复，因为这一行的任务还没完成
storeSum(d, &totalD, row, col, &sum);//把值存入目标矩阵
i = rowBegin;
for (; newB[j].row==col; j++) {;}
col = newB[j].row;
}else if (newB[j].row!=col) {//如果b的这一行取完了
storeSum(d, &totalD, row, col, &sum);//把值存入目标矩阵
i = rowBegin; col = newB[j].row;
}else{
switch (compare(a[i].col, newB[j].col)) {
case -1://go to next term in a
i++; break;
case 0://add terms,go to next term in a and b
sum += (a[i++].value*newB[j++].value);
break;
case 1://advance to next term in b
j++;
}
}
}
for (; a[i].row==row; i++) {;}    //开始取a的下一行
rowBegin = i; row = a[i].row;
}
d[0].row = rowsA; d[0].col = colsB; d[0].value = totalD;
}

void main()
{
output(a);
printf("\n");
output(b);
printf("\n");
mmult(a,b,d);
output(d);
system("pause");
}