矩阵的LU分解求解线性方程组(C++实现）

继续解方程组！
高斯消元法：http://blog.csdn.net/qq_26025363/article/details/53027926
列主元素高斯消元法：http://blog.csdn.net/qq_26025363/article/details/53044843，
这次使用LU分解求解方程组的解，该方法思想就是将一个矩阵分解为一个单位下三角矩阵L和一个上三角矩阵U，属于矩阵的三角分解法，又称杜利特尔（Doolittle）分解。其实高斯消元法的进行的每一步消元操作，都是等于左乘一个初等置换阵。具体的该方法的证明在数值分析或者线性代数的课本上应该都有，在这就不给出证明和说明了。本篇文章主要给出具体的C++实现过程。
还是老样子，先贴代码吧！

#include<iostream>
using namespace std;

const int n = 3;
//矩阵的ALU分解
void ALU(double a[n][n], double b[n])
{
    double l[n][n] = { 0 };
    double u[n][n] = { 0 };
    int i, r, k;
    //进行U的第一行的赋值
    for (i = 0; i<n; i++)
    {
        u[0][i] = a[0][i];
    }

    //进行L的第一列的赋值
    for (i = 1; i<n; i++)
    {
        l[i][0] = a[i][0] / u[0][0];
    }

    //计算U的剩下的行数和L的剩下的列数
    for (r = 1; r<n; r++)
    {
        for (i = r; i <n; i++)
        {
            double sum1 = 0;
            for (k = 0; k < r; k++)
            {
                sum1 += l[r][k] * u[k][i];
                //cout << "" << r << "" << sum1 << endl;
            }
            u[r][i] = a[r][i] - sum1;
        }


        if(r!=n)
        for(i=r+1;i<n;i++)
        {
            double sum2 = 0;
              for (k = 0; k<r; k++)
            {
                  sum2 += l[i][k] * u[k][r];
            }
                l[i][r] = (a[i][r] - sum2) / u[r][r];
        }

    }

    double y[n] = { 0 };
    y[0] = b[0];
    for (i = 1; i<n; i++)
    {
        double sum3 = 0;
        for (k = 0; k<i; k++)
            sum3 += l[i][k] * y[k];
        y[i] = b[i] - sum3;
    }

    double x[n] = { 0 };
    x[n - 1] = y[n - 1] / u[n - 1][n - 1];
    for (i = n - 2; i >= 0; i--)
    {
        double sum4 = 0;
        for (k = i + 1; k<n; k++)
            sum4 += u[i][k] * x[k];
        x[i] = (y[i] - sum4) / u[i][i];
    }
    for (i = 0; i<n; i++)
        cout << "x[" << i + 1 << "]=" << x[i] << endl;
    return;
}


int main()
{
    double a[3][3] = { 1,2,3,2,5,2,3,1,5 };
    double b[3] = { 14,18,20 };
    ALU(a, b);
    return 0;
}

来分析一下实现过程吧，首先进行A=LU分解，按照分解的方法来，U的第一行和L的第一列算是初始化，可以直接用循环简单计算赋值，接下来计算剩下的行和列。这里有个计算的次序问题，必须先算U的一行后，才能计算L的一列，不然是错误的计算方法，我一开始实现的时候就犯了这个大错，挑了好久的bug，还有一个注意点就是循环变量的值域问题，也是比较难以掌握的。LU分解完后，就可以进行直接的求解了，先按LY=B，求解Y，然后按UX=Y，求解X。具体就是上述这样。即求AX=B;先将A=LU,然后LUX=B,先LY=B,求出Y，然后UX=Y，求出X。
求解三角矩阵的时间复杂度是O(N^2)，但是分解过程是O(N^3)，但是总体来说，时间复杂度跟高斯消元差不多，但是避免了高斯消元法的主元素为0的过程。
空间性能分析：我的实现方法不是最省空间的，我每个中间矩阵都开辟了新的空间，无疑浪费了很多不用的空间。给出改进方法，A=LU分解的元素都可以存放在A的矩阵中，因为当用过了aij那个元素后，便不在用了，所以可以占用原来空间。但是各有利弊，考虑如果是上千个元素的矩阵，引用传参，这样
无疑就改变原矩阵了，如果不想改变原矩阵的话，还是老老实实建一个矩阵，然后把LU都放进去，我的是建了2个矩阵，极大的浪费。但是实现简单，逻辑也简单，所以就贴上了。
如果文章有错误或不妥之处，请各位指出，谢谢！