稀疏矩阵的实现

什么是稀疏矩阵，怎么定义才是稀疏矩阵?假设在m x n的矩阵中，有t个不为0的元素。令z=t/(m x n),若z<=0.05则称此矩阵为稀疏矩阵。由于稀疏矩阵的非0元素较少，所以如果用大容量存储必定会造成内存浪费，因此只需存储非0元素值即可，以下列出了常用的三种存储稀疏矩阵的方法。

三元组顺序表

三元组顺序表的实现用顺序存储结构来实现。定义一个数组的结构体，存储的是三元组（即由 3 部分数据组成的集合），组中数据分别表示（行标，列标，元素值）。结构体的实现如下：

//三元组结构体
typedef struct {
    int i,j;//行标i，列标j
    int data;//元素值
}triple;
//矩阵的结构表示
typedef struct {
    triple data[number];//存储该矩阵中所有非0元素的三元组
    int n,m,num;//n和m分别记录矩阵的行数和列数，num记录矩阵中所有的非0元素的个数
}TSMatrix;

如上图所示即为依次再数组中插入的 三元表。具体实现过程如下所示：

#include <stdio.h>
#define MAXSIZE 100
#define OK 1
#define ERROR 0
#define ELEMTYPE int

typedef struct {
	int i;   //行
	int j;   //列
	ELEMTYPE e; //元素
}Triple;
typedef struct {
	Triple data[MAXSIZE+1];
	int mu, nu, tu;
}TSmatrix;
void display(TSmatrix *s);
int main(void)
{
	TSmatrix s;
	s.mu = 3;
	s.nu = 3;
	s.tu = 3;
	s.data[0].i = 3;
	s.data[0].j = 3;
	s.data[0].e = 6;
	s.data[1].i = 2;
	s.data[1].j = 3;
	s.data[1].e = 8;
	s.data[2].i = 2;
	s.data[2].j = 1;
	s.data[2].e = 4;
	display(&s);
	getchar();
	return 0;
}

void display(TSmatrix *s)
{
	for (int i = 1; i <= s->mu; i++)
	{
		for (int j = 1; j <= s->nu; j++)
		{
			int value = 0;            //用来判断是否到了指定的（i,j）位置
			for (int k = 0; k < s->tu; k++)//遍历数组中的三元表值
			{
				if (s->data[k].i == i && s->data[k].j == j) //若遍历至指定（i,j）就打印元素值
				{
					value = 1;
					printf("%d ", s->data[k].e);
					break;
				}
			}
			if(value==0)  //若不为三元表中存储的值就打印0
			printf("%d ", 0);
		
		}
		printf("\n");
	}
}

结果如上图所示，检查矩阵可知打印的位置是没有问题的。
如上所示的方法能够实现稀疏矩阵的存储，但是在数据提取方面效率比较低，现在只是3x3的矩阵可能看不出来，但若是m和n都以千万计的情况下稀疏矩阵中的非0元素个数近似于n，那么时间复杂度就会达到O（m x n xn），这将非常浪费时间。因此在此方法上改进有了下面的方法。

行逻辑链接的顺序表

行逻辑链接顺序表示基于上种方法的改进，在上面的基础上再增加一个数组用来存储用来用于存储三元表的数组中每一行第一个非0元素的位置，这样就无需遍历整个三元表数组，只需要遍历对应行的数据就可以了，大大提高了效率。


使用数组 rpos 记录矩阵中每行第一个非 0 元素在一维数组中的存储位置。遍历的具体过程如下所示：
由 rpos 数组可知，第一行首个非 0 元素位于data[1]，因此在遍历此行时，可以直接从第 data[1] 的位置开始，一直遍历到下一行首个非 0 元素所在的位置（data[3]之前；
同样遍历第二行时，由 rpos 数组可知，此行首个非 0 元素位于 data[3]，因此可以直接从第 data[3] 开始，一直遍历到下一行首个非 0 元素所在的位置（data[4]）之前；
遍历第三行时，由 rpos 数组可知，此行首个非 0 元素位于 data[4]，由于这是矩阵的最后一行，因此一直遍历到 rpos 数组结束即可（也就是 data[tu]，tu 指的是矩阵非 0 元素的总个数）。

