题解 P3389 【【模板】高斯消元法】
YellowBean_Elsa · · 题解
截至2020年1月26日13：00，此题所有题解都多少用到了高斯消元
那我就来一发非高斯消元的解法吧！
LU分解！！！
1. 定义及求法
即把系数矩阵A分解成
一个对角线上元素都是1的下三角矩阵L和
一个上三角矩阵U的乘积
利用矩阵乘法法则可推出求LU分解的递推公式（如上图）
要注意递推顺序：
-
L的对角线元素全为1
-
用矩阵乘法法则直接求出U第一行，再是L第一列
-
利用递推公式，一行一列交错，U和L交错地求
2. LU分解的作用
-
求行列式：
det(A)=det(L*U)=det(L)*det(U) 由上下三角矩阵的行列式公式（对角元素之积）可快速求得A的行列式
-
解线性方程组
AX=B (X为存方程解的矩阵)
建一个临时矩阵
我们根据
因为L（U）为下（上）三角矩阵，关于它们的方程很好解
3. 注意事项
看图中的递推式，有个
若
直接 No Solution !
因为这意味着U行列式值为0，A也是
而行列式值为0，就是 No Solution.
代码
//coder: Feliks Hacker of IOI == YB an AKer of IMO
//LU分解
//https://blog.csdn.net/yang10560/article/details/77948021
#include<bits/stdc++.h>
#define db double
#define ll long long
using namespace std;
const int N=216;
int n;
int a[N][N],b[N];
db l[N][N],u[N][N];
//递推前两步
inline void init(){
for(int i=1;i<=n;i++){
l[i][i]=1;
u[1][i]=(db)a[1][i];
//特判
if(u[1][1]==0){
puts("No Solution");
exit(0);
}
l[i][1]=(db)a[i][1]/u[1][1];
}
}//行列交错递推
inline void lu_div(){
init();
//递推式
//U[i][j] = A[i][j] - (for (k = 1 to i-1) L[i][k] * U[k][j])
//L[j][i] = (A[j][i] - (for (k = 1 to i-1) L[j][k] * U[k][i])) / U[i][i]
for(int i=2;i<=n;i++){
for(int j=i;j<=n;j++){
//递推U第i行
u[i][j]=a[i][j];
for(int k=1;k<=i-1;k++)
u[i][j]-=l[i][k]*u[k][j];
//递推L第i列
l[j][i]=a[j][i];
for(int k=1;k<=i-1;k++)
l[j][i]-=l[j][k]*u[k][i];
//特判
if(u[i][i]==(db)0){
puts("No Solution");
exit(0);
}
l[j][i]/=u[i][i];
}
}
}//求行列式的值
inline db det(){
//det(A) = det(L * U) = det(L) * det(U) = det(U)
db res=1;
for(int i=1;i<=n;i++)res*=u[i][i];
return res;
}//解 Z相关的那两个方程
//直接用代入消元一个个求即可
db z[N];
//L * Z = B
inline void solve_z(){
for(int i=1;i<=n;i++){
z[i]=b[i];
for(int j=1;j<i;j++)
z[i]-=l[i][j]*z[j];
}
}
//U * X = Z
db x[N];
inline void solve_x(){
for(int i=n;i>=1;i--){
x[i]=z[i];
for(int j=n;j>i;j--)
x[i]-=u[i][j]*x[j];
x[i]/=u[i][i];
}
}
int main(){
scanf("%d",&n);
for(int i=1;i<=n;i++){
for(int j=1;j<=n;j++){
scanf("%d",&a[i][j]);
}scanf("%d",&b[i]);
}lu_div();
//行列式 = 0，方程无唯一解
if(det()==0){puts("No Solution");return 0;}
solve_z();solve_x();
for(int i=1;i<=n;i++)
printf("%.2lf\n",x[i]);
return 0;
}