目录
代码
例子

代码
import numpy as np

def jiedianmatrix(B, X, n):
    n1 =B.shape[0]
    Y = np.mat(np.zeros((n, n)))
    Y = complex(0)*Y
    for i in range(n1):
        if X[i, 1] != 0:
            p = int(X[i, 0])-1
            Y[p, p] = 1/X[i, 1]
    for k in range(n1):
        if B[k, 5] == 0:
            k1 = int(B[k, 0].real)-1
            k2 = int(B[k, 1].real)-1
        else:
            k2 = int(B[k, 0].real)-1
            k1 = int(B[k, 1].real)-1
        Y[k1, k2] = Y[k1, k2] - 1/(B[k, 2]*B[k, 4])
        Y[k2, k1] = Y[k1, k2]
        Y[k2, k2] = Y[k2, k2] + 1/(B[k, 2]*B[k, 4]**2) + B[k, 3]/2
        Y[k1, k1] = Y[k1, k1] + 1/B[k, 2] + B[k, 3]/2
    for i in range(n1):
        for j in range(n1):
            Y[i, j] = np.round(Y[i, j], 4)
    return Y


例子
if __name__ == "__main__":
    n = 5
    B = np.mat([[1, 2, 0.03j,      0,    1.05, 0],
                [2, 3, 0.08+0.3j,  0.5j, 1,    0],
                [3, 4, 0.015j,     0,    1.05, 1],
                [2, 5, 0.1+0.35j,  0,    1,    0],
                [3, 5, 0.04+0.25j, 0.5j, 1,    0]])
    X = np.mat([[1, 0],
                [2, 0],
                [3, 0],
                [4, 0],
                [5, 0]])
    Y = jiedianmatrix(B, X, n)
    print('{:12}{:10}{:10}{:8}{:7}{:8}'.format('节点1', '节点2',
                                '阻抗值', '线路电钠',
                                '变压器支路变比',
                                '第一个节点是否为高压端（是1，否0）'))
    print(B)
    print('形成的节点导纳矩阵为：')
    print(Y)


形成的节点导纳矩阵如下图所示：

