高斯消元法应该没什么好说的，看维基百科上还分高斯消元法和高斯-若尔当消元法，不知道有什么区别。也没有仔细去看，该因该法实在简单，去扣细节没什么必要。所谓高斯消元法，其实就是我们在小学还是初中学的求解线性代数方程组的消元法。程序算法要做的就是把这个方法用计算机语言实现。
做法很简单，首先给出系数矩阵，然后通过矩阵行变换，用矩阵的第0行消去其他行的第0个元素，然后用第一行的第一个元素消去其他行的第一个元素，当然同时也要对等式右端做同样的变换。依次消除，最后会得到一个对角矩阵。对角矩阵和未知数列向量相乘后，各个未知数之间相互独立，和等式右端的列向量建立一一对应关系，就可以解出所有的未知数。
如果等式右端是一个单位矩阵的话，根据线性代数的理论，我们解出的未知数矩阵很明显就是系数矩阵的逆，所以这也是求解逆矩阵的一种方法。即给系数矩阵右边添加一个单位矩阵，变成一个增广矩阵，然后只进行行变换，当左边的部分变换成单位阵时，右边的部分就是原矩阵的逆了。
这里有个问题就是，在消元法中，当前的那个用来消去其他行中同列元素的“元”可能不够好，比如这个元是0，或者是一个非常小的数。这个时候如果我们选择这个元乘一个系数来消除其他行中同列元素的话，或者会导致算法失效，或者会使结果差的很大。大学期间上过计算物理这门课，当时学过高斯消元法。根据我的经验，当一个矩阵大到一百多个元素的话，基本上总是会出现（当然这可能跟我当时处理的问题有关）。总之，直接的高斯消元法不够安全，即使系数矩阵是非奇异的，也可能发生算法失效的情况，或者产生非常大的误差。本人当时脑洞一开，每次消元前都寻找要消元的列中最大的元素，然后进行行变换，用这个最大元素来消元，结果效果很好。这个办法本教材中同样采用，被叫做选择主元法（系数矩阵和右侧列向量同时做行交换，不会影响未知数列向量），看来本人还是很有先见之明嘛。 -_-
本书的选择主元法高斯消元函数有一点值得吐槽，它真的对矩阵进行了行交换，这样很没有必要。我们可以把每一行的首元素地址存在一个指针数组里，进行行交换的时候，只是交换指针数组里对应的两个指针就行，这样能节省很多没必要的复制。TNT的Array2D就是这么干的。
相应函数的代码：
[code title=”gaussj.h” lang=”cpp”]void gaussj(MatDoub_IO &a, MatDoub_IO &b)
{
Int i,icol,irow,j,k,l,ll,n=a.nrows(),m=b.ncols();
Doub big,dum,pivinv;
VecInt indxc(n),indxr(n),ipiv(n);
for (j=0;j<n;j++) ipiv[j]=0;
for (i=0;i<n;i++) {
big=0.0;
for (j=0;j<n;j++)
if (ipiv[j] != 1)
for (k=0;k<n;k++) {
if (ipiv[k] == 0) {
if (abs(a[j][k]) >= big) {
big=abs(a[j][k]);
irow=j;
icol=k;
}
}
}
++(ipiv[icol]);
if (irow != icol) {
for (l=0;l<n;l++) SWAP(a[irow][l],a[icol][l]);
for (l=0;l<m;l++) SWAP(b[irow][l],b[icol][l]);
}
indxr[i]=irow;
indxc[i]=icol;
if (a[icol][icol] == 0.0) throw("gaussj: Singular Matrix");
pivinv=1.0/a[icol][icol];
a[icol][icol]=1.0;
for (l=0;l<n;l++) a[icol][l] *= pivinv;
for (l=0;l<m;l++) b[icol][l] *= pivinv;
for (ll=0;ll<n;ll++)
if (ll != icol) {
dum=a[ll][icol];
a[ll][icol]=0.0;
for (l=0;l<n;l++) a[ll][l] -= a[icol][l]*dum;
for (l=0;l<m;l++) b[ll][l] -= b[icol][l]*dum;
}
}
for (l=n-1;l>=0;l–) {
if (indxr[l] != indxc[l])
for (k=0;k<n;k++)
SWAP(a[k][indxr[l]],a[k][indxc[l]]);
}
}
void gaussj(MatDoub_IO &a)
{
MatDoub b(a.nrows(),0);
gaussj(a,b);
}[/code]
Visits: 1514