#include <stdio.h>
#define MAXSIZE 100
#define MAXSD  100
#define ELEMTYPE int

typedef struct {
	int i;
	int j;
	ELEMTYPE e;
}Triple;
typedef struct {
	Triple data[MAXSIZE];
	int rpos[MAXSD];  //定义rpos数组用于存储每一行第一个非0元素
	int mu, nu, tu;
}RLSmatrix;
void display(RLSmatrix *s);
int main(void)
{
	RLSmatrix s;
	s.mu = 4;
	s.nu = 3;
	s.tu = 4;

	s.data[1].i = 1;
	s.data[1].j = 1;
	s.data[1].e = 3;
	s.data[2].i = 2;
	s.data[2].j = 2;
	s.data[2].e = 4;
	s.data[3].i = 3;
	s.data[3].j = 3;
	s.data[3].e = 6;
	s.data[4].i = 4;
	s.data[4].j = 3;
	s.data[4].e = 3;

	s.rpos[1] = 1;
	s.rpos[2] = 2;
	s.rpos[3] = 3;
	s.rpos[4] = 4;
	display(&s);
	getchar();
	return 0;

}

void display(RLSmatrix *s)
{
	for (int i = 1; i <= s->mu; i++)
	{
		for (int j = 1; j <= s->nu; j++)
		{
			int value = 0;
			if(i+1<=s->mu) //假设一共i行则前面的i-1行遍历对应行的元素
			{
				for (int k = s->rpos[i]; k < s->rpos[i + 1]; k++) //遍历数组中对应行的非0元素即可
				{
					if (s->data[k].i == i && s->data[k].j == j) 
					{
						value = 1;
						printf("%d ", s->data[k].e);
						break;
					}
					if (value == 0)
						printf("0 ");
				}
			}
			else  //最后一行只需遍历到最后一个非0元素即可
			{
				for (int k = s->rpos[i]; k <=s->tu; k++)
				{
					if (s->data[k].i == i && s->data[k].j == j)
					{
						value = 1;
						printf("%d ", s->data[k].e);
						break;
					}
					if (value == 0)
						printf("0 ");
				}
			}
			
		}
		printf("\n");
	}
}

十字链表法存储稀疏矩阵

上面两种存储方式用于固定存储时可以很好的起作用，但是当要涉及非0元素的插入或者删除的时候回引起元素值的移动，例如两矩阵相加的操作，这种时候用链式存储表示三元组更为恰当，该存储方式采用的是 “链表+数组” 结构。


由上图可以看到，非0元素用一个含有五个域的节点表示，两个指针与分别用来指向同一列和同一行的元素。再用两个存储行链表和列链表的一维数组存储这些链表。每一个非0元既是某行链表的一个节点，又是列链表的一个节点。整个矩阵构成了一个十字交叉的链表。

