我們從一個高斯消元開始：

struct row{
	double a[105];
	row operator -(const row& other)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]-other.a[i];
		return b;
	}
	row operator *(const double& div)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]*div;
		return b;
	}
	row operator /(const double& div)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]/div;
		return b;
	}
	bool iszero(){
		bool ret=1;
		for(int i=1;i<m;i++) if(a[i]) ret=0;
		return ret;
	}
}r[105];
int Gauss_Elimination() 
{
	for(int i=1;i<m;i++){
		if(!r[i].a[i]){
			int f=0;
			for(int j=i+1;j<=n;j++)
				if(r[j].a[i])
					f=j;
			if(f) swap(r[i],r[f]);
		}
		if(r[i].iszero()){
			for(int j=i;j<=n;j++)
				if(r[j].a[m]!=0) return NO_SOLUTION;
			return UNLIMITED_SOLUTIONS;
		}
		int fir=0;
		for(int j=1;j<=m;j++)
			if(fabs(r[i].a[j])>1e-9) {fir=j;break;}
		for(int j=i+1;j<=n;j++) r[j]=r[j]-r[i]*(r[j].a[fir]/r[i].a[fir]);
		r[i]=r[i]/(r[i].a[fir]);
	}
	for(int i=m-1;i>=1;i--) 
		for(int j=m;j>i;j--)
			r[i]=r[i]-r[j]*(r[i].a[j]);
	return SOLVED;
}

值得注意的是，我發現很多人的高斯消元板子比我短很多，仔細觀看后發現主要有兩原因：
1.我是完全按照傳統線代教材：先化行階梯矩陣，最后最簡化；但是實際可以優化成直接最簡；
2.實際沒有必要判斷該行是不是頭非零，如果是零去找人更換，直接從這行開始找，非零更換即可（等價）
3.實際可以按列遍歷，這列沒有行就看下一列，此時行還是這一行。如果找到行，才在下一列換下一行。（這個倒沒太大所謂）

于是我們最后的模板為

struct row{
	double a[105];
	row operator -(const row& other)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]-other.a[i];
		return b;
	}
	row operator *(const double& div)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]*div;
		return b;
	}
	row operator /(const double& div)const{
		row b;
		for(int i=1;i<=m;i++) b.a[i]=a[i]/div;
		return b;
	}
	bool iszero(){
		bool ret=1;
		for(int i=1;i<m;i++) if(a[i]) ret=0;
		return ret;
	} 
}r[105];
int Gauss_Elimination() 
{
	for(int i=1;i<m;i++)
	{
		for(int j=i;j<=n;j++)
			if(r[j].a[i])
			{
				swap(r[i],r[j]);
				break;
			}
		if(r[i].iszero())
		{
			for(int j=i;j<=n;j++)
				if(r[j].a[m]!=0) return NO_SOLUTION;
			return UNLIMITED_SOLUTIONS;
		}
		int fir=0;
		for(int j=1;j<=m;j++)
			if(fabs(r[i].a[j])>1e-6)
			{
				fir=j;
				break;
			}
		for(int j=1;j<=n;j++)
		{
			if(j==i) continue;
			r[j]=r[j]-r[i]*(r[j].a[fir]/r[i].a[fir]);
		}
		r[i]=r[i]/(r[i].a[fir]);
	}
	return SOLVED;
}