节点导纳矩阵作为电力系统分析计算中广泛采用的一种数学模型，是电力系统分析计算的基础，其计算的正确与否有着重要的意义。文中详细分析了在实际工程应用中遇到的三序网络节点和支路编号、YN/d接线方式的变压器在零序导纳矩阵计算中的处理及多回线路的处理三个问题，并对节点导纳矩阵的程序实现方法进行了改进。经过验证，该程序能正确计算电网模型的三序节点导纳矩阵，并已成功应用于多个实际工程。
电力网络的数学模型是现代电力系统分析的基础，电力网络通常是由相应的节点导纳矩阵或节点阻抗矩阵来描述的。在电力系统分析中，我们需要面对成千上万个节点及电力网络所连接的电力系统，因此，对电力网络的描述和处理往往成为解决问题的关键。
以阻抗矩阵为基础的节点方程可以由注入电流直接求出电力网络各节点的电压，曾在电力系统计算中获得广泛应用。但由于节点阻抗矩阵是满矩阵，其存储规模及形成所需的计算量都随电网规模增大而呈平方增加，要求计算机有较大的内存容量，同时也增加了运算次数，降低了计算速度，因而当系统规模变大时，节点阻抗矩阵的应用受到一定限制。
而电力网络的节点导纳矩阵具有良好的稀疏特性，其存储规模和计算量一般与电网的规模呈线性关系，通过稀疏技术的引入，计算速度相比于阻抗矩阵模型大幅度降低，同时考虑到节点导纳矩阵在处理复杂电网结构和复杂故障方面的灵活性，可以高效地处理电力网络方程，因此是电力系统分析中广泛应用的数学模型。
1节点导纳矩阵的程序实现
节点导纳矩阵是用来描述电力系统网络模型的最基本的矩阵，由网络的结构和网络元件的参数决定。它是电力系统分析计算及相关高级应用软件的计算基础，因此，导纳矩阵的正确与否对电力系统分析计算及相关高级应用软件的计算结果准确性有着举足轻重的意义。
通过扫描网络中的支路，并根据支路在网络中的连接关系直接形成自导纳和互导纳，由于节点导纳矩阵为稀疏矩阵且是对称矩阵，可以利用稀疏矩阵技术进行存储，只存储导纳矩阵的上三角部分。节点导纳矩阵的程序实现如图1所示。
图中的YII、YIJ数组分别用来存放各节点的自导纳和互导纳；NYseq数组按导纳矩阵行号顺序存放各行非对角元素的首地址；NYsum数组用来存放导纳矩阵中各行非对角元素的个数。支路状态属性0代表停运，在计算中不考虑；1代表线路支路；2代表YN/Yn接法的变压器；3代表YN/d且运行的变压器或其他接法的变压器；4代表对地支路；5代表发电机支路；6代表Yn/d接法且停运的变压器。其中，支路状态为3和6的支路计算时比较特殊，其正、负序导纳矩阵和零序导纳矩阵的形成不同。
实际上，在零序网络中，YN/d接线的变压器支路的零序导纳只对YN侧节点的自导纳有影响，也就是YN/d变压器相当于一条对地支路。对于支路状态为6的支路在正、负序导纳矩阵计算中由于处于停运状态不进行计算，但在零序导纳矩阵计算中却要考虑。
图1节点导纳矩阵的程序实现
2 节点导纳矩阵计算中应注意的问题
2.1三序网络节点和支路编号
当电力系统处于不对称状态下，会出现正序、负序和零序分量，此时电力网络方程就是三序网络方程。三序网络方程是正序网络、负序网络、零序网络以及它们之间相互关系的电气特性的数学描述，它是电力系统分析和计算的理论基础。三序网络方程的系数矩阵就是三序导纳矩阵。因此，三序网络分析是导纳矩阵的计算的前提。
对于正序网和负序网，其网络结构相同，只是在计及发电机和负荷阻抗时，正序和负序导纳矩阵的对角元素不相同而已；对于正序网和零序网，其不同之处有两点：一是零序参数不同于正序参数，二是由于零序网络与变压器绕组接线方式和中性点接地方式有密切的关系，网络结构与正序网存在较大差异，零序网络通常比正序网规模小，因此，存在序网间节点和支路编号不对应问题。
为了解决这个问题，通常有两种解决方法：一种是正序网和零序网的编号各自独立进行，在编号完成后建立序间节点、支路编号的检索，节省储存零序导纳矩阵的内存开销，但由于正序网和零序网节点号不同，会给计算带来许多麻烦；一种是零序网采用与正序网完全一致的编号，这样，零序网络会出现一些空节点，即不连接任何支路的节点，在编写网络方程求解程序时，应增加判别空节点的功能，避免运算出错。
2.2 含YN/d接线方式变压器在零序导纳矩阵计算中的处理
在电力网络拓扑分析时，根据电网中各设备之间的连接关系将设备模型处理为两类：节点类和支路类，通常将网络中的输电线路、变压器和对地支路统称为支路。输电线路及对地支路在正序、负序和零序网络结构中是相同的，但变压器在三序网络结构中却受自身接线方式和中性点接地方式的影响。
文中采用的网络模型中支路类包含以下属性：首节点编号、末节点编号、支路状态、支路正序、负序、零序阻抗，其中支路状态中0代表停运，1代表投运的线路， 2为YN/yn接线的变压器，3为YN/d接线的变压器，其他接线的变压器填2或3均可，但零序阻抗为0。
通过在实际中的应用发现，变压器支路的状态属性用0、2、3这三种状态表示还不完全。
如图2(a)所示的简单网络，按照各设备元件的正序、零序阻抗的等值电路分析[8]，其正序、零序网络结构如图2(b，c)所示。
以图中的变压器T4为例，对于YN/yn/d接线方式的变压器可以等效为三个YN/yn，YN/yn，YN/d接线方式的双绕组变压器，若T43支路处于停运状态，将该支路的状态属性中置为0，此时形成的正序网络会有所改变，如图3，而零序网络不变。
按照上述计算导纳矩阵的方法，在形成零序网络时，在判断支路状态为0时，将停运支路剔除，因此，对于节点4会少计算一条支路形成的自导纳，造成零序导纳矩阵计算的不准确。
为了解决这个问题，为此，在进行网络分析形成支路类数据时，对于变压器的支路状态描述只有0、2、3三种状态是不够的，需要新增一个状态类型代表该支路为YN/d接线方式的变压器且处于停运状态。
这样在形成节点导纳矩阵时，通过判断支路的状态属性就不会造成支路丢失，保证计算的正确性。
图2序网图
2.3两节点间多回线路的处理
在实际电网模型中，常常存在两节点间有多回线路相连的情况，手工计算时可以轻易的找出两节点间的多回线路，而在计算机程序实现过程中，需要通过比较电网支路模型的首、末节点属性来判断是否为多回线路。
当为多回线路时，首、末节点的自导纳及首末节点间的互导纳需要修改，但在统计节点导纳矩阵对应行的非对角元个数时需要注意，由于多回线路只是影响自导纳和互导纳的计算，并没有造成非对角元个数的增加，因此，在程序中当判断存在多回线路时，NYsum数组中对应该行的非对角元个数不增加。
3 实例
如图2(a)所示的网络，各元件的参数如下(各参数均为标幺值，100MVA为基准值)：
线路L1：正序电抗=j0.0343，零序电抗=j0.0796；线路L2：正序电抗=j0.255，零序电抗=j0.0591；线路L3：正序电抗=j0.218，零序电抗=j0.6407；线路L4：正序电抗=j0.178，零序电抗=j0.5304；线路L5：正序电抗=j0.19，零序电抗=j0.5658；线路L6：正序电抗=j0.0131，零序电抗=j0.03768；线路L7：正序电抗=j0.037，零序电抗=j0.113；线路L8：正序电抗=j0.037，零序电抗=j0.113；变压器T1：正序电抗=j0.0217，零序电抗=j0.0503；变压器T2：正序电抗=j0.064，零序电抗=j0.134；变压器T3：正序电抗=j0.0124，零序电抗=j0.049；变压器T5：正序电抗=j0.0124，零序电抗=j0.049；变压器T4：正序电抗分别为j0.018，-j0.002，j0.0553，零序电抗分别为j0.018，-j0.002，j0.0553；变压器T6：正序电抗分别为j0.01，j0.001，j0.0337，零序电抗分别为j0.01，j0.001，j0.0337；发电机A：正序电抗=j0.15，零序电抗=j0.05；发电机C：正序电抗=j0.075；零序电抗=j0.025。
其中，各元件的负序参数与正序参数相同。通过程序计算得出的结果如下表1-3。
表1 导纳矩阵对角元
表2 导纳矩阵非对角元
表3 数组NYsum、NYseq
4 结论
节点导纳矩阵作为电力系统分析计算中广泛应用的数学模型，其计算结果的正确性有着举足轻重的意义。根据在多个实际工程中的应用，发现三序节点导纳矩阵计算中应注意的问题，通过对这些问题的分析，提出了一种实用化的节点导纳矩阵的程序实现方法。经过多次的验证，该程序能够正确的计算电网的节点导纳矩阵，为电力系统分析计算做好基础。
(本文选编自《电气技术》，作者为张洪波、付小鹏。)

