北航 数值分析 矩阵的拟上三角化

#include<stdio.h> #include<math.h> #define N 4 void main() { double sgn(double a);//符号函数 void angle(double A[N+1][N+1]);//拟上三角化 double A[N+1][N+1]={{0,0,0,0,0},{0,2,1,3,4},{0,1,-1,2,1},{0,-1,2,1,2},{0,1,0,-1,3}}; int i,j,k; for(i=1;i<=N;i++) { for(j=1;j<=N;j++) printf("%f",A[i][j]); printf("\n"); } angle(A); /* for(i=1;i<=N;i++) { for(j=1;j<=N;j++) printf("%f",A[i][j]); printf("\n"); }*/ } double sgn(double a)//符号函数 { if(a>0) return 1; else if(a=0) return 0; else return -1; } //拟上三角化 void angle(double A[N+1][N+1]) { double h,t,c,d;double u[N+1],p[N+1],q[N+1],w[N+1];int r,i,j,flag; for(r=1;r<=N-2;r++) { flag=0; //(1) for(i=r+2;i<=N;i++) if(A[i][r]==0) flag=flag+1; if(flag==N-r-1) continue; //(2) d=0; for(i=r+1;i<=N;i++) d=d+A[i][r]*A[i][r]; d=sqrt(d); //printf("%f",d); //printf("\n"); if(A[r+1][r]==0) c=d; else c=(-1)*sgn(A[r+1][r])*d; h=c*c-c*A[r+1][r];//printf("%f",c);printf("\n");// //(3) for(i=1;i<=r;i++) u[i]=0; u[r+1]=A[r+1][r]-c; for(i=r+2;i<=N;i++) u[i]=A[i][r]; for(i=1;i<=N;i++) //{printf("%f",u[i]);printf("\n");} //(4) for(i=1;i<=N;i++) { p[i]=0;q[i]=0; for(j=r+1;j<=N;j++) { {p[i]=p[i]+A[j][i]*u[j];} {q[i]=q[i]+A[i][j]*u[j];} } p[i]=p[i]/h;q[i]=q[i]/h; //printf("%f",p[i]);printf("\n"); } t=0; for(i=r+1;i<=N;i++) t=t+p[i]*u[i]; t=t/h; for(i=1;i<=N;i++) w[i]=q[i]-t*u[i]; for(i=1;i<=N;i++) for(j=1;j<=N;j++) A[i][j]=A[i][j]-w[i]*u[j]-u[i]*p[j]; } for(i=1;i<=N;i++) { for(j=1;j<=N;j++) printf("%f",A[i][j]); printf("\n"); } }