高斯消元

Posted 2020-08-02 脑洞的分析与证明

 1 #include<iostream>
 2 #include<math.h>
 3 using namespace std;
 4 #define MAXN 100
 5 #define fabs(x) ((x)>0?(x):-(x))
 6 #define eps 1e-10
 7 //列主元gauss消去求解a[][]x[]=b[]
 8 //返回是否有唯一解,若有解在b[]中
 9 int gauss_cpivot(int m,int n, double **a, double b[])
10 //m,n为系数矩阵的行列，a[][]系数矩阵
11 //b[]常数列向量，将解x[]存入b[]中
12 {
13     int i, j, k, row;
14     double maxp, t;
15     for (k = 0; k<n; k++) {
16         for (maxp = 0, i = k; i<m; i++)
17             if (fabs(a[i][k])>fabs(maxp))
18                 maxp = a[row = i][k];
19         if (fabs(maxp)<eps)
20             return 0;
21         if (row != k) {
22             for (j = k; j<n; j++)
23                 t = a[k][j], a[k][j] = a[row][j], a[row][j] = t;
24             t = b[k], b[k] = b[row], b[row] = t;
25         }
26         for (j = k + 1; j<n; j++) {
27             a[k][j] /= maxp;
28             for (i = k + 1; i<m; i++)
29                 a[i][j] -= a[i][k] * a[k][j];
30         }
31         b[k] /= maxp;
32         for (i = k + 1; i<m; i++)
33             b[i] -= b[k] * a[i][k];
34     }
35     for (i = m - 1; i >= 0; i--)
36         for (j = i + 1; j<n; j++)
37             b[i] -= a[i][j] * b[j];
38     return 1;
39 }
40 void main()
41 {
42     int r, c,i,j;
43     double **a, *b;
44     cin >> r >> c;
45     a = new double*[r];
46     for (i = 0; i < r; i++)
47         a[i] = new double[c];
48     b = new double[r];
49     for (i = 0; i < r; i++)
50     {
51         b[i] = rand();
52         for (j = 0; j < c; j++)
53             a[i][j] = rand();
54     }
55     gauss_cpivot(r,c, a, b);
56     for (i = 0; i < r; i++)
57         cout << b[i] << endl;
58 }