高斯消元

 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
");
13     for(i = 0 ; i < n ; i++){
14         // 选一行r并与第i行交换
15         //printf("do
");
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
",A[i][3]);
48     return 0;
49 }
原文地址:https://www.cnblogs.com/cyb123456/p/5824244.html