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 }