#include<stdio.h>
#include<stdlib.h>
typedef struct OLNode
{
    int i, j, e; //矩阵三元组i代表行 j代表列 e代表当前位置的数据
    struct OLNode *right, *down; //指针域 右指针 下指针
}OLNode, *OLink;
typedef struct
{
    OLink *rhead, *chead; //行和列链表头指针
    int mu, nu, tu;  //矩阵的行数,列数和非零元的个数
}CrossList;
CrossList CreateMatrix_OL(CrossList M);
void display(CrossList M);
int main()
{
    CrossList M;
    M.rhead = NULL;
    M.chead = NULL;
    M = CreateMatrix_OL(M);
    printf("输出矩阵M:\n");
    display(M);
    return 0;
}
CrossList CreateMatrix_OL(CrossList M)
{
    int m, n, t;
    int i, j, e;
    OLNode *p, *q;
    printf("输入矩阵的行数、列数和非0元素个数：");
    scanf("%d%d%d", &m, &n, &t);
    M.mu = m;
    M.nu = n;
    M.tu = t;
    if (!(M.rhead = (OLink*)malloc((m + 1) * sizeof(OLink))) || !(M.chead = (OLink*)malloc((n + 1) * sizeof(OLink))))
    {
        printf("初始化矩阵失败");
        exit(0);
    }
    for (i = 1; i <= m; i++)
    {
        M.rhead[i] = NULL;
    }
    for (j = 1; j <= n; j++)
    {
        M.chead[j] = NULL;
    }
    for (scanf("%d%d%d", &i, &j, &e); 0 != i; scanf("%d%d%d", &i, &j, &e)) {
        if (!(p = (OLNode*)malloc(sizeof(OLNode))))
        {
            printf("初始化三元组失败");
            exit(0);
        }
        p->i = i;
        p->j = j;
        p->e = e;
        //链接到行的指定位置
        if (NULL == M.rhead[i] || M.rhead[i]->j > j)
        {
            p->right = M.rhead[i];
            M.rhead[i] = p;
        }
        else
        {
            for (q = M.rhead[i]; (q->right) && q->right->j < j; q = q->right);
            p->right = q->right;
            q->right = p;
        }
        //链接到列的指定位置
        if (NULL == M.chead[j] || M.chead[j]->i > i)
        {
            p->down = M.chead[j];
            M.chead[j] = p;
        }
        else
        {
            for (q = M.chead[j]; (q->down) && q->down->i < i; q = q->down);
            p->down = q->down;
            q->down = p;
        }
    }
    return M;
}
void display(CrossList M) {
    for (int i = 1; i <= M.nu; i++)
    {
        if (NULL != M.chead[i])
        {
            OLink p = M.chead[i];
            while (NULL != p)
            {
                printf("%d\t%d\t%d\n", p->i, p->j, p->e);
                p = p->down;
            }
        }
    }
}

【数据结构】稀疏矩阵乘法

1.传统矩阵相乘的算法使用三个嵌套循环实现，算法复杂度为O（m * n1 * n2）
2.使用三元组顺序表存储稀疏矩阵时，实现 Q= M * N，对于M中M（i，j）元素来说，只需要与N中第j行元素N（j,q）相乘，再存入Q（i，q）中。为了实现这一操作，增加一个向量rpos，表示每一行的第一个非零元在三元组中的位置，rpos作用相当于快速转置中的cpot向量。
这种结构叫做 行链接的顺序表

typedef struct {
	Triple data[MAXSIZE + 1]; //非零三元组表，data0未用
	int rpos[MAXRC + 1];	  //各行第一个非零元的位置表
	int mu, nu, tu;			  //三元组行数，列数，非零元素值
}RLSMatrix;

在创建矩阵时，需要求得矩阵的rpos向量，求法与cpot一致。

Status CreateRLSMatrix(RLSMatrix &M) {
	cout << "输入矩阵行数，列数，非零元素个数：" << endl;
	cin >> M.mu >> M.nu >> M.tu;
	cout << "依次输入" << M.tu << "个元素的行数，列数，元素值:" << endl;
	for (int i = 1; i <= M.tu; i++)
	{
		cin >> M.data[i].i >> M.data[i].j >> M.data[i].e;
	}
	int* num = (int*)malloc(M.mu * sizeof(int));
	int arow;
	int q;
	if (M.tu) {
		for ( arow = 1; arow <= M.mu; ++arow) num[arow] = 0;
		for (int t = 1; t <= M.tu; ++t) ++num[M.data[t].i];//求M中每一行含非零元个数
		M.rpos[1] = 1;
		//求第col列中第一个非零元在b.data 中的序号
		for (arow = 2; arow <= M.mu; ++arow)M.rpos[arow] = M.rpos[arow - 1] + num[arow - 1];
	}
		cout << "矩阵构造完成" << endl;
	return OK;
}

稀疏矩阵乘法
算法思路：
大循环是对M中的每一行进行。
rpos向量用来确定每行元素个数
再对该行每个元素M（i，j）进行处理，找到N中对应的第 j 行，将这一行所有元素与M（i，j）相乘，结果存入ctemp累加器中，该行元素都处理完成后，ctemp中所存储的便是Q中第 i 行数据

Status MultSMatrix(RLSMatrix M, RLSMatrix N, RLSMatrix &Q) {
	//求稀疏矩阵乘积 Q = M x N  采用行逻辑链接存储表示
	if (M.nu != N.mu)return ERROR;
	Q.mu = M.mu; Q.nu = N.nu; Q.tu = 0;
	int arow,tp,p,brow,t,q,ccol;
	int* ctemp = (int*)malloc((M.mu+1) * sizeof(int));
	if (M.tu * N.tu != 0) {
		for (arow = 1; arow <= M.mu; ++arow) {
		//对M的每行进行处理
			for (int i = 0; i <= M.mu; i++)ctemp[i] = 0;//累加器清零
			Q.rpos[arow] = Q.tu + 1;
			if (arow < M.mu)tp = M.rpos[arow + 1];
			else tp = M.tu + 1;//确定M中该行元素个数
			for ( p = M.rpos[arow]; p < tp; p++)
			{//M中对该行的每个元素进行处理
				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) {
				//将该元素与矩阵N中对应行的元素进行相乘，结果存入累加器
					ccol = N.data[q].j;
					ctemp[ccol] += M.data[p].e *N.data[q].e;
				}
			}
			for(ccol =1;ccol<=Q.nu;++ccol)
				if (ctemp[ccol]) {
				//第arow行处理完毕，将累加器中值存入Q中
					if (++Q.tu > MAXSIZE)return ERROR;
					Q.data[Q.tu] = { arow, ccol, ctemp[ccol] };
				}//if
		}//for arow
	}//if
}

该算法总的时间复杂度为O（M.mu * N.nu + M.tu * N.tu/N.mu）
ctemp初始化 O（M.mu * N.nu）
求Q中所有非零元 O（M.tu * N.tu / N.mu）
压缩存储 O（M.mu * N.nu）

头文件，宏定义，重命名

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

#define OK 1
#define ERROR -1
#define OVERFLOW -2
#define MAXSIZE 12500
 
typedef int ElemType;
typedef int Status;
//-----稀疏矩阵的三元组顺序表存储表示----
typedef	struct 
{
	int i,j;		//非零元的行下表和列下表 
	ElemType e;		//非零元的值 
}Triple;//三元组
typedef struct
{
	Triple data[MAXSIZE+1];//三元组表
	int mu,nu,tu;//矩阵行数，列数，非零元个数
}TSMatrix;

创建矩阵

挺长的
嗯，这个函数的功能
1.先输入矩阵的行数，列数，非零元个数，判断行数*列数是否大于非零元个数
2.三元组表的第0行用来存储矩阵的行数，列数，非零元个数
3.是否在相同位置重复输入
4.是否是按行非递减按列非递减输入
5.当输入次数大于等于非零元个数时退出循环，矩阵创建成功

Status CreateSMatrix(TSMatrix *M)
{
	int m=0;
	int k=0;
	while(1)
	{
		printf("请输入行数 列数 非零元个数\n");
		scanf("%d%d%d",&M->mu,&M->nu,&M->tu);
		if((M->mu*M->nu) < M->tu)
		{
			printf("输入错误，非零元素个数要小于等于行数乘列数，请重新输入。\n");
		}else
		{
			break;
		}	
	}
	(*M).data[k].i=M->mu;
	(*M).data[k].j=M->nu;
	(*M).data[k].e=M->tu;
	k++;
	while(1)
	{
		
		printf("请输入行下标 列下标 元素值\n");
		scanf("%d%d%d",&(*M).data[k].i,&(*M).data[k].j,&(*M).data[k].e);
		
		if(k>=2)
		{
			for(m=1;m<k;m++)
			{
				if((*M).data[k].i==(*M).data[m].i&&(*M).data[k].j==(*M).data[m].j)
				{
					printf("输入错误，输入的下标重复，请重新输入\n");
					k--;
					break;
				}
			}
			for(m=1;m<k;m++)
			{
				if(((*M).data[k].i==(*M).data[m].i&&(*M).data[k].j<(*M).data[m].j)||(*M).data[k].i<(*M).data[m].i)
				{
					printf("输入错误，下标输入时要递增输入，请重新输入！\n");
					k--;
					break;
				}
			}
		}
		k++;
		if(k>=(*M).tu+1)				
		{
			printf("输入完成\n");
			break;
		}	
	}
	return OK;
}

