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

代码如下：

#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[n].row)&&(j==c[n].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");
}