电力系统分析实验课接近尾声，整理一下大三这一年做过的电力系统分析实验程序。 

这是第一个实验–生成节点导纳矩阵： 

实验数据如下 

 

第一、二列为节点编号，第三列为线路电阻，第四列为线路电抗，第五列为线路对地电导，如果是变压器支路，第五列为0，第六列为电压器变比。 

实验图如下：



function Y = BuildY(linedata)
    if nargin < 1
        linedata = [
            1   4   0.1     0.4     0.01528 1;
            1   3   0       0.3     0       1.1;
            1   2   0.12    0.5     0.01920 1;
            4   2   0.08    0.40    0.01413 1];
    end
    nf = linedata(:,1);
    nt = linedata(:,2);
    r = linedata(:,3);
    x = linedata(:,4);
    b = linedata(:,5);
    k = linedata(:,6);
    branchnum = length(nf);
    busnum =  max([nf',nt']);
    y = ones(branchnum,1)./(r+1j*x);
    Y = zeros(busnum);

    for n = 1:branchnum
        Y(nf(n),nt(n)) = Y(nf(n),nt(n))-k(n)*y(n);
        Y(nt(n),nf(n)) = Y(nf(n),nt(n));

        Y(nf(n),nf(n)) = Y(nf(n),nf(n))+k(n)^2*y(n)+1j*b(n);
        Y(nt(n),nt(n)) = Y(nt(n),nt(n))+y(n)+1j*b(n);

    end
end

最近学习电力系统分析这门专业课，发现计算机分析在这门课上非常重要。

大电网的等值电路参数计算、导纳矩阵存储、潮流计算等都需要在电脑上编程实现。

打算试试用拿手的纯C来实现导纳矩阵的存储，也就是实现一个稀疏矩阵的数据结构，可以完成建立、运算、修改、删除等操作。

最先想到的方法是用图论算法里用的很多的邻接表，但是发现邻接表在删除和增加行列上明显不方便，翻阅数据结构书籍发现十字链表很适合存储表示稀疏矩阵。

尝试一波…2018/11/22

 

开发中遇到的问题：

2018/11/23

关于calloc函数的段错误/溢出问题，申请时将矩阵类型和矩阵指针类型弄反，导致申请的空间小于预期。

对申请的堆空间进行了越界操作，并没有立即引发错误，但之后再次通过calloc函数申请堆空间时出现了段错误。

第一次遇到这个问题，标记一下，申请空间时一定要注意大小是否足够，以及要及时释放！

2018/11/24

矩阵变换牵动三部分，行/列指针连接，元素信息，矩阵信息。缺一不可。

2018/11/26

尺取法对链表遍历进行各种操作有很大帮助。（起到双向链表的部分功能）

如果使用双向链表，好处是删除链表节点时，行遍历到目标节点，列指针维护连通性就会很方便，可以直接得到上方元素的地址，这会减少一部分列指针维护的时间。

坏处是每个非零元素的存储量会增加，每次添加新节点时需要额外的赋值，逻辑更为复杂些。暂且不考虑换用双向十字链表。

 

1.头文件 SparseMatrix.h

预设函数返回状态、存储结构体类型、想要实现的函数。

/*稀疏矩阵——采用十字链表*/
/*Author:Belous*/
/*Date:2018/11/22*/
/*Date:2018/11/26*/

#ifndef SparseMatrix_H
#define SparseMatrix_H

#define TRUE 1
#define FALSE 0 
#define OK 1
#define ERROR 0
#define INF -1
#define OVERFLOW -2
#define UNKNOWN -3
/*函数结果状态代码（预设）*/ 

//#define SMatrix_Help
/*开启提示*/ 


typedef int Status;
/*函数返回状态类型*/
typedef int ElemType;
#define ElemSp "%i"
/*节点数据类型*/

typedef struct OLNode
{
	int x,y;
	ElemType value;
	struct OLNode *right,*down;
}OLNode;
typedef OLNode *OLink;
/*稀疏矩阵的元素体*/

typedef struct
{
	OLink *rhead,*chead;
	int row,column,unzero;
}CrossList;
/*稀疏矩阵结构体*/ 

typedef CrossList *CrossLink;
/*****矩阵指针类型 CrossLink *****/ 

Status CreateSMatrix(CrossLink *M);
/*创建，通过stdin读取数据*/
Status DestorySMatrix(CrossLink *M);
/*销毁*/
Status PrintSMatrix(CrossLink M);
/*打印，通过stdout输出数据，矩阵指针只读模式*/ 
Status CopySMatrix(CrossLink *M,CrossLink T);
/*复制,将稀疏矩阵T复制给M*/
Status ReferSMatrix(CrossLink *M,CrossLink N);
/*引用，将 M 矩阵指针指向 N 矩阵*/ 
Status AddMatrix(CrossLink M,CrossLink N);
Status AddSMatrix(CrossLink *M,CrossLink N,CrossLink *T);
/*矩阵加法*/ 
Status SubMatrix(CrossLink M,CrossLink N);
Status SubSMatrix(CrossLink *M,CrossLink N,CrossLink *T);
/*矩阵减法*/
Status MulSMatrix(CrossLink M,CrossLink N,CrossLink *T);
/*矩阵乘法*/
Status TranSMatrix(CrossLink M,CrossLink *T);
/*矩阵转置*/
Status InvSMatrix(CrossLink M,CrossLink *T);
/*矩阵求逆*/
Status ChangeMatrix(CrossLink M,int row,int column,ElemType value);
/*修改、添加矩阵元素*/ 
Status ReadMatrix(CrossLink M,int row,int column,ElemType *value);
/*读取矩阵元素,矩阵指针只读模式*/ 
Status DelMatrix(CrossLink M,int row,int column);
/*删除矩阵元素*/ 
Status DelRowMatrix(CrossLink M,int row);
/*删除矩阵行*/ 
Status DelColMatrix(CrossLink M,int column);
/*删除矩阵列*/
Status NewRowMatrix(CrossLink M,int row);
/*添加矩阵行*/
Status NewColMatrix(CrossLink M,int column);
/*添加矩阵列*/
#endif


