矩阵的重要性常常被在校学生忽视。实际上，电力系统分析和仿真，离开矩阵这个工具，可以说寸步难行。电力系统所用到的矩阵，绝大多数场合都是稀疏的。也就是说，这些矩阵，大多数元素都是零，只有少部分元素是非零的。存储稀疏矩阵时，只需要存储非零元；计算稀疏矩阵时，零值元素都可以跳过去。从而大大节省存储空间、加快计算速度。电力系统计算，绝大多数时间都是消耗在稀疏矩阵的计算和处理上。

只要是实际存在的网络，都具有稀疏特性。除了电力系统，交通系统、互联网搜索(PageRank算法)、通信系统、人际交流网络，都必然是稀疏的。另外，在有限元分析、天气预报、石油探测等领域也主要应用稀疏矩阵。这是因为实际系统节点度（即连接某个节点的支路数量）都是相对来说是固定的，例如电力系统，连接节点的线路数量平均大约是6个。假定网络规模增大n倍，那么矩阵元素数量变为n的平方。但是非零元的数量只是与n同数量级增长。也就是说，网络规模增大n倍，稀疏度（非零元占比）往往下降n倍。大规模电力系统所对应的矩阵是高度稀疏的，当然，交通系统、通信系统等网络都是类似的。

正因为稀疏矩阵的广泛存在和重要性，数值分析和计算机专业的学者投入了大量精力来优化算法。电力系统分析软件的开发应该优先使用国际先进的稀疏矩阵软件包，而不是选择自己开发。同样是稀疏矩阵计算，国际先进的优化算法与外行编写的算法，有超过1个数量级的性能差别。我自己进行稀疏矩阵测试时，常常被不同程序性能的差距吓一跳。我参加的一次会议上，某顶尖高校的电力系统专业教授对自己写的稀疏矩阵程序有充分自信，觉得已经充分考虑稀疏性了，性能已经是足够优化了。但是电力系统的专家并不必然是数学和计算机专家，专业事情最好还是交给专业人士去做。

稀疏矩阵的求解从大的方面来说，有直接法和迭代法。其中直接法是通过把矩阵分解为多个特殊矩阵乘积（例如分解为两个三角矩阵乘积、或者分解为正交矩阵及上三角矩阵，还有其它分解方法），再进行后续处理。迭代法是不进行分解，而是通过多次迭代来求解方程组，迭代法适用于维度非常大的矩阵。但是对于电力系统来说，最大的电网也不超过10万阶，相对来说，属于不那么大的矩阵，不能充分发挥迭代法的优势。因此，电力系统的稀疏矩阵求解，都使用直接法。

在国际上，最受公认的开源稀疏矩阵直接求解库是德州农工大学TimDavis教授开发的SuiteSparse软件包，该软件包有许多稀疏矩阵求解程序，适应不同的场合，我后文提到的程序大部分都是属于SuiteSparse软件包的一部分。最著名的大概就是其中的UMFPACK库了。据说TimDavis教授此人和蔼可亲，有问必答，与工业界的关系极好。包括MATLAB中的稀疏矩阵计算都默认使用TimDavis编写的程序。我认识一位电力系统专业的老师想去德州农工大学进修访问，因为TimDavis在那里，当然德州农工的电力系统专业也很强。有兴趣的同行或同学可以读一下他写的《Direct Methods for Sparse Linear Systems》，深入浅出，包含很多示例代码。

高性能的稀疏矩阵直接求解，有两种主要思路：

1）不调用稠密矩阵的直接法，例如Left-looking方法，SuiteSparse软件包中的KLU属于这一类型；

2）调用稠密矩阵BLAS库的直接法。主要有：超节点方法（SuperLU软件包）、多波前法（SuiteSparse软件包中UMFPACK库）。核心思路是通过变换操作获取密集子阵，例如超节点、波前阵，再调用高性能的BLAS库来加快计算速度。

BLAS库是专门处理稠密矩阵基本计算的软件包，BLAS库是相当接近底层硬件的程序，因此硬件厂家针对特定设备编写的BLAS库是效率最优的，例如Intel的MKL，AMD的ACML，还包括一些针对GPU的BLAS库实现等（NVIDA的cuBLAS）。但是这样做就会导致安装过程繁琐，安装维护都无法统一。另外这些库大部分都不是开源的，或者商业上的使用需要授权。

