Posted on 2007-05-28 19:26
oyjpart 阅读(2359)
评论(0) 编辑 收藏 引用 所属分类:
ACM/ICPC或其他比赛
一直都没写过高斯消元 今天数值分析老师说要写一个交上去 我就写了一个 呵呵
//by OyjpArt
//使用方法:
//ok = 0;
//solve(0);
const int N = 1010;
double mat[N][N]; //增广矩阵
double ans[N];
int n; //未知数个数
const double EPS = 1e-7;
bool ok;
int dblcmp(double a) { if(fabs(a) < EPS) return 0; if(a < 0) return -1; return 1; }
void solve(int x) {
if(x == n-1) {
if(dblcmp(mat[x][x]) == 0) ok = 0;
else ans[x] = mat[x][x+1] / mat[x][x];
return;
}
int i, j;
for(i = x; i < n && dblcmp(mat[i][x]) == 0; i++);
if(i == n) { ok = 0; return; }
if(i != x) {
double tmp[N];
memcpy(tmp, mat[x], (n+1) * sizeof(double));
memcpy(mat[x], mat[i], (n+1) * sizeof(double));
memcpy(mat[i], tmp, (n+1) * sizeof(double));
}
for(i = x+1; i < n; i++) {
if(dblcmp(mat[i][x]) == 0) continue;
double m = mat[x][x] / mat[i][x];
for(j = x; j < n + 1; j++)
mat[i][j] = mat[i][j] * m - mat[x][j];
}
solve(x+1);
double sum = mat[x][n];
for(i = x+1; i < n; i++) sum -= mat[x][i] * ans[i];
ans[x] = sum / mat[x][x];
}