2. InitSMatrix函数（内部函数）

用于为矩阵指针申请一个矩阵空间。

static Status InitSMatrix(CrossLink *M)
{
	if(!(*M=(CrossLink)calloc(1,sizeof(CrossList))))
		return OVERFLOW;
	/*通过指针给矩阵指针申请空间*/ 
	return OK;
}
/*内部函数：矩阵指针初始化函数 InitSMatrix */

3. CreateSMatrix 函数

得到一个矩阵指针地址，如果矩阵指针不为空，先通过 DestorySMatrix 函数销毁矩阵指针指向的空间。

通过 InitSMatrix 函数为矩阵指针申请一个新的稀疏矩阵空间。

从stdin标准输入流得到矩阵数据，输入零元素直接跳过，出现越界等错误立即返回异常代码。

元素x、y值可不按顺序输入。

Status CreateSMatrix(CrossLink *M)
{
	if(*M!=NULL)
		DestorySMatrix(M);
	/*矩阵指针不为空，通过DestorySMatrix函数销毁矩阵指针所指空间*/ 
	if(InitSMatrix(M)==OVERFLOW)
		return OVERFLOW;
	/*通过 InitSMatrix 函数创建新的矩阵指针*/ 
	#ifdef SMatrix_Help
		printf("输入稀疏矩阵的行数、列数\n");
	#endif
	scanf("%d%d",&((*M)->row),&((*M)->column));
	if(!((*M)->rhead=(OLink*)calloc((*M)->row+1,sizeof(OLink))))
		return OVERFLOW;
	if(!((*M)->chead=(OLink*)calloc((*M)->column+1,sizeof(OLink))))
		return OVERFLOW; 
	(*M)->unzero=0;
	int tx,ty;
	ElemType tv;
	#ifdef SMatrix_Help
	printf("输入非零元，x y value ，x是0时结束输入\n");
	#endif
	while(scanf("%d%d"ElemSp,&tx,&ty,&tv),tx!=0)
	{
		if(tv==(ElemType)0)
			continue;
		/*只有非零元素才需要存储！*/ 
		if(tx<=0||ty<=0||tx>((*M)->row)||ty>((*M)->column))
			return ERROR;
		/*判断输入非零元的x、y有效性*/ 
		OLink p;
		if(!(p=(OLink)calloc(1,sizeof(OLNode))))
			return OVERFLOW;
		p->x=tx,p->y=ty,p->value=tv;
		if((*M)->rhead[tx]==NULL||(*M)->rhead[tx]->y>ty)
		/*检查横向链表的首地址y是否在当前ty的右侧*/
		{
			p->right=(*M)->rhead[tx];
			(*M)->rhead[tx]=p;
		}
		else
		{
			OLink q;
			for(q=(*M)->rhead[tx];(q->right)&&q->right->y<ty;q=q->right);
			/*寻找行表中插入的位置*/
			p->right=q->right;
			q->right=p;
			/*插入到行表当前位置的右侧*/ 
		}
		if((*M)->chead[ty]==NULL||(*M)->chead[ty]->x>tx)
		/*检查横向链表的首地址x是否在当前tx的下侧*/
		{
			p->down=(*M)->chead[ty];
			(*M)->chead[ty]=p;
		}
		else
		{
			OLink q;
			for(q=(*M)->chead[ty];(q->down)&&q->down->x<tx;q=q->down);
			/*寻找列表中插入的位置*/
			p->down=q->down;
			q->down=p;
			/*插入到列表当前位置的下侧*/ 
		}
		(*M)->unzero++;
		/*完成一次非零元插入*/ 
	}
	return OK;
}

4. DestorySMatrix 函数

得到一个矩阵指针地址，检查矩阵指针是否已为空，为空则返回。

行指针遍历所有非空元素，逐一释放。

通过free函数释放矩阵指针所指的堆空间。

Status DestorySMatrix(CrossLink *M)
{
	if(*M==NULL)
		return OK;
	/*矩阵指针已空，无需操作*/ 
	for(int i=1;i<=(*M)->row;i++)
	{
		OLink temp=(*M)->rhead[i];
		while(temp!=NULL)
		{
			OLink Death=temp;
			temp=temp->right;
			free(Death);
		}
	}
	/*释放所有矩阵元素*/ 
	free((*M)->rhead);
	free((*M)->chead);
	free(*M);
	*M=NULL;
	/*清除矩阵指针所指的原地址*/ 
	return OK;
}

5. ReadMatrix 函数

对矩阵元素读取返回，因此对矩阵是只读操作，得到矩阵指针、行号、列号、返回变量地址。

先检测行号、列号是否有头指针，若头指针为空，直接返回零元素。

采取行指针线性查找，查找到对应的列号则返回该非零元素，否则返回元素0。

Status ReadMatrix(CrossLink M,int row,int column,ElemType *value)
{
	if(M==NULL)
		return ERROR;
	/*判断矩阵指针是否有效*/ 
	if(row<=0||column<=0||row>M->row||column>M->column)
		return ERROR;
	/*判断参数是否越界*/
	if(M->chead[column]==NULL||M->rhead[row]==NULL)
	{
		*value=(ElemType)0;
		return OK;
	}
	/*行/列指针为空，不存在元素*/
	OLink temp=M->rhead[row];
	/*得到行指针*/
	while(temp->y<column&&temp->right!=NULL)
		temp=temp->right;
	/*线性查找行指针*/
	if(temp->y==column)
	{
		*value=temp->value;
		return OK;
	}
	/*找到矩阵非零元素*/
	*value=(ElemType)0;
	/*未找到，返回零元素*/ 
	return OK;
}

