高斯消元法、LU分解法与克莱姆法则解方程组的C++实现

数值分析老师布置了这样的实习作业：

随机生成一个 $n$ 阶方阵与 $n$ 阶向量构成 $Ax=b$

• 至少生成 $5$ 个方程并取不同的 $n$
• 编程实现克莱姆法则求 $x$
• 编程实现高斯消去法基本方法求 $x$
• 编程实现高斯消去法的 $LU$ 分解求 $x$
• 编程实现列主元高斯消去法求 $x$
• 列表比较并分析以上方法的运算时间与效率

看到这个作业后，我内心其实是拒绝的。要是用 C++ 的话，估计得写好久。

但作业也不能不交是吧，于是就写了下面的代码。

#include<bits/stdc++.h>
using namespace std;
const int N=110;
const double eps=1e-6;
int n,flag,bg,ed;
double d[N],a[N][N],acp[N][N],l[N][N],u[N][N];

void randinit(){ for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) acp[i][j]=(rand()%200)/10.0-10; }
void init(){ for(int i=1;i<=n;i++) for(int j=1;j<=n+1;j++) a[i][j]=acp[i][j]; }
void out(){ printf("x = [ "); for(int i=1;i<=n;i++) printf("%lf ",a[i][n+1]); printf("]\n"); }
void pt(double a[N][N]){ for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) printf("%lf%c",a[i][j]," \n"[j==n]); }

void calc_U(){	//求解上三角方程组
	for(int i=n;i>=1;i--){
		for(int j=n;j>i;j--) a[i][n+1]-=a[i][j]*a[j][n+1];
		a[i][n+1]/=a[i][i];
	}
}

void calc_D(){	//求解下三角方程组
	for(int i=1;i<=n;i++){
		for(int j=1;j<i;j++) a[i][n+1]-=a[i][j]*a[j][n+1];
		a[i][n+1]/=a[i][i];
	}
}

void divide(){	//将A分解为L和U;
	for(int i=1;i<=n;i++){
		if(abs(a[i][i])<eps) return (void)(flag=true);
		for(int k=i+1;k<=n;k++){
			l[k][i]=a[k][i]/a[i][i];
			if(abs(a[k][i])>eps) for(int j=n+1;j>=i;j--)
				a[k][j]-=a[i][j]*(a[k][i]/a[i][i]);
		}
	}
	for(int i=1;i<=n;i++) for(int j=i;j<=n;j++) l[i][j]=i==j;
	for(int i=1;i<=n;i++) for(int j=i;j<=n;j++) u[i][j]=a[i][j];
}

void gauss_base(){	//高斯消去法（基本）
	divide();		//化解为上三角
	calc_U();	//回代
}

void gauss(){	//高斯消去法（列主元、无回代）
	for(int i=1;i<=n;i++){
		int mx=i;
		for(int j=i+1;j<=n;j++) if(a[j][i]>a[mx][i]) mx=j;
		if(abs(a[mx][i])<eps) return (void)(flag=true);
		for(int j=i;j<=n+1;j++) swap(a[i][j],a[mx][j]);
		for(int j=n+1;j>=i;j--) a[i][j]/=a[i][i];
		for(int k=1;k<=n;k++) if(i!=k)
			for(int j=n+1;j>=i;j--) a[k][j]-=a[k][i]*a[i][j];
	}
}

double det(){	//计算det
	double ret=1;
	for(int i=1;i<=n;i++){
		int mx=i;
		for(int j=i+1;j<=n;j++) if(a[j][i]>a[mx][i]) mx=j;
		if(abs(a[mx][i])<eps){ flag=true; return 0; }
		if(mx!=i) ret=-ret;
		ret*=a[mx][i];
		for(int j=i;j<=n+1;j++) swap(a[i][j],a[mx][j]);
		for(int j=n+1;j>=i;j--) a[i][j]/=a[i][i];
		for(int k=1;k<=n;k++) if(i!=k)
			for(int j=n+1;j>=i;j--) a[k][j]-=a[k][i]*a[i][j];
	}
	return ret;
}

