高斯赛德尔迭代压缩稀疏矩阵存储后的求解算法

#include "math.h" #include "stdio.h" #include "stdlib.h" int agsdl(a,b,n,x,eps) int n; double a[],b[],x[],eps; { int i,j,u,v,k; double p,t,s,q, aii; int d=1;//记录非零元素的个数 int l; int real_j;//记录真实的列号 int dd=1; int *kd, *d_row, *d_colum; FILE *fd_a, *fd_row , *fd_colum, *fkd; double *d_a; kd = (int*)malloc(n*sizeof(int)); kd[0] = 0; for(i=1;i<=n-1;i++) { int temp = 0; for(j=0;j<=i;j++) { if(a[i*n+j]==0) continue; else//(a[i*n+j]!=0) { //temp = temp+1; //kd[i] = kd[i-1] +temp; kd[i] = kd[i-1]+i-j+1;//建立kd矩阵 break; } } } fkd = fopen("kd.dat", "w+"); for(k=0; k<n; k++) { fprintf(fkd, "%d \n", kd[k]); } // for(i = 1; i<n; i++)//得到下三角矩阵中a非零元素的个数 // { // l = kd[i]-kd[i-1]; // for( j=0; j<l; j++) // { // if (a[kd[i-1]+j+1]!=0.0) // { // d = d+1; // } // } // } for(i = 1; i<n; i++)//得到下三角矩阵中a非零元素的个数 { for( j=0; j<=i; j++) { if (a[i*n+j]!=0.0) { d = d+1; } } } d_a = (double*)malloc(d*sizeof(double)); d_row = (int*)malloc(d*sizeof(int)); d_colum = (int*)malloc(d*sizeof(int)); d_a[0] = a[0]; d_row[0] = 0; d_colum[0] =0; //for( i = 1; i<n; i++) // { // l = kd[i]-kd[i-1]; // real_j = 0; // for( j=0; j<l; j++) // { // int temp = kd[i-1]+j+1; // real_j = real_j + 1; // // if (a[kd[i-1]+j+1]!=0.0) // { // dd = dd+1; // //d_a[dd-1] = a[kd[i-1]+j+1+i*n]; // d_a[dd-1] = a[j+i*n]; // d_row[dd-1] = i; // d_colum[dd-1] = i-l+j+1; // } // } // } for( i = 1; i<n; i++) { for( j=0; j<=i; j++) { if (a[i*n+j]!=0.0) { dd = dd+1; //d_a[dd-1] = a[kd[i-1]+j+1+i*n]; d_a[dd-1] = a[j+i*n]; d_row[dd-1] = i; d_colum[dd-1] = j; } } } fd_a = fopen("fd_a.txt", "w+"); for (i=0; i<d; i++) { fprintf(fd_a, "%f \n", d_a[i]); } fd_row = fopen("fd_row.dat", "w+"); for (i=0; i<d; i++) { fprintf(fd_row, "%d \n", d_row[i]); } fd_colum = fopen("fd_colum.dat", "w+"); for (i=0; i<d; i++) { fprintf(fd_colum, "%d \n", d_colum[i]); } for (i=0; i<=n-1; i++) { u=i*n+i; p=0.0; x[i]=0.0; for (j=0; j<=n-1; j++) if (i!=j) { v=i*n+j; p=p+fabs(a[v]); } if (p>=fabs(a[u])) { printf("fail\n"); return(-1);//不是严格对角占优矩阵，计算失败 } } while (p>=eps) { p=0.0; for (i=0; i<=n-1; i++) { t=x[i]; s = 0.0; for(j=0; j<d; j++) { if(d_row[j]==i) { if(d_colum[j]!=i) s = s + d_a[j]* x[d_colum[j]]; } if(d_colum[j]==i) { if(d_row[j]!=i) s = s +d_a[j]* x[d_row[j]]; if(d_row[j]==i) aii =d_a[j]; } } x[i] = (b[i]-s)/aii; q=fabs(x[i]-t)/(1.0+fabs(x[i])); if (q>p) p=q; } } // while (p>=eps) // { // p=0.0; // for (i=0; i<=n-1; i++) // { // t=x[i]; s=0.0; // for (j=0; j<=n-1; j++) // { // if (j!=i) // s=s+a[i*n+j]*x[j]; // } // x[i]=(b[i]-s)/a[i*n+i]; // q=fabs(x[i]-t)/(1.0+fabs(x[i])); // if (q>p) // p=q; // } // } // p=eps+1.0; return(1); }