销毁矩阵

把三元组表的存储空间给释放了，行数，列数，非零元个数清零

//销毁矩阵
Status DestroySMatrix(TSMatrix *M)
{
	free((*M).data);
	(*M).mu=0;
	(*M).nu=0;
	(*M).tu=0;
	return ERROR;
}

输出矩阵

双重循环，当三元组表的行列都分别与循环变量m,n相等时，打印出三元组的该位置的非零元，否则，输出0

//输出矩阵
Status PrintSMatrix(TSMatrix M)
{
	int m,n,k=1;
	
	for(m=1;m<=M.mu;m++)
	{
		for(n=1;n<=M.nu;n++)
		{
			if(M.data[k].i==m&&M.data[k].j==n)
			{
				printf("%d ",M.data[k].e);
				k++;
			}else
			{
				printf("0 ");
			}
		}
		printf("\n");
	}
	return OK;	
}

普通转置

M为原矩阵，T为转置矩阵
M,T的行列互换，令T->tu = M.tu；
如果T->tu不等于0，q为T的三元组表当前的存储位置
算法描述如下
遍历M的三元组表（内循环），如果存在M.data[p].j 等于1（1是外循环的col），那么，给T赋值，即第一次循环，找出M中列数等于1的元素，赋值给T，第二次循环，找出M中列数等于2的元素，赋值给T，第三次，找出M中列数等于3的元素，赋值给T······

Status TransposeSMatrix(TSMatrix M, TSMatrix *T) 
{	
	T->mu = M.nu;
	T->nu = M.mu;
	T->tu = M.tu;
	int col,p,q;
	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;
}

快速转置

num[col]数组表示矩阵M中的第col列中非零元的个数
cpot[col]表示M中第col列的第一个非零元在T的三元组表中的位置
如果读者不懂快速转置的算法，建议去找该算法学习，笔者水平有限，只能写出代码，却无法讲出具体的算法，略为愚钝了

Status FastTransposeSMatrix(TSMatrix M, TSMatrix *T)
{
	T->mu=M.nu;
	T->nu=M.mu;
	T->tu=M.tu;
	int col,t,p,q;
	int num[MAXSIZE];
	int cpot[MAXSIZE];
	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.tu;++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;
}

完整源码

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

#define OK 1
#define ERROR -1
#define OVERFLOW -2
#define MAXSIZE 12500
 
typedef int ElemType;
typedef int Status;
//-----稀疏矩阵的三元组顺序表存储表示----
typedef	struct 
{
	int i,j;		//非零元的行下表和列下表 
	ElemType e;		//非零元的值 
}Triple;
typedef struct
{
	Triple data[MAXSIZE+1];
	int mu,nu,tu;
}TSMatrix;