6. PrintSMatrix 函数

打印矩阵指针所指的稀疏矩阵所有元素、行数、列数、非空元素个数、打印所耗CPU时间（s）

Status PrintSMatrix(CrossLink M)
{
	if(M==NULL)
		return ERROR;
	/*检查矩阵指针是否有效*/ 
	clock_t begin=clock();
	/*记录函数启动时间*/
	printf("\nrow = %d ,column = %d ,unzero = %d\n\n",M->row,M->column,M->unzero);
	/*输出目标矩阵的行，列，非零元个数*/
	int i,j;
	ElemType temp;
	for(i=1;i<=M->row;i++)
		for(j=1;j<=M->column;j++)
		{
			ReadMatrix(M,i,j,&temp);
			/*读取i行j列的矩阵元素*/ 
			printf(ElemSp"%c",temp,j==M->column?'\n':' ');
		}
	printf("\nElapsed Time： %.3lf s\n\n",(clock()-begin)/(double)CLOCKS_PER_SEC);
	/*输出函数耗费CPU时间，单位s*/ 
	return OK;
}

7. CopySMatrix 函数

复制T矩阵给M矩阵。

T矩阵指针不得为空，立即终止函数并返回ERROR标识，此时复制也是完成的，都是空指针。

按行指针遍历逐一复制，并通过辅助指针维护列指针。

复制完成后须及时释放辅助指针。

Status CopySMatrix(CrossLink *M,CrossLink T)
{
	if(*M!=NULL)
		DestorySMatrix(M);
	/*矩阵指针不为空，通过DestorySMatrix函数销毁矩阵指针所指空间*/
	if(T==NULL)
		return ERROR;
	/*待复制矩阵为空，函数立即终止，返回ERROR*/ 
	if(InitSMatrix(M)==OVERFLOW)
		return OVERFLOW;
	/*通过 InitSMatrix 函数创建新的矩阵指针*/ 
	(*M)->row=T->row;
	(*M)->column=T->column;
	(*M)->unzero=T->unzero;
	if(!((*M)->rhead=(OLink*)calloc((*M)->row+1,sizeof(OLink))))
		return OVERFLOW;
	if(!((*M)->chead=(OLink*)calloc((*M)->column+1,sizeof(OLink))))
		return OVERFLOW;
	/*为rhead和chead申请空间*/
	OLink Rtemp=NULL;
	bool isRhead=false;
	OLink *Ctemp=NULL;
	if(!(Ctemp=(OLink*)calloc(T->column+1,sizeof(OLink))))
		return OVERFLOW;
	bool *isChead=NULL;
	if(!(isChead=(bool*)calloc(T->column+1,sizeof(bool))))
		return OVERFLOW;
	/*临时辅助列指针*/ 
	for(int i=1;i<=T->row;i++)
	/*行指针遍历矩阵T*/ 
	{
		isRhead=false;
		OLink temp=T->rhead[i];
		Rtemp=NULL;
		if(temp==NULL)
			continue;
		/*当前行头指针为空，没有元素*/
		while(temp!=NULL)
		{
			OLink p;
			if(!(p=(OLink)calloc(1,sizeof(OLNode))))
				return OVERFLOW; 
			p->x=temp->x;
			p->y=temp->y;
			p->value=temp->value;
			if(isRhead==false)
			{
				(*M)->rhead[i]=p;
				Rtemp=p;
				isRhead=true;
			}
			else
			{
				Rtemp->right=p;
				Rtemp=p;
			}
			if(isChead[p->y]==false)
			{
				(*M)->chead[p->y]=p;
				Ctemp[p->y]=p;
				isChead[p->y]=true;
			}
			else
			{
				Ctemp[p->y]->down=p;
				Ctemp[p->y]=p;
			}
			/*更新行、列辅助指针*/
			temp=temp->right;
			/*继续遍历右侧元素地址*/ 
		} 
	}
	free(Ctemp);
	free(isChead);
	/*释放辅助列指针*/ 
	return OK;
}

8. DelMatrix 函数

删除稀疏矩阵中指定位置的非零元素。

首先判断矩阵指针是否有效，row、column值是否正确。检测对应行、列头指针是否存在。

逐行寻找，未发现对应的非零元素则退出。找到后存下该节点的地址，行、列指针表均更新完毕后，释放该节点，矩阵非零元个数unzero减一。

Status DelMatrix(CrossLink M,int row,int column)
{
	if(M==NULL)
		return ERROR;
	/*判断矩阵指针是否有效*/ 
	if(row<=0||column<=0||row>M->row||column>M->column)
		return ERROR;
	/*判断参数是否越界*/
	OLink Rtemp=M->rhead[row];
	OLink Ctemp=M->chead[column];
	OLink temp=NULL;
	if(Rtemp==NULL||Ctemp==NULL)
		return ERROR;
	/*行/列指针为空，元素不存在*/
	if(Rtemp->y==column)
	/*在行头指针处发现非零元素*/ 
	{
		temp=Rtemp;
		M->rhead[row]=temp->right;
		/*行指针表中删除*/
	}
	else
	{
		while(Rtemp->right!=NULL&&Rtemp->right->y<column)
			Rtemp=Rtemp->right;
		if(Rtemp->right!=NULL&&Rtemp->right->y==column)
		/*在行指针表中找到目标位置的非零元素*/ 
		{
			temp=Rtemp->right;
			Rtemp->right=temp->right;
			/*行指针表中删除*/
		}
		else
			return ERROR;
	}
	if(Ctemp->x==row)
	/*在列头指针处！*/ 
	{
		M->chead[column]=temp->down;
		M->unzero--;
		free(temp);
		return OK;
	}
	while(Ctemp->down!=NULL&&Ctemp->down->x<row)
		Ctemp=Ctemp->down;
	Ctemp->down=temp->down;
	M->unzero--;
	free(temp);
	return OK;
}

