1 #include <iostream> 2 #include <cstdio> 3 #include <vector> 4 #include <cmath> 5 using namespace std; 6 const double EPS=1e-8; 7 8 typedef vector<double> vec; 9 typedef vector<vec> mat; 10 11 vec Gauss(const mat& A,const vec& b){ 12 int n=A.size(); 13 mat B(n,vec(n+1)); 14 for(int i=0;i<n;++i) 15 for(int j=0;j<n;++j) B[i][j]=A[i][j]; 16 for(int i=0;i<n;++i) B[i][n]=b[i]; 17 for(int i=0;i<n;++i){ 18 int pivot=i; 19 for(int j=i+1;j<n;++j)if(abs(B[j][i]>abs(B[pivot][i]))) pivot=j; 20 swap(B[i],B[pivot]); 21 if(abs(B[i][i])<EPS) return vec(); 22 for(int j=i+1;j<=n;++j) B[i][j]/=B[i][i]; 23 for(int j=i+1;j<n;++j) 24 if(j!=i) for(int k=i+1;k<=n;++k) B[j][k]-=B[j][i]*B[i][k]; 25 } 26 vec x(n); 27 for(int i=0;i<n;++i) x[i]=B[i][n]; double ans; 28 for(int i=n-1;i>=0;--i){ 29 ans=x[i]; 30 for(int j=i+1;j<n;++j) ans-=B[i][j]*x[j];x[i]=ans; 31 } 32 return x; 33 } 34 double ai[3][3]={{1,-2,3},{4,-5,6},{7,-8,10}}; 35 double bi[]={6,12,21}; 36 int main(){ 37 mat A(3,vec(3)); 38 vec b(3); 39 for(int i=0;i<3;++i){ 40 for(int j=0;j<3;++j){ 41 A[i][j]=ai[i][j]; 42 } 43 } 44 for(int i=0;i<3;++i) b[i]=bi[i]; 45 vec x(3); 46 x=Gauss(A,b); 47 for(int i=0;i<x.size();++i){ 48 printf("%f,",x[i]); 49 } 50 printf("\n"); 51 return 0; 52 }
可能有一个令你疑惑的地方是为什么22行，j从i+1开始除，30行，j从i+1开始除
因为这些地方其实是被搞成系数为1的，我们在计算的时候不会用到这些值，所以我们不用去管当前消去的这一变量
最后一个步骤容易忽略，那就是回带的步骤，因为这个方法算完得到的B矩阵(不含增广，是一个上三角行列式，所以我们不断迭代就好了
这个方法的复杂度是n^3的
对于1000规模的，我们就需要剪枝了....
下面附上剪枝后的模板...sdut2878这个题需要剪枝
1 #include <iostream> 2 #include <cstdio> 3 #include <vector> 4 #include <cmath> 5 using namespace std; 6 const double EPS=1e-8; 7 8 typedef vector<double> vec; 9 typedef vector<vec> mat; 10 11 vec Gauss(const mat& A,const vec& b){ 12 int n=A.size(); 13 mat B(n,vec(n+1)); 14 for(int i=0;i<n;++i) 15 for(int j=0;j<n;++j) B[i][j]=A[i][j]; 16 for(int i=0;i<n;++i) B[i][n]=b[i]; 17 for(int i=0;i<n;++i){ 18 int pivot=i; 19 for(int j=i+1;j<n;++j)if(abs(B[j][i]>abs(B[pivot][i]))) pivot=j; 20 if(pivot!=i) swap(B[i],B[pivot]); 21 if(abs(B[i][i])<EPS) return vec(); 22 for(int j=i+1;j<=n;++j) B[i][j]/=B[i][i]; 23 for(int j=i+1;j<n;++j) 24 if(j!=i) { 25 if(abs(B[j][i])<EPS) continue; 26 for(int k=i+1;k<=n;++k) B[j][k]-=B[j][i]*B[i][k]; 27 } 28 } 29 vec x(n); 30 for(int i=0;i<n;++i) x[i]=B[i][n]; double ans; 31 for(int i=n-1;i>=0;--i){ 32 ans=x[i]; 33 for(int j=i+1;j<n;++j) ans-=B[i][j]*x[j];x[i]=ans; 34 } 35 return x; 36 } 37 double ai[3][3]={{1,-2,3},{4,-5,6},{7,-8,10}}; 38 double bi[]={6,12,21}; 39 int main(){ 40 mat A(3,vec(3)); 41 vec b(3); 42 for(int i=0;i<3;++i){ 43 for(int j=0;j<3;++j){ 44 A[i][j]=ai[i][j]; 45 } 46 } 47 for(int i=0;i<3;++i) b[i]=bi[i]; 48 vec x(3); 49 x=Gauss(A,b); 50 for(int i=0;i<x.size();++i){ 51 printf("%f,",x[i]); 52 } 53 printf("\n"); 54 return 0; 55 }
20,25行是剪枝的地方