Status CreateSMatrix(TSMatrix *M);//创建矩阵
Status DestroySMatrix(TSMatrix *M);//销毁矩阵
Status PrintSMatrix(TSMatrix M);//输出矩阵
Status TransposeSMatrix(TSMatrix M,TSMatrix *T);//转置矩阵
Status FastTransposeSMatrix(TSMatrix M, TSMatrix *T);//快速转置矩阵
int main()
{
	TSMatrix M,T,Qt;
	int i;
	int flag;
	printf("1.创建矩阵\n");
	printf("2.销毁矩阵\n");
	printf("3.输出矩阵\n");
	printf("4.转置矩阵\n");
	printf("5.快速转置矩阵\n");
	printf("输入0退出程序\n");
	while(1)
	{
		printf("请输入数字\n");
		scanf("%d",&i);
		switch(i)
		{
			case 1:
				flag=CreateSMatrix(&M);
				if(flag==1)
					printf("创建成功\n");
				else
					printf("创建失败\n");	
				break;
			case 2:
				if(flag==1)
				{
					flag=DestroySMatrix(&M);
					if(flag==ERROR)
					{
						printf("销毁成功\n");
						break;
					}else
					{
						printf("销毁失败\n");
					}	
				}else
				{
					printf("矩阵不存在\n");
				}
				break;
			case 3:
				if(flag==1)
				{
					printf("矩阵如下\n");
					PrintSMatrix(M);
				}else
				{
					printf("矩阵不存在\n");
				}
				break;
			case 4:
				if(flag==1)
				{
					printf("转置矩阵如下\n");
					TransposeSMatrix(M,&T);
					PrintSMatrix(T);
				}else
				{
					printf("矩阵不存在\n");
				}
				break;
			case 5:
				if(flag==1)
				{
					printf("快速转置矩阵如下\n");
					FastTransposeSMatrix(M,&Qt);
					PrintSMatrix(Qt);
				}else
				{
					printf("矩阵不存在\n");
				}
				break;
			default:
				printf("输入错误,请重新输入\n");
				break;					
		}
	}
}
//创建矩阵
Status CreateSMatrix(TSMatrix *M)
{
	int m=0;
	int flag=1;
	int k=0;
	while(1)
	{
		printf("请输入行数 列数 非零元个数\n");
		scanf("%d%d%d",&M->mu,&M->nu,&M->tu);
		if((M->mu*M->nu) < M->tu)
		{
			printf("输入错误，非零元素个数要小于等于行数乘列数，请重新输入。\n");
		}else
		{
			break;
		}	
	}
	(*M).data[k].i=M->mu;
	(*M).data[k].j=M->nu;
	(*M).data[k].e=M->tu;
	k++;
	while(flag)
	{
		
		printf("请输入行下标 列下标 元素值\n");
		scanf("%d%d%d",&(*M).data[k].i,&(*M).data[k].j,&(*M).data[k].e);
		
		if(k>=2)
		{
			for(m=1;m<k;m++)
			{
				if((*M).data[k].i==(*M).data[m].i&&(*M).data[k].j==(*M).data[m].j)
				{
					printf("输入错误，输入的下标重复，请重新输入\n");
					k--;
					break;
				}
			}
			for(m=1;m<k;m++)
			{
				if(((*M).data[k].i==(*M).data[m].i&&(*M).data[k].j<(*M).data[m].j)||(*M).data[k].i<(*M).data[m].i)
				{
					printf("输入错误，下标输入时要递增输入，请重新输入！\n");
					k--;
					break;
				}
			}
		}
		k++;
		if(k>=(*M).tu+1)				
		{
			printf("输入完成\n");
			break;
		}	
	}
	return OK;
}
//销毁矩阵
Status DestroySMatrix(TSMatrix *M)
{
	free((*M).data);
	(*M).mu=0;
	(*M).nu=0;
	(*M).tu=0;
	return ERROR;
}
//输出矩阵
Status PrintSMatrix(TSMatrix M)
{
	int m,n,k=1;
	
	for(m=1;m<=M.mu;m++)
	{
		for(n=1;n<=M.nu;n++)
		{
			if(M.data[k].i==m&&M.data[k].j==n)
			{
				printf("%d ",M.data[k].e);
				k++;
			}else
			{
				printf("0 ");
			}
		}
		printf("\n");
	}
	return OK;	
}
//转置矩阵
Status TransposeSMatrix(TSMatrix M, TSMatrix *T) 
{	
	T->mu = M.nu;
	T->nu = M.mu;
	T->tu = M.tu;
	int col,p,q;
	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)
{
	T->mu=M.nu;
	T->nu=M.mu;
	T->tu=M.tu;
	int col,t,p,q;
	int num[MAXSIZE];
	int cpot[MAXSIZE];
	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.tu;++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;
}

笔者无法讲解出这些算法的玄妙之处，无法为各位提供有价值的讲解，不能够进行答疑解惑，只能提供笔者所写的完整代码供读者借鉴，略感遗憾