1 #include <cstdio> 2 #include <cmath> 3 #include <iostream> 4 const int maxn = 100; 5 typedef double Matrix[maxn][maxn]; 6 // 要求系数矩阵可逆 7 // 这里A是增广矩阵，A[i][n]表示第i个方程右边的常数bi 8 // 运行结束后A[i][n]是第i个未知数的值 9 void gauss_elimination(Matrix A,int n){10     int i,j,k,r;11     // 消元过程12     //printf("do\n");13     for(i = 0 ; i < n ; i++){14         // 选一行r并与第i行交换15         //printf("do\n");16         r = i;17         for(j = i + 1 ; j < n ; j++)18             if(fabs(A[j][i]) > fabs(A[r][i])) r = j;19         if(r != i) for(j = 0 ; j <= n ; j++) std::swap(A[r][j],A[i][j]);20 21         // 与第i+1~n行进行消元22         for(k = i + 1 ; k < n ; k++){23             double f = A[k][i] / A[i][i];24             for(j = i ; j <= n ; j++) A[k][j] -= f * A[i][j];25         }26         /* // 如果要求高精度27         for(j = n ; j >= i ; j--) // 必须逆序枚举28             for(k = i + 1 ; k < n ; k++)29                 A[k][j] -= A[k][i] / A[i][i] * A[i][j];30         */31     }32     // 回带过程33     for(i = n - 1 ; i >= 0 ; i--){34         for(j = i + 1 ; j < n ; j++)35             A[i][n] -= A[j][n] * A[i][j];36         A[i][n] /= A[i][i];37     }38 }39 int main(){40     Matrix A;41     A[0][0] = 2,A[0][1] = 1,A[0][2] = -1;42     A[1][0] = -3,A[1][1] = -1,A[1][2] = 2;43     A[2][0] = -2,A[2][1] = 1,A[2][2] = 2;44     A[0][3] = 8,A[1][3] = -11,A[2][3] = -3;45     gauss_elimination(A,3);46     for(int i = 0 ; i < 3 ; i++)47         printf("%lf\n",A[i][3]);48     return 0;49 }