高斯-约旦消元法中的一个细节

首先给出几个定义:

非零行:矩阵至少有一个非零元素的行
矩阵中非零行的先导元素指的是该行最左边的非零元素（感觉这个名称有点多余）
一个矩阵为(行)阶梯型((row) echelon form),若它有以下三个性质:

1. 每一非零行都在每一零行的上面
2. 某一行的最左边的非零元素位于前一行最左边的非零元素的右侧
3. 某一最左边的非零元素所在列下方元素都是0

简化阶梯型矩阵还需满足
每一非零行的先导元素为1,先导元素是该列唯一非零元素

阶梯型矩阵中的先导元素位置被称为主元位置，(对于任意矩阵)主元位置上的非零元素称为主元

\[\begin{bmatrix} 2 \ *\ *\ * \\0\ 3\ *\ * \\ 0\ 0\ 0\ 0\ \end{bmatrix} \]

注意阶梯型矩阵第\(i\)行的主元不一定在第\(i\)列,而是在该行第一个非零的位置，即所谓的“阶梯型”不一定是连续的

这是一个阶梯型矩阵。

这是OI中常使用的高斯消元板子。该板子可以在有\(n\)个方程\(n\)个未知数，保证方程有唯一解或无解的情况下使用
有唯一解的时候,消元完后增广矩阵第\(i\)行第\(i\)列一定为1,而第\(i\)行第\(n+1\)列则是\(x_i\)的解

//LuoguP3389 n个方程,n个未知数
bool gaussJordan(int n){
	for(int i=1;i<=n;i++){//枚举主元(共n个未知数)
		int id=i;
		for(int j=i+1;j<=n;j++) if(fabs(a[j][i])>fabs(a[id][i])) id=j;   
		for(int j=1;j<=n+1;j++) swap(a[i][j],a[id][j]);
		if(a[i][i]==0) return 0;
		for(int j=1;j<=n;j++){//从1开始,把前面的都消掉,这样最后只剩对角线 
			if(j==i) continue;//第i行不动 
			double d=a[j][i]/a[i][i];
			for(int k=i;k<=n+1;k++) a[j][k]=(a[j][k]-a[i][k]*d);
		}
	}
	for(int i=1;i<=n;i++)a[i][n+1]/=a[i][i];
	
	return 1; 
} 

但如果不是\(n\)个方程\(n\)个未知数，方程可能有无数组解时，这种写法就会出问题。
考虑一个简单的方程组 \(\begin{cases}x+y+z=5 \\ x+y=0 \end{cases}\)

它对应的增广矩阵为

\[\begin{bmatrix} 1 \ 1 \ 1 \ 5 \\1 \ 1 \ 0 \ 0 \end{bmatrix} \]

首先处理第1列，第一列把主元下方的位置变成0, 那么把第二行减去第一行，增广矩阵变成了

\[\begin{bmatrix} 1 \ 1 \ 1 \ 5 \\0 \ 0 \ -1 \ -5 \end{bmatrix} \]

上述板子里的高斯消元法默认增广矩阵第\(i\)行第\(i\)列一定不为0，如果为0，则判定方程无解。板子会把无数解判定为无解，这必须要修改。因此，对于第2行，我们应保持操作行不变，继续处理第3列，如果还有更多方程的话，应该把第3列主元下方的位置变成0。这样就可以维持阶梯形的性质。

修改后的代码如下

inline bool isZero(vtype v){return fabs(v)<EPS;}
int gauss(int varN,int eqN){//varN是变量数,eqN是方程数
	int nwl=1;//当前处理的方程(行) 
	for(int i=1;i<=varN;i++){//i是当前的主元序号(列)，注意把i和nwl区分开来
		int tmp=nwl;//选取一列中绝对值最大的作为主元 
		for(int j=i+1;j<=eqN;j++) if(fabs(mat[j][i])>fabs(mat[tmp][i])) tmp=j;
		if(isZero(mat[tmp][i])) continue;//找不到主元，转而处理下一列, 操作行不变 
		for(int j=1;j<=varN+1;j++) swap(mat[tmp][j],mat[nwl][j]);
		for(int j=1;j<=eqN;j++){
			if(j==nwl) continue;
			double d=mat[j][i]/mat[nwl][i];//第nwl行,把第i列除nwl外都变成0
			for(int k=1;k<=varN+1;k++) mat[j][k]-=mat[nwl][k]*d;
			//这里令j=1,k=1可以真正消出简化阶梯型矩阵。如果有唯一解:i+1之前我们知道一定会被消成0,就不用计算了。但如果要处理其他的情况,就要从1开始 
		}
		nwl++; 
	}//nwl之后的行都是0 
	return nwl; 
}

结束消元后，如果不是无解，对nwl之前的行的每个主元都可以给出它的表达式（可用自由变量表示）

结果判定的方法如下

if(nwl<=varN){//此时一定存在varN-nwl+1个自由变量或无解(加1是因为循环里多++了一次)
		for(int i=nwl;i<=varN;i++){//依次判断每个自由变量
			if(!isZero(mat[i][varN+1])){
				printf("方程组无解");
				return 0;
			}//如果都是0 0就不用管
		}
		printf("方程组有无数解\n");
	}else{//有唯一解
		for(int i=1;i<=varN;i++){
			cout<<vars[i]<<"=";
			printf("%.3lf",mat[i][varN+1]/mat[i][i]); 
			cout<<"\n";
		}
	} 

重点是:第\(i\)行的主元不一定在第\(i\)列,而是在该行第一个非零的位置，即所谓的“阶梯型”不一定是连续的

LuoguP2455 [SDOI2006]线性方程组就考察了这种情况

手写的方程组求解器