1. 程式人生 > 其它 >[模板]高斯消元

[模板]高斯消元

#include<cstdio>
#include<cstring>
#include<string>
#include<cmath>
#define WR WinterRain
using namespace std;
const int WR=1010;
const double eps=1e-6;//焯!!!!!!!!精度不是int!!!!!!! 
int equ,var;//有equ個方程,var個變元
double a[WR][WR];//增廣矩陣 
double x[WR];//存解集 
bool fre[WR];//標記是否為自由元,即可隨意取值 
bool
f=true; int read(){ int s=0,w=1; char ch=getchar(); while(ch>'9'||ch<'0'){ if(ch=='-') w=-1; ch=getchar(); } while(ch>='0'&&ch<='9'){ s=(s<<3)+(s<<1)+ch-48; ch=getchar(); } return s*w; } void gauss(){ /*轉化成階梯矩陣
*/ int row=1; for(int col=1;col<=var;col++){ //row:當前處理的行 //col:當前處理的列 int mxrow=row;//尋找當前列係數絕對值最大的行 for(int i=row+1;i<=equ;i++){//沒有必要去上面找,上面的都已經處理好了 if(fabs(a[i][col])>fabs(a[mxrow][col])){ mxrow=i; } }
if(fabs(a[mxrow][col])<eps){//col列從第row行以下全都為0,處理當前行的下一列 // printf("111"); continue; } //交換 for(int i=row;i<=var+1;i++){//手動模擬高斯消元可以發現 //前面的幾項(從1到row-1)都為0了,不用交換 swap(a[row][i],a[mxrow][i]); } for(int i=row+1;i<=equ;i++){//手模一遍相當於消去該行的第col列的元 if(fabs(a[i][col])>eps){ double tmp=a[i][col]/a[row][col];//乘數 for(int j=col+1;j<=var+1;j++){ a[i][j]-=a[row][j]*tmp; } //a[i][col]=0.0; } } row++; } //printf("%d\n",row); //for(int i=1;i<=equ;i++){ // for(int j=1;j<=equ+1;j++){ // printf("%.2lf ",a[i][j]); // } // printf("\n"); //} if(row<=equ){ while(row<=equ){ if(a[row++][var+1]!=0){ printf("-1"); f=false; return; } } printf("0"); f=false; return; } for(int i=equ;i>=1;i--){ double tmp=a[i][var+1]; for(int j=i+1;j<=var;j++){ tmp-=a[i][j]*x[j]; } x[i]=tmp/a[i][i]; //printf("%d %lf\n",i,x[i]); } } int main(){ equ=var=read();//這裡讓var=equ(看題面) for(int i=1;i<=equ;i++){ for(int j=1;j<=equ+1;j++){ int t=read(); a[i][j]=(double)t; } } gauss(); if(!f) return 0; for(int i=1;i<=var;i++){ if(fabs(x[i])<eps) printf("x%d=0\n",i); else printf("x%d=%.2lf\n",i,x[i]); } return 0; }