1 #include<iostream>
2 #include<cstdio>
3 #include<iomanip>
4 using namespace std;
5 #define e 0.00000001
6 #define maxn 50
7
8 int n;//规模nXn
9 double a[maxn][maxn];//系数矩阵
10 double b[maxn];//b矩阵
11 double m[maxn][maxn];//中间变量矩阵
12 double x[maxn];//最终解
13 int H=1;//扩大H被结算（优化）
14 /*
15 读取数据
16 */
17 void read(){
18 cout<<"请输入系数矩阵规模n:= ";
19 cin>>n;
20 cout<<"|-----------------------------\n";
21 cout<<"|请输入系数矩阵,如：\n";
22 cout<<"|1.1348 3.8326 1.1651 3.4017\n";
23 cout<<"|0.5301 1.7875 2.5330 1.5435\n";
24 cout<<"|3.4129 4.9317 8.7643 1.3142\n";
25 cout<<"|1.2371 4.9998 10.6721 0.0147\n";
26 cout<<"|-----------------------------\n";
27 for(int i=1;i<=n;i++)
28 for(int j=1;j<=n;j++){
29 cin>>a[i][j];
30 a[i][j]*=H;
31 }
32 cout<<"|-----------------------------\n";
33 cout<<"|请输入b矩阵，如：\n";
34 cout<<"|9.5342 6.3941 18.4231 16.9237\n";
35 cout<<"|-----------------------------\n";
36 for(int i=1;i<=n;i++){
37 cin>>b[i];
38 b[i]*=H;
39 }
40 }
41
42 /*
43 中间矩阵输出
44 参数：消元次数
45 */
46 void PrintProc(int cases){
47 printf("--------第%d次消元结果如下：\n",cases);
48 for(int i=1;i<=n;i++){
49 for(int j=1;j<=n;j++){
50 cout<<setw(10)<<a[i][j]<<' ';
51 }
52 cout<<setw(10)<<b[i]<<'\n';
53 }
54 cout<<"END THIS SHOW-------------\n";
55 }
56
57 /*
58 显示结果
59 */
60 void Print(){
61 cout<<"|-----------------------------\n";
62 cout<<"|结果为：\n";
63 for(int i=1;i<=n;i++){
64 printf("x[%d]= %lf\n",i,x[i]);
65 }
66 cout<<"|-----------------------------\n\n";
67 }
68
69 /*
70 顺序消元法
71 */
72 void ShunXuXiaoYuan(){
73 //消元计算
74 for(int k=1;k<n;k++){
75 for(int i=k+1;i<=n;i++){
76 m[i][k]=a[i][k]/a[k][k];
77 for(int j=k+1;j<=n;j++){
78 a[i][j]-=m[i][k]*a[k][j];
79 }
80 }
81 for(int i=k+1;i<=n;i++){
82 b[i]-=m[i][k]*b[k];
83 }
84 PrintProc(k);//输出中间计算过程
85 }
86 //回代求解
87 x[n]=b[n]/a[n][n];
88 for(int i=n-1;i>0;i--){
89 x[i]=b[i];
90 for(int j=i+1;j<=n;j++)
91 x[i]-=a[i][j]*x[j];
92 x[i]/=a[i][i];
93 }
94 //输出结果
95 Print();
96 }
97
98 /*
99 列主消元
100 */
101 void LieZhuXiaoYuan(){
102 for(int k=1;k<n;k++){
103 //选主元[这一列的绝对值最大值]
104 double ab_max=-1;
105 int max_ik;
106 for(int i=k;i<=n;i++){
107 if(abs(a[i][k])>ab_max){
108 ab_max=abs(a[i][k]);
109 max_ik=i;
110 }
111 }
112 //交换行处理[先判断是否为0矩阵]
113 if(ab_max<e){//0矩阵情况
114 cout<<"det A=0\n";
115 break;
116 }else if(max_ik!=k){//是否是当前行，不是交换
117 double temp;
118 for(int j=1;j<=n;j++){
119 temp=a[max_ik][j];
120 a[max_ik][j]=a[k][j];
121 a[k][j]=temp;
122 }
123 temp=b[max_ik];
124 b[max_ik]=b[k];
125 b[k]=temp;
126 }
127 //消元计算
128 for(int i=k+1;i<=n;i++){
129 a[i][k]/=a[k][k];
130 for(int j=k+1;j<=n;j++){
131 a[i][j]-=a[i][k]*a[k][j];
132 }
133 b[i]-=a[i][k]*b[k];
134 }
135 PrintProc(k);//输出中间计算过程
136 if(k<n-1)continue;
137 else{
138 if(abs(a[n][n])<e){
139 cout<<"det A=0\n";
140 break;
141 }else{//回代求解
142 x[n]=b[n]/a[n][n];
143 for(int i=n-1;i>0;i--){
144 x[i]=b[i];
145 for(int j=i+1;j<=n;j++)
146 x[i]-=a[i][j]*x[j];
147 x[i]/=a[i][i];
148 }
149 //输出结果
150 Print();
151 }
152 }
153 }
154 }
155
156 /*
157 主函数
158 */
159 int main(){
160 while(1){
161 read();
162 LieZhuXiaoYuan();
163 //ShunXuXiaoYuan();
164 }return 0;
165 }
166 /*
167 书上高斯消元的例子：
168 1 1 1
169 1 3 -2
170 2 -2 1
171
172 6 1 1
173 */
174 /*
175 书上列主消元的例子：
176 -0.002 2 2
177 1 0.78125 0
178 3.996 5.5625 4
179
180 0.4 1.3816 7.4178
181 */
