矩阵LU分解分块算法实现

本文主要描述实现LU分解算法过程中遇到的问题及解决方案，并给出了全部源代码。

1. 什么是LU分解？

         矩阵的LU分解源于线性方程组的高斯消元过程。对于一个含有N个变量的N个线性方程组，总可以用高斯消去法，把左边的系数矩阵分解为一个单位下三角矩阵和一个上三角矩阵相乘的形式。这样，求解这个线性方程组就转化为求解两个三角矩阵的方程组。具体的算法细节这里不做过多的描述，有很多的教材和资源可以参考。这里推荐的参考读物如下：

 Numerical recipes C++,还有包括MIT的线性代数公开课。

2. LU分解有何用？

    LU分解来自线性方程组求解，那么它的直接应用就是快速计算下面这样的矩阵乘法

       A^(-1)*B，这是矩阵方程 AX=B 的解
       A^(-1)*b，这是线性方程组 Ax=b 的解

       A^(-1),       这是矩阵方程AX=E的解，E是单位矩阵。

    

      另外，LU分解之后还可以直接计算方阵的行列式。
 

3.  分块LU分解算法

        如果矩阵很大，采用分块计算能有效减小系统cache miss，这也是很多商业软件的实现方法。分块算法需要根据非分块算法本身重新设计算法流程，而不是简单在代码结构上用分块内存直接去改。线性代数的开源软件有很多，这里我就不枚举了。我主要测试了MATLAB和openCv的实现。MATLAB的矩阵运算的效率是及其高效的，openCv里面调用了著名的LAPACK。大概看了LAPACK的实现，用的也是分块算法。

       LU分解的分块算法的文献比较多，我主要参考了下面的两篇文献：

       LU分解分快算法的研究与实现

       LU分解递归算法的研究

       我作了两张图，可以详细的描述算法，这里以应用比较广泛的部分选主元LU块分解算法的执行过程。

 

        图中的画斜线的阴影部分，表示要把当前块LU分解得到的排列矩阵左乘以这部分数据组成的子矩阵，以实现行交换。从上图可以看出，在第一块分解之后，只需要按照排列矩阵交换A12，A22组成的子矩阵，而后面的每一次，则需要交换两个子矩阵。

         块LU分解算法主要由4部分构成：

         非块的任意瘦型矩阵的LU分解, 行交换，下三角矩阵方程求解, 矩阵乘法.

         LU分解来自方阵的三角分解。实际上，任意矩阵都有LU分解。但这里一般需要求解非分块的瘦型矩阵的LU分解，可以采用任意的部分选主元的LU分解算法。但是实现起来仍然有讲究，如果按照LAPACK实现的算法仍然不会快，而采用crout算法实现的结果是很快的。在我的测试中，采用crout算法的1024大小的矩阵非分块的LU分解和LAPACK实现的分块大小为64时的性能相当。LAPACK实现的算法本身是很高效的，但是其代码本身没有做太多的优化。实际上，没有经过任何优化的LAPACK的代码仍然比较慢。

  &nbs