9. ChangeMatrix 函数

检测矩阵指针、参数有效性后，判断要改变的值是否是 0 ，如为 0 则为删除对应位置的元素，调用 DelMatrix 函数后正常退出即可。

有两种可能，一种是更新已有的非零元素，另一个是在目标位置建立新的节点。

先默认是第二个情况，根据参数建立p节点。行指针遍历查找合适位置，存在已有元素则释放p节点，更新元素值后退出。未找到则更新行、列指针表，非零元个数unzero加一。

Status ChangeMatrix(CrossLink M,int row,int column,ElemType value)
{
	if(M==NULL)
		return ERROR;
	/*判断矩阵指针是否有效*/ 
	if(row<=0||column<=0||row>M->row||column>M->column)
		return ERROR;
	/*判断参数是否越界*/
	if(value==(ElemType)0)
	{
		DelMatrix(M,row,column);
		return OK;
	}
	/*待更改元素的值为 0 ，即清除该位置*/ 
	OLink Rtemp=M->rhead[row];
	OLink Ctemp=M->chead[column];
	OLink temp=Rtemp;
	/*得到行、列头指针，辅助指针*/ 
	OLink p;
	if(!(p=(OLink)calloc(1,sizeof(OLNode))))
		return OVERFLOW;
	p->x=row;
	p->y=column;
	p->value=value;
	/*预设需要建立新的元素节点*/ 
	if(temp==NULL)
		M->rhead[row]=p;
	/*行头指针为空*/ 
	else
	{
		if(temp->y==column)
		{
			temp->value=value;
			free(p);
			return OK;
		}
		/*在头指针找到元素值,更改后释放p节点退出*/
		else if(temp->y>column)
		{
			p->right=temp;
			M->rhead[row]=p;
		}
		/*放在头指针前*/
		else
		{
			while(temp->right!=NULL&&temp->right->y<column)
				temp=temp->right;
			/*右侧查找*/ 
			if(temp->right!=NULL&&temp->right->y==column)
			{
				temp->right->value=value;
				free(p);
				return OK;
			}
			/*找到元素值，更新后释放p节点退出*/
			p->right=temp->right;
			temp->right=p;
		}
	}
	temp=Ctemp;
	/*改为列辅助指针*/
	if(temp==NULL)
	{
		M->chead[column]=p;
		M->unzero++;
		return OK;	
	}
	if(temp->x>row)
	{
		p->down=temp;
		M->chead[column]=p;
		M->unzero++;
		return OK;
	}
	while(temp->down!=NULL&&temp->down->x<row)
		temp=temp->down;
	/*列指针表中找到正确位置插入*/ 
	p->down=temp->down;
	temp->down=p;
	M->unzero++;
	return OK;
}

10. DelRowMatrix 函数

删除指定行。

删除大体需要完成三部分操作，更新行指针表；维护列指针正确连接；更新变换坐标后元素的行号、矩阵信息（row、unzero）

Status DelRowMatrix(CrossLink M,int row)
{
	if(M==NULL)
		return ERROR;
	/*判断矩阵指针是否有效*/ 
	if(row<=0||row>M->row)
		return ERROR;
	/*判断参数是否越界*/
	OLink* cpRhead=M->rhead;
	/*用于更新行指针的辅助指针*/
	if(!(M->rhead=(OLink*)calloc(M->row+1,sizeof(OLink))))
		return OVERFLOW;
	/*申请新的行指针列表*/ 
	for(int i=1;i<row;i++)
		M->rhead[i]=cpRhead[i];
	if(cpRhead[row]!=NULL)
	{
		OLink temp=cpRhead[row];
		while(temp!=NULL)
		{
			OLink Death=temp;
			OLink connect=M->chead[temp->y];
			if(connect->x==row)
				M->chead[temp->y]=temp->down;
			else
			{
				while(connect->down!=NULL&&connect->down->x<row)
					connect=connect->down;
				connect->down=temp->down;
			}
			/*维护列指针正确连通*/ 
			temp=temp->right;
			free(Death);
		}
	}
	for(int i=row+1;i<=M->row;i++)
	{
		M->rhead[i-1]=cpRhead[i];
		OLink temp=M->rhead[i-1];
		while(temp!=NULL)
		{
			temp->x--;
			temp=temp->right;
		}
		/*更新元素行号*/ 
	}
	M->unzero-=M->column;
	M->row--;
	/*更新矩阵信息*/ 
	free(cpRhead);
	/*释放旧的行指针列表*/
	return OK; 
}

11. DelColMatrix 函数

删除指定列。

与 DelRowMatrix 函数原理类似。