开源的BLAS库主要有：

1）ATLAS库，是比较优化的开源BLAS实现，采用了自动优化技术（根据系统硬件情况自动选择或者生成源代码）。早期版本的MATLAB就采用ATLAS库。ATLAS虽然不能提供最优的效率，但在大多数平台上，都可以获得不错的性能。ATLAS的编译安装过程较复杂。

2）GotoBLAS库，针对一些特定的设备进行了手工汇编优化。在一些测试中，取得了非常高的成绩。缺点是GotoBLAS很久都不更新了，在一些平台上编译不能通过（因为采用很多汇编代码）。好在中科院在GotoBLAS基础上开发了OpenBLAS库，比GotoBLAS的适用范围广。

下面着重介绍一下OpenBLAS库，这个库是中国人开发维护的。其作者张先轶是中科院博士，到美国做博士后研究，现在创办了PerfXLab澎峰科技公司。下面的内容摘选自2014年张先轶博士参加的一次国内会议。

张先轶：谢谢各位，能够有机会介绍我们OpenBLAS项目。这BLAS到底是什么意思呢？他就是这四个英文，Basic、Linear、Algebra、Subprograms。这BLAS做什么呢？就是向量和矩阵的一些基本运算。

   BLAS库最重要的东西就是三点，性能、性能、性能！就是要快，在所有平台上比别人快我们就成功了。以目前来说，这位仁兄是一个日本人（KazushigeGoto），他做了一个GotoBLAS特别出名。这位先生，在2010年的时候去了微软，因为微软那边钱比较多。在微软没有多么久，去年就去了另外一家公司。

因为2010年这个库就停止了，我们在2011年又发起了OpenBLAS，我们的长远目标是最好的BLAS库开源实现。

主要做的事情是三个事情，一、X86平台，二、龙芯3a，三、ARM平台。虽然我们只有3个人，我们的结果还真不错。

在BLAS性能优化里边，有两个大流派，一个是手工优化，手工优化要做大量的汇编。比如说GotoBLAS，它需要写大量的汇编。我不知道大家有没有写过很多汇编的。

   另一个流派，这个典型库就是ATLAS，采用自动优化，自动优化的时候要花费一天的时间。

   我们能否把手工和自动结合起来，走第三条路，把手工的经验尽量自动化，这里边实际上就解决了自动的生成高效汇编。

   结论：我们性能是最好的，大家只有记住这一点就可以了。我们减轻了开发成本。我们在这里用了一些自动搜索的技术，很容易可以把参数调好。好，就到这里，谢谢各位！

回到电力系统稀疏矩阵计算，前文已经提过，随着系统规模的扩大，越来越稀疏。对于大型电力系统来说，很难生成密集子阵，也就是说，调用BLAS库加速的机会变小。因此，电力系统比较适用KLU这样的不依赖于BLAS稀疏矩阵库。对于一些算例，KLU的性能超出UMFPACK一倍以上。

现代处理器已经是CPU多核并行的时代了，通用图形处理器(GPU)也得到了广泛应用。怎样利用现代处理器的并行特性来加速稀疏矩阵计算，也得到了广泛研究。当然，SuiteSparse软件包在最新的更新中，也逐渐加入了适用GPU的代码。当然，SuiteSparse毕竟诞生于单核时代，并行化改造才刚刚开始。

关于利用并行特性加速稀疏矩阵求解，我有以下看法：

1）还是优先使用库作者的代码，而不是自己重新写。例如SuperLU_MT就是适用于多线程共享内存环境的程序包。例如SuiteSparse中的QR分解已经加入了GPU特性。

2）对于特殊场合，不能直接把库拿来用的情况下。可自己在库的基础上进行并行化修改，例如KLU是单线程程序，用OpenMP把它改造成多线程执行并不困难。但是，如果充分发挥并行优势，还是需要对计算机和应用问题都有深刻的理解才行。例如如何将大矩阵分块、并行求解、结果的同步协调汇总。