高斯消元

#include <cstdio>
#include <cmath>
#include <iostream>
const int maxn = 100;
typedef double Matrix[maxn][maxn];
// 要求系数矩阵可逆
// 这里A是增广矩阵，A[i][n]表示第i个方程右边的常数bi
// 运行结束后A[i][n]是第i个未知数的值
void gauss_elimination(Matrix A,int n){
    int i,j,k,r;
    // 消元过程
    //printf("do
");
    for(i = 0 ; i < n ; i++){
        // 选一行r并与第i行交换
        //printf("do
");
        r = i;
        for(j = i + 1 ; j < n ; j++)
            if(fabs(A[j][i]) > fabs(A[r][i])) r = j;
        if(r != i) for(j = 0 ; j <= n ; j++) std::swap(A[r][j],A[i][j]);

        // 与第i+1~n行进行消元
        for(k = i + 1 ; k < n ; k++){
            double f = A[k][i] / A[i][i];
            for(j = i ; j <= n ; j++) A[k][j] -= f * A[i][j];
        }
        /* // 如果要求高精度
        for(j = n ; j >= i ; j--) // 必须逆序枚举
            for(k = i + 1 ; k < n ; k++)
                A[k][j] -= A[k][i] / A[i][i] * A[i][j];
        */
    }
    // 回带过程
    for(i = n - 1 ; i >= 0 ; i--){
        for(j = i + 1 ; j < n ; j++)
            A[i][n] -= A[j][n] * A[i][j];
        A[i][n] /= A[i][i];
    }
}
int main(){
    Matrix A;
    A[0][0] = 2,A[0][1] = 1,A[0][2] = -1;
    A[1][0] = -3,A[1][1] = -1,A[1][2] = 2;
    A[2][0] = -2,A[2][1] = 1,A[2][2] = 2;
    A[0][3] = 8,A[1][3] = -11,A[2][3] = -3;
    gauss_elimination(A,3);
    for(int i = 0 ; i < 3 ; i++)
        printf("%lf
",A[i][3]);
    return 0;
}