Status DelColMatrix(CrossLink M,int column)
{
	if(M==NULL)
		return ERROR;
	/*判断矩阵指针是否有效*/ 
	if(column<=0||column>M->column)
		return ERROR;
	/*判断参数是否越界*/
	OLink* cpChead=M->chead;
	/*用于更新列指针的辅助指针*/
	if(!(M->chead=(OLink*)calloc(M->column+1,sizeof(OLink))))
		return OVERFLOW;
	/*申请新的列指针列表*/ 
	for(int i=1;i<column;i++)
		M->chead[i]=cpChead[i];
	if(cpChead[column]!=NULL)
	{
		OLink temp=cpChead[column];
		while(temp!=NULL)
		{
			OLink Death=temp;
			OLink connect=M->rhead[temp->x];
			if(connect->y==column)
				M->rhead[temp->x]=temp->right;
			else
			{
				while(connect->right!=NULL&&connect->right->y<column)
					connect=connect->right;
				connect->right=temp->right;
			}
			/*维护行指针连通*/ 
			temp=temp->down;
			free(Death);
		}
	}
	for(int i=column+1;i<=M->column;i++)
	{
		M->chead[i-1]=cpChead[i];
		OLink temp=M->chead[i-1];
		while(temp!=NULL)
		{
			temp->y--;
			temp=temp->down;
		}
		/*更新元素列号*/ 
	}
	M->unzero-=M->row;
	M->column--;
	/*更新矩阵信息*/ 
	free(cpChead);
	/*释放旧的行指针列表*/
	return OK; 
}

12. NewRowMatrix 函数

与 DelRowMatrix 函数作用相反，添加一个新的空行，参数就是新行的行号。

输入参数是 0 时，等效为在行末添加一行。

Status NewRowMatrix(CrossLink M,int row)
{
	if(M==NULL)
		return ERROR;
	/*判断矩阵指针是否有效*/ 
	if(row<0||row>M->row+1)
		return ERROR;
	/*判断参数是否越界*/
	OLink* cpRhead=M->rhead;
	/*用于更新行指针的辅助指针*/
	if(!(M->rhead=(OLink*)calloc(M->row+2,sizeof(OLink))))
		return OVERFLOW;
	/*申请新的行指针列表*/
	if(row==0)
		row=M->row+1;
	/*参数输入 0 ，等效添加末尾行*/ 
	for(int i=1;i<row;i++)
		M->rhead[i]=cpRhead[i];
	M->rhead[row]=NULL;
	for(int i=row+1;i<=M->row+1;i++)
	{
		M->rhead[i]=cpRhead[i-1];
		OLink temp=M->rhead[i];
		while(temp!=NULL)
		{
			temp->x++;
			temp=temp->right;
		}
		/*更新元素行号*/ 
	}
	M->row++;
	/*更新矩阵信息*/ 
	free(cpRhead);
	/*释放旧的行指针列表*/
	return OK; 
}

13. NewColMatrix 函数

与 NewRowMatrix 函数作用类似，添加新的空列。

Status NewColMatrix(CrossLink M,int column)
{
	if(M==NULL)
		return ERROR;
	/*判断矩阵指针是否有效*/ 
	if(column<0||column>M->column+1)
		return ERROR;
	/*判断参数是否越界*/
	OLink* cpChead=M->chead;
	/*用于更新列指针的辅助指针*/
	if(!(M->chead=(OLink*)calloc(M->column+2,sizeof(OLink))))
		return OVERFLOW;
	/*申请新的列指针列表*/
	if(column==0)
		column=M->column+1;
	/*参数输入 0 ，等效添加末尾列*/ 
	for(int i=1;i<column;i++)
		M->chead[i]=cpChead[i];
	M->chead[column]=NULL;
	for(int i=column+1;i<=M->column+1;i++)
	{
		M->chead[i]=cpChead[i-1];
		OLink temp=M->chead[i];
		while(temp!=NULL)
		{
			temp->y++;
			temp=temp->down;
		}
		/*更新元素列号*/ 
	}
	M->column++;
	/*更新矩阵信息*/ 
	free(cpChead);
	/*释放旧的列指针列表*/
	return OK; 
}

14. AddMatrix 函数

终于到了重头戏！！开始各种运算函数！！

功能：将矩阵 N 加到矩阵 M。

加法运算函数，矩阵加法很简单，对应元素相加即可。

十字链表形式的稀疏矩阵想要实现对应位元素相加可以很简单的做到，调用上面那些已经编写好的函数即可。

但是，为单个元素就调用一次函数会造成很大的时间浪费！不断的遍历同样的东西是不明智的，超大规模稀疏矩阵运算时这种问题就会显现而来。所以使用十字链表形式，就要用出这种形式的优势！

优势就在于只需要每行遍历一遍，并时刻维护列指针连通即可。

这么做自然算的飞快，但是逻辑会复杂很多。需要处理指针头、相加后变为零元素须释放、原本没有元素的位置需要插入新节点等等。

我采用了辅助指针锁定上个操作的地址，像尺取法，一次一次的将矩阵 N 中的非零元素通过各种操作 加 / 插入 / 删除 到M矩阵中。各类数据测试完毕——2018/11/26

