高斯消元其實在算法競賽中算是一個十分常見的算法。它的大致思想就和初中階段學(xué)到的加減消元法差不多。這個算法的時間復(fù)雜度為\(O(n^3)\)，是一個相當簡單的算法，但是具體實現(xiàn)需要一些思考。

這里給出模板題的鏈接：
洛谷P3389
P4035

1.1 問題引入

給定方程組\(\begin{cases}x+3y+4z=5\quad(1)\\x+4y+7z=3\quad(2)\\9x+3y+2z=2\quad(3)\end{cases}\),求解之。

這是一個相當簡單的三元一次方程組，直接用加減消元就可以得出解。當然，這里給出一個比較不討巧的消元順序，為了方便后面的理解。

首先（2）式減去（1）式：\(y+3z = -2\quad(2)'\)

接下來\(9\times(1)-(3)\):\(24y+34z=43\quad(3)'\)

此時\(y\)和\(z\)就構(gòu)成了二元一次方程組了。

計算\(24\times(2)'-(3)'\): \(38z=-91\).

解得\(z=-\frac{91}{38}\).此時不斷的回代可以得到\(y=\frac{197}{38}\), \(x=-\frac{37}{38}\).

我們可以把整個計算過程得到的新的方程式列出來：

仔細觀察一下這個方程組，這對于之后的解題有幫助。

思考一下，對于任意的三元一次方程組，或者更進一步的，對于任意的\(n\)元一次方程組，我們能不能設(shè)計一個通用的算法，求解方程組呢？

1.2 初等列變換

我們把\(x\),\(y\),\(z\)的系數(shù)提取出來，將它們列成一個3*3的表：\(\begin{bmatrix}1&3&4\\1&4&7\\9&3&2\end{bmatrix}\)

我們還可以把方程的右側(cè)常數(shù)列在右邊：\(\begin{bmatrix}1&3&4&|5\\1&4&7&|3\\9&3&2&|2\end{bmatrix}\)
這個3*4的表格叫做方程組的增廣矩陣。規(guī)定矩陣的第\(i\)行為\(r_i\),如\(r_1=\begin{bmatrix}1&3&4&|5\end{bmatrix}\).
我們可以對它進行這樣幾個操作：

• 用一個非零常數(shù)乘上某一行。這個操作記為\(r_i=cr_i\)，其中\(c\)為非零常數(shù)
• 把其中一行的若干倍加至另一行上。這個操作記為\(r_i=cr_i+c'r_j\)
• 交換兩行的位置。直接記作\(\text{swap}(r_i,r_j)\)就好了。

這樣一來，我們之前的操作可以這樣表示：
1.\(r_2=r_2-r_1\)
2.\(r_3=9r_1-r_3\)
3.\(r_3=24r_2-r_3\)

這個時候增廣矩陣變成了這個樣子：\(\begin{bmatrix}1&3&4&|5\\0&1&3&|-2\\0&0&38&|-91\end{bmatrix}\)

往往矩陣元素等于0時我們省略不寫，上面的矩陣還可以這樣寫：\(\begin{bmatrix}1&3&4&|5\\ &1&3&|-2\\ & &38&|-91\end{bmatrix}\)

排版有點丑，但是我們還是可以看出來，這個矩陣左邊的形狀有點像一個倒過來的階梯。人們就叫它簡化階梯型矩陣。把原增廣矩陣變?yōu)楝F(xiàn)在的簡化階梯型矩陣的過程就叫做高斯消元。
算出了簡化階梯型矩陣之后，直接一步一步回代就可以了。

1.3 基本實現(xiàn)

總結(jié)一下，高斯消元的思路：
對于每一行\(r_i\)，我們要讓\(r_{i,i}=1\)，并且\(r_{i,i}\)下面的元素均變?yōu)?span id="w0obha2h00" class="math inline">\(0\)。我們這樣做：

1. 選取\(|r_{j,i}|\)最大的一行\(r_j\)，將它與\(r_i\)交換。之后\(r_{i,i}\)就會變?yōu)楫斍傲凶畲蟮脑亍?/li>
2. 把\(r_{i,i}\)變?yōu)?span id="w0obha2h00" class="math inline">\(1\)。很簡單，直接讓\(r_i=\frac{1}{r_{i,i}}r_i\)就行了。
3. 把\(r_{i,i}\)下面的元素都變?yōu)榱恪τ?span id="w0obha2h00" class="math inline">\(i+1 \leq j \leq n\)，我們讓\(r_j=r_j-\frac{r_{j,i}}{r_{i,i}}r_i\)，調(diào)整\(r_i\)的比例并把\(r_{j,i}\)消去。因為\(r_{i,i}=1\)，這個操作可以寫成\(r_j=r_j-r_{j,i}r_i\)就行了。

這一段的大概代碼如下：

	RP(i,1,n)
	{
		rg int cur=i;
		RP(j,1,n)
		if(fabs(r[cur][i])<fabs(r[j][i]))
			cur=j;
		if(i!=cur)
			swap(r[cur],r[i]);
		db div=r[i][i];
		RP(j,i,n+1)
			r[i][j]/=div;
		RP(row,i+1,n)
		{
			db rat=r[row][i];
			RP(col,i,n+1)
				r[row][col]-=rat*r[i][col];
		}
	}

在最后回代的過程中，由于最后一行的形式一定是一個一元一次方程，最后一行的解就是\(x_n=r_{n+1}\)。
之后的每一行\(r_i\):\(\sum\limits_{j=i}^n{x_j\cdot r_{i,j}}=r_{i,n+1}\)，移項就得到了解：$$x_i=\dfrac{r_{i,n+1}-\sum\limits_{j=i+1}^n{x_j\cdot r_{i,j}}}{r_{i,i}}=r_{i,n+1}-\sum\limits_{j=i+1}^n{x_j\cdot r_{i,j}}$$

x[n]=r[n][n+1];
    DRP(i,n-1,1)
    {
        x[i]=r[i][n+1];
        RP(j,i+1,n)
            x[i]-=r[i][j]*x[j];
    }

1.4 一些特判細節(jié)和其他注意事項

如果方程的某一列元素均為\(0\)，而常數(shù)項\(\neq 0\)，那么方程組是一定無解的。這里直接特判就可以了。
另外，在交換行時可以特判一下兩行是否相等。減掉沒必要的操作可以稍微提速。
這個算法還有一個加強版本:高斯-約旦消元法。這個算法其實就是把回代的操作直接用初等行變換實現(xiàn)，不過這個算法在矩陣求逆上有著很好的應(yīng)用。