一、高斯消元的原理
對於n元的m個線性方程組成的方程組，我們將其以矩陣的形式記錄下來：
a11 a12 a13 ...... a1n b1
a21 a22 a23 ...... a2n b2
...
...
...
an1 an2 an3 ...... ann bn
然后進行初等行列變換，嘗試構造出一個上三角矩陣，逐步使系數不為零的項減少；
等最后只剩下一個系數不為零時，進行回代，逐步求出已知解。（詳解過程咨詢小學老師）
二、高斯消元的實現
老老實實的回代代碼參見其他人的博客，這里介紹一種比較毒瘤的不回代暴力消元法：
1.Process
對於每個方程，按照一定的規則（后話）挑選一個主元，記錄該主元對應第幾個方程，然后用初等行列變換消去其他所有該元的系數；
最后我們得到的是一個每行只有一個數不為零，每列只有一個數不為零的鬼畜矩陣（自己腦補）
此時令ans向量對應的數字出去該行的非零系數，即為對應該元的解。
2.Code
設a數組為已經記錄系數的數組（格式見上方），那么a應該是n行n+1列的，最后一列代表方程的常數項；
設w數組記錄每個方程的主元是第幾項，v數組記錄答案；
當多解時輸出“Multiple solutions”，無解時輸出”No solution”；
double a[max_n][max_n+1],v[max_n];
int w[max_n];
void gauss(){
double eps=1e-6;
for(int i=1;i<=n;++i){ //Enumerate the equation;
int p=0; //Record the position of the largest number;
double mx=0; //Recording the largest number;
for(int j=1;j<=n;++j)
if(fabs(a[i][j])-eps>mx){
mx=fabs(a[i][j]);p=j; //fabs() returns the absolute value of float;
}
if(!p){
if(fabs(a[i][n+1]<eps))printf("Multiple solutions");
else printf("No solution");
return;
}
w[i]=p;
for(int j=1;j<=n;++j)
if(i!=j){ //other equations
double t=a[j][p]/a[i][p];
for(int k=1;k<=n+1;++k) //n+1 is important
a[j][k]-=a[i][k]*t;
}
}
for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];
}
3.notice
（1）精度的設置
眾所周知浮點數是有精度丟失的，在高斯消元中，精度的選擇要依題目而定，精度過低會導致系數較小的數被誤認為系數為零，而精度過高也有可能會導致誤差大於精度，導致本應該系數消為0的項誤認為系數不為零，所以精度的選擇是很哲學的問題，要注意。
推薦范圍：1e-4到1e-10
（2）主元的選取原則
接着（1）說，我們知道，用浮點數是有精度丟失的，那么讓一個較大的數除以一個極小的數誤差自然大的可想而知，所以我們想得到在精度允許的條件下系數最大的主元，所以對於每個方程，我們都應該選擇最大系數的元作為主元。
（3）在模2意義下的高斯消元
使用bitset優化運行時間，詳見相關應用中第三個例題的講解；
三、相關應用
這里給出高斯消元的幾道基礎題目，難度適合初學者。
1.[Luogu P3389]【模板】高斯消元
Description
給定一個線性方程組，對其求解
輸入格式：
第一行，一個正整數 n
第二至 n+1行，每行 n+1個整數，為 a1,a2⋯an和 b，代表一組方程。
輸出格式：
共n行，每行一個數，第 i行為 xi(保留2位小數）
如果不存在唯一解，在第一行輸出"No Solution".
Solution
如上所述。
Code
#include<iostream>
#include<cstdio>
#include<cmath>
#include<cstring>
#include<algorithm>
using namespace std;
const int max_n=110;
double a[max_n][max_n+1],v[max_n];
int n,w[max_n];
inline int rd(){
int x=0;
bool f=0;
char c=getchar();
while(!isdigit(c)){
if(c=='-') f=1;
c=getchar();
}
while(isdigit(c)){
x=(x<<1)+(x<<3)+(c^48);
c=getchar();
}
return f?-x:x;
}
void gauss(){
double eps=1e-6;
for(int i=1;i<=n;++i){//enumerate the equation;
int p=0; //Record the position of the largest number;
double mx=0; //Recording the largest number;
for(int j=1;j<=n;++j)
if(fabs(a[i][j])-eps>mx){
mx=fabs(a[i][j]);p=j;//fabs() returns the absolute value of float;
}
if(!p){
printf("No Solution");
return;
}
w[i]=p;
for(int j=1;j<=n;++j)
if(i!=j){ //other equations
double t=a[j][p]/a[i][p];
for(int k=1;k<=n+1;++k)//n+1 is important
a[j][k]-=a[i][k]*t;
}
}
for(int i=1;i<=n;++i) v[w[i]]=a[i][n+1]/a[i][w[i]];
for(int i=1;i<=n;++i) printf("%.2lf\n",v[i]);
}
int main(){
n=rd();
for(int i=1;i<=n;++i)
for(int j=1;j<=n+1;++j)
a[i][j]=rd();
gauss();
return 0;
}
2.[BZOJ 1013][JSOI 2008] 球形空間產生器sphere
詳解參考我的隨筆：http://www.cnblogs.com/COLIN-LIGHTNING/p/8982341.html