SIRT.rar_M?n_NSIRT119_sirt_sirt迭代

#include<stdio.h> #include<math.h> #define L 12 #define W 9 #define M L*W //108 #define N 144 //方程个数，可变657 #define e 0.000001 //精度 void SIRT(double A[N][M], double b[N], double x[M]); //SIRT算法函数 void main() { FILE *fp1,*fp2; int i,j; double x[M],b[N]; double A[N][M],V[L][W]; /* 读入系数矩阵 */ fp1=fopen("S_T.txt","r"); for(i=0;i<N;i++) { for(j=0;j<=M;j++) { if(j!=M) { fscanf(fp1,"%lf",&A[i][j]); //printf("%f",A[i][j]); //if(fmod(double(j+1),double(W))==0.0) // printf("\n"); } else { fscanf(fp1,"%lf",&b[i]); //printf("b[%d]=%f\n",i,b[i]); } } } fclose(fp1); /* 调用SIRT算法函数 */ printf("程序运行中......\n"); SIRT(A,b,x); /* 打印反演结果并存盘 */ fp2=fopen("R_SIRT.txt","w"); for(i=0;i<L;i++) { for(j=0;j<W;j++) { V[i][j]=1.0/x[i*W+j]; fprintf(fp2,"%-5.2f ",V[i][j]); printf("%5.2f ",V[i][j]); } fprintf(fp2,"\n"); printf("\n"); } /* for(i=0;i<100;i++) { fprintf(fp2,"%f\n",x[i]); printf("%f\n",x[i]); }*/ fclose(fp2); } ////////////////////////////////////////////////////// void SIRT(double A[N][M], double b[N], double x[M]) { int i,j,k,l,m,flag; double r[N]; double ex[M]; double al,bt,sum; for(i=0;i<M;i++) x[i]=0.0; for(i=0;i<N;i++) { sum=0.0; for(j=0;j<M;j++) sum+=A[i][j]*x[j]; r[i]=b[i]-sum; } l=0; do { flag=0; for(j=0;j<M;j++) { //第j列非零元素的个数m m=0; for(k=0;k<N;k++) { if(A[k][j]!=0.0) m++; } sum=0.0; for(i=0;i<N;i++) { al=A[i][j]*r[i]; bt=0.0; for(k=0;k<M;k++) bt+=A[i][k]*A[i][k]; sum+=al/bt; } ex[j]=sum/double(m); /////////////// if(ex[j]>e) flag=1; x[j]=x[j]+ex[j]; } for(i=0;i<N;i++) { sum=0.0; for(j=0;j<M;j++) sum+=A[i][j]*x[j]; r[i]=b[i]-sum; } l++; }while(flag==1); printf("迭代次数:l=%d\n",l); }