算法学习之一：矩阵计算——高斯消元法

在许多应用中，我们需求解出一个有n个方程n个未知数的方程组

a11 x1 + a12x2 …… +a1nxn =
b1
a21x1 +a22x2 …… +a2nxn =
b2
.
.
.
.
an1 x1 +an2x2 …… +annxn =
bn

那么，学过线性代数的就都知道只要解出 Ax=b这个非齐次的解即可。
其所用的方法就是高斯消元法大致的思想便是由固定一行，由此行逐行向下运算，将固定行的第一个非零数据列在该行以下数据全部消去为0。这样执行至倒数第二行时便可以将原方阵变换为一个上三角方阵，如果将此方阵变为其增广矩阵在进行运算的话便可以根据变换后的矩阵求出这个方程组的n个解。
在此只先考虑该系数矩阵是一个满秩矩阵的情况。
对于一个矩阵首先要为其开辟存放该矩阵的空间，由于变换时是操作其增广矩阵，那么不妨开一个 n * （n + 1）大小的矩阵，并且需要一个一维矩阵来存放计算出的n个未知数，这些空间在使用完之后要将其释放掉，实现代码如下：

double **mat = NULL, *x = NULL;

void applyRoom(int n){
    int i = 0;
    x = (double *)malloc(sizeof(double) * n);

    mat = (double **)malloc(sizeof(double *) * n);
    for(; i < n; i++)
        mat[i] = (double *)malloc(sizeof(double) * (n + 1));
}

void freeRoom(int n){
    int i = 0;
    for(; i < n; i++)
        free(mat[i]);
    free(mat);
    free(x);
}

存放实现时候接下来就是输入的实现了，避免麻烦。此处将整个增广矩阵一起输入：

void inputMat(int n){
    int i = 0, j = 0;
    for(i = 0; i < n; i++)
        for(j = 0; j < n + 1; j++)
            scanf("%lf", &mat[i][j]);
}

在具体实现这个增广变换的时候，需要进行行之间的交换，因为如果作为固定行是一个很小的浮点数，本身在机器中浮点数的存放就是不精确的，那么别的数再除以它就会将这个误差放大，最后可能会得到我们所不想要的结果，那么再每次固定行之前，我们就要向下进行一次比较，找到开头数据是最大的那行，将其交换上来。所以，我们需要一个实现交换的函数：

void swapLine(int a, int b, int n){
    double mid = 0.0;
    int i = 0;
    for(; i < n + 1; i++){
        mid = mat[a][i];
        mat[a][i] = mat[b][i];
        mat[b][i] = mid;
    }
}

当固定行确定之后，接下来我们就要开始执行消去了。我们需要记录是哪两行进行操作，消去操作是从那一列开始（虽然开始列的列数跟固定行的行数是相同，但还是找一个参数来传递吧，这样看的清楚一些）。这样我们就能在两行之间执行操作了：

void binOperat(const int a, const int b, int i, const int n){
    double temp = mat[b][i] / mat[a][i];
    for(; i < n + 1; i++)
        mat[b][i] -= mat[a][i] * temp;
}

当我们讲系数矩阵变成一个上三角矩阵后，接下来就是对这个变换过后的增广矩阵从最后一行开始，逐行向上操作进而计算出每一个未知数的值啦：

void calculateX(int n){
    int m = n, i = 0, j = 0;
    double mid = 0.0;
    while(m--){
        mid = mat[m]
;
        for(i = n - 1; i > m; i--)
            mid -= x[i] * mat[m][i];
        x[m] = mid / mat[m][m];
    }
}

每一个小的模块都可以确定了，可开心了呢~将这些函数组合起来就可以得到一个基于高斯消元法的解出n个方程的n个未知数的小程序了。

int main(){

    int n = 0, max = 0;
    int i = 0, j = 0;
    scanf("%d", &n);

    applyRoom(n);
    inputMat(n);

    for(i = 0; i < n; i++){
        max = i;
        for(j = i + 1; j < n; j++)
            max = fabs(mat[j][i]) > fabs(mat[max][i]) ? j : max;
        swapLine(max, i, n);

        for(j = i + 1; j < n; j++)
            binOperat(i, j, i, n);
    }

    calculateX(n);

    printf("%.6lf", x[0]);
    for(i = 1; i < n; i++)
        printf(" %lf", x[i]);
    printf("\n");

    freeRoom(n);

    return 0;
}

由这个方法变换一下便可以得到求行列式的值得方法了，比起刚学线性代数辛辛苦苦异常头疼的求解行列式的值，这样做是不是觉得高端了不少呢~很酷有没有~哈哈。。。

PS：还是个学生，刚开始学习算法课，这个程序的正确性还不太确定，自己尝试了几个方程组，还都是对的。求各位大神拍砖！我要在跌跌撞撞中进步！谢！！！

其中 是一个大整数.理论上, 我们可通过推广上述求解系统(*1)的迭代法来求解系统(*2)
; 然而