void cramer(){	//克莱姆法则求解
	for(int i=0;i<=n;i++){
		init();
		if(i!=0) for(int j=1;j<=n;j++) swap(a[j][i],a[j][n+1]);
		d[i]=det();
	}
	for(int i=1;i<=n;i++) a[i][n+1]=d[i]/d[0];
}


void LU(){	//LU分解法
	divide();
	for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) a[i][j]=l[i][j];
	for(int i=1;i<=n;i++) a[i][n+1]=acp[i][n+1];
	calc_D();
	for(int i=1;i<=n;i++) for(int j=1;j<=n;j++) a[i][j]=u[i][j];
	calc_U();
}


int main(){
	srand(time(0));
	printf("input n (2<=100) :\n");
	scanf("%d",&n);
	do{
		flag=false;
		randinit();
		init();
		divide();
	}while(flag);	//生成有唯一解的方程组
	init();
	bg=clock(); gauss_base();	ed=clock();
	printf("gauss_base ans is : \n"); out();
	printf("gauss_base run time is %d ms\n\n",ed-bg);

	init();
	bg=clock(); gauss();	ed=clock();
	printf("gauss ans is : \n"); out();
	printf("gauss run time is %d ms\n\n",ed-bg);

	init();
	bg=clock(); LU();	ed=clock();
	if(n<=10){
		printf("L is :\n"); pt(l);
		printf("u is :\n"); pt(u);
	}
	printf("LU ans is : \n"); out();
	printf("LU run time is %d ms\n\n",ed-bg);

	bg=clock(); cramer();	ed=clock();
	printf("cramer ans is : \n"); out();
	printf("cramer run time is %d ms\n\n",ed-bg);
}

当 $n$ 为 $100$ 时，四种方法的运行时间如下表：

基本Gauss	列主元Gauss	LU分解法	克莱姆法则
3ms	4ms	3ms	632ms

克莱姆法则的复杂度为 $O(n^4)$，明显慢于其他 $O(n^3)$ 的算法。

刚开始很不理解，$LU$ 分解法的动作明显比高斯消元更冗余。然后搜了搜它们的优缺点，如下：

$LU$ 和 $Gauss$ 都是 $O(n^3)$ 。更精确地讲，乘除法大概都是 $\frac{n^3}{3}$ 次，时间上差别不大，不过：
$LU$ 具有承袭性，这是 $LU$ 的优点。
$LU$ 只适用于解所有顺序主子式都大于 $0$ 的，通用性欠缺，这是 $LU$ 的缺点。
$LU$ 法不保证具有数值稳定性，这是 $LU$ 的缺点。（ $Gauss$ 法可以用选取列主元技巧保证数值稳定性）
集合 $LU$ 与 $Gauss$ 优点，同时规避掉这些缺点的，是 $LUP$ 分解法。

参考链接 LU分解法与Gauss消元法两者的比较。

所谓承袭性，就是说当计算 $Ax=y_i$ ($1\le i\le n$) 时，若使用 $Gauss$ ，那么只能依次调用 $n$ 次，总复杂度为 $O(k\times n^3)$。而使用 $LU$ 分解法时，先用 $O(n^3)$ 的复杂度将原式转化为 $LUx=y_i$ ，之后进行 $k$ 次回代就能求出所有的 $x_i$ ，总复杂度为 $O(n^3+k\times n^2)$。

另外，在调试代码的时候，我最初将数组开的很大，试了试 $200\times 200$ 的矩阵，结果其他方法都能正确求解，唯独克莱姆法则输出全为 $nan$ 。而对于 $100\times 100$ 的矩阵，四种方法都能正确求解。
后来发现，使用克莱姆法则，需要求解行列式的值。而我随机生成的矩阵，这个行列式的值会随着矩阵的规模增大而急剧增大，溢出 $double$。C++党已经输麻了