Status AddMatrix(CrossLink M,CrossLink N)
{
	if(M==NULL||N==NULL)
		return ERROR; 
	if((M->row!=N->row)||(M->column!=N->column))
		return ERROR;
	/*矩阵不同维度，无法加法运算*/
	for(int i=1;i<=N->row;i++)
	{
		if(N->rhead[i]==NULL)
			continue;
		/*待加数当前行为空，无需计算*/
		OLink MRhead=M->rhead[i];
		OLink NRhead=N->rhead[i];
		OLink Last=MRhead;
		bool isRhead=false;/*是否对头指针进行更新操作*/ 
		/*两个矩阵的辅助行指针，准备开始当前行运算*/ 
		while(NRhead!=NULL)
		{
			if(M->rhead[i]==NULL)
			/*1.行指针整体为空！*/ 
			{
				OLink p;
				if(!(p=(OLink)calloc(1,sizeof(OLNode))))
					return OVERFLOW;
				/*申请非零元素空间*/
				p->x=i;
				p->y=NRhead->y;
				p->value=NRhead->value;
				M->rhead[i]=p;
				OLink temp=M->chead[p->y];
				/*列辅助指针*/
				if(temp==NULL)
				{
					M->chead[p->y]=p;
				}
				else if(temp->x>p->x)
				{
					p->down=temp;
					M->chead[p->y]=p;
				}
				else
				{
					while(temp->down!=NULL&&temp->down->x<p->x)
						temp=temp->down;
					/*列指针表中找到正确位置插入*/ 
					p->down=temp->down;
					temp->down=p;
				}
				Last=MRhead=p;
				M->unzero++;
				/*非零元素增加！*/
			}
			else if(M->rhead[i]->y>NRhead->y)
			{
				OLink p;
				if(!(p=(OLink)calloc(1,sizeof(OLNode))))
					return OVERFLOW;
				/*申请非零元素空间*/
				p->x=i;
				p->y=NRhead->y;
				p->value=NRhead->value;
				p->right=M->rhead[i];
				M->rhead[i]=p;
				OLink temp=M->chead[p->y];
				/*列辅助指针*/
				if(temp==NULL)
				{
					M->chead[p->y]=p;
				}
				else if(temp->x>p->x)
				{
					p->down=temp;
					M->chead[p->y]=p;
				}
				else
				{
					while(temp->down!=NULL&&temp->down->x<p->x)
						temp=temp->down;
					/*列指针表中找到正确位置插入*/ 
					p->down=temp->down;
					temp->down=p;
				}
				Last=MRhead=p;
				/*标记此次运算最后操作的地址*/ 
				M->unzero++;
				/*非零元素增加！*/
			}
			else
			{
				while(MRhead->right!=NULL&&MRhead->y<NRhead->y)
				{
					if(isRhead==false)
						isRhead=true;
					else
						Last=Last->right;
					MRhead=MRhead->right;
				}
				if(MRhead->y==NRhead->y)
				{
					MRhead->value+=NRhead->value;
					if(MRhead->value==(ElemType)0)
					{
						/*快速删除当前元素*/ 
						Last->right=MRhead->right;
						OLink temp=M->chead[MRhead->y];
						if(temp->x==MRhead->x)
						/*在列头指针处！*/ 
							M->chead[MRhead->y]=MRhead->down;
						else
						{
							while(temp->down!=NULL&&temp->down->x<MRhead->x)
								temp=temp->down;
							temp->down=MRhead->down;
						}
						M->unzero--;
						free(MRhead);
						MRhead=Last;
						/*指针左移，MRhead原内容会被释放*/ 
					}
					//相加导致零元素！
					Last=MRhead;
				}
				else
				{
					OLink p;
					if(!(p=(OLink)calloc(1,sizeof(OLNode))))
						return OVERFLOW;
					/*申请非零元素空间*/
					p->x=i;
					p->y=NRhead->y;
					p->value=NRhead->value;
					p->right=MRhead->right;
					OLink temp=M->chead[p->y];
					/*列辅助指针*/
					if(temp==NULL)
					{
						M->chead[p->y]=p;
					}
					else if(temp->x>p->x)
					{
						p->down=temp;
						M->chead[p->y]=p;
					}
					else
					{
						while(temp->down!=NULL&&temp->down->x<p->x)
							temp=temp->down;
						/*列指针表中找到正确位置插入*/ 
						p->down=temp->down;
						temp->down=p;
					}
					/************/
					MRhead->right=p;
					Last=MRhead=p;
					M->unzero++;
					/*非零元素增加！*/ 
				}
				isRhead=false;
				/*重新置位历史值*/ 
			}
			NRhead=NRhead->right;
		}
	}
	return OK;
}

15. ReferSMatrix 函数

引用函数，令矩阵指针 M 指向矩阵 N 。

使用时必须小心。这会释放 矩阵指针 M 原来所指的矩阵。矩阵 N 被释放时须同时清空 M ，否则再对 M 进行操作会导致越界（段错误）。

Status ReferSMatrix(CrossLink *M,CrossLink N)
{
	if(*M!=NULL)
		DestorySMatrix(M);
	/*矩阵指针不为空，通过DestorySMatrix函数销毁矩阵指针所指空间*/
	if(N==NULL)
		return ERROR;
	/*N矩阵为空，函数中止*/
	(*M)=N;
	return OK;
}

16. AddSMatrix 函数

将矩阵 M 矩阵 N 的和存入 矩阵 T 中，不改变矩阵 M 和 N 的值。

采用 ReferSMatrix 和 CopySMatrix 函数相配合，调用 AddMatrix 后恢复M的原值，相加后的矩阵给 T 。

也可以逐行挨个相加，但较为复杂。简单的实现计算速度太慢，还不如 AddMatrix +CopySMatrix 。

矩阵相加如果不再需要原 M 值就用 AddMatrix ，AddSMatrix 会比 AddMatrix 多出一个复制矩阵的时间。

Status AddSMatrix(CrossLink *M,CrossLink N,CrossLink *T)
{
	if(*T!=NULL)
		DestorySMatrix(T);
	/*存储和的矩阵不为空，调用 DestorySMatrix 函数释放*/ 
	CrossLink Temp=NULL;
	CopySMatrix(&Temp,*M);
	/*复制 M 矩阵*/
	AddMatrix(*M,N);
	/*调用 AddMatrix 函数，将矩阵 N 加在矩阵 M 上*/
	ReferSMatrix(T,*M);
	*M=Temp;
	/*为矩阵 M 恢复原始值*/
	return OK; 
}

.

17. SubMatrix 函数

18. SubSMatrix 函数

19. MulSMatrix 函数

20. TranSMatrix 函数

21. InvSMatrix 函数

 

 

END