您的位置:首页 > 其它

HDU 3359 提醒下自己!!

2011-08-27 19:13 246 查看
写这道题报告就是提醒下自己,要精心去写代码,去思考,不然检查的时候就要哭了。。

高斯消元

A[][] * x[] = B[];

B是转换后的,x就是一开是的,求出系数阵 然后高斯消元模板就行了

输入是浮点行的,输出要控制,先输入的矩阵的列值,然后是行值和曼哈顿距离

HDU 3359

int n,m,t;
double a[15][15];
double A[105][105],B[105];

//列主元gauss消去求解a[][]x[]=b[]
//返回是否有唯一解,若有解在b[]中
int gauss_cpivot(int n,double a[][MAXN],double b[]){
int i,j,k,row;
double maxp,t;
for (k=0;k<n;k++){
for (maxp=0,i=k;i<n;i++)
if (fabs(a[i][k])>fabs(maxp))
maxp=a[row=i][k];
if (fabs(maxp)<eps)
return 0;
if (row!=k){
for (j=k;j<n;j++)
t=a[k][j],a[k][j]=a[row][j],a[row][j]=t;
t=b[k],b[k]=b[row],b[row]=t;
}
for (j=k+1;j<n;j++){
a[k][j]/=maxp;
for (i=k+1;i<n;i++)
a[i][j]-=a[i][k]*a[k][j];
}
b[k]/=maxp;
for (i=k+1;i<n;i++)
b[i]-=b[k]*a[i][k];
}
for (i=n-1;i>=0;i--)
for (j=i+1;j<n;j++)
b[i]-=a[i][j]*b[j];
return 1;
}

int get_v(int a,int b) {
return a * m + b;    //铭记啊铭记 一开是写的是   a * n + b   ! !!!!!
}

void init() {
int x, y,stmp;

SET(A,0);SET(B,0);
for(int i = 0;i<n;i++) {
for(int j = 0;j<m;j++) {
scanf(" %lf",&a[i][j]);
B[get_v(i,j)] = a[i][j];
}
}

//构造系数阵
for(int i = 0;i<n;i++) {
for(int j = 0;j<m;j++) {
int ith = get_v(i,j);
stmp = 0;

//两种转换方式

/*for(int p = -t ; p <= t;p++) {
for(int q = -t ; q <= t;q++) {
if(abs(p) + abs(q) <= t) {
x = i + p;
y = j + q;
if(x >=0 && x < n && y >=0 && y < m)
A[ith][get_v(x,y)] = 1.0,stmp++;
}

}
}*/
for(int p = 0 ; p <n;p++) {
for(int q = 0 ; q <m;q++) {
if(abs(i - p) + abs(j - q) <= t) {
A[ith][get_v(p,q)] = 1;
stmp ++;
} else {
A[ith][get_v(p,q)] = 0;
}

}
}
//double tmp = 1.0 / (double) stmp;
B[ith] *= stmp;
//for(int k = 0;k<n*m;k++) {
//    if(A[ith][k]) A[ith][k] = tmp;
//}

}
}

gauss_cpivot(n*m,A,B);

for(int i = 0;i < n*m ;i++) {
if(i % m == m - 1) {
printf("%8.2lf\n",B[i]);
} else {
printf("%8.2lf",B[i]);
}
}
}

int main() {
int kk = 1;
while(scanf(" %d%d%d",&m,&n,&t) != EOF) {
if(n == 0 && m == 0 && t == 0) break;
else {

if(kk != 1) puts("");
kk++;
init();
}

}
}


  
内容来自用户分享和网络整理,不保证内容的准确性,如有侵权内容,可联系管理员处理 点击这里给我发消息
标签: