/* ATK IV - Numerical Programming (2006) */ /* Excercise 8.1 - Gaussian Elimination */ #include #include #define NRANSI #include #include void gaussel(float **a, float *b, float *x, int n); void gausspivot(float **a, float *b, float *x, int n); int main(void){ float **a, **c; float *b, *x, *e; int *indx; int n; float d; FILE *file; int i,j; file=fopen("matrix_01.dat","rt"); fscanf(file,"%d\n",&n); a=matrix(1,n,1,n+1); b=vector(1,n); x=vector(1,n); indx = ivector(1,n); c=matrix(1,n,1,n); e=vector(1,n); for(i=1;i<=n;i++){ for(j=1;j<=n;j++){ fscanf(file,"%f",&a[i][j]); } } for(i=1;i<=n;i++){ fscanf(file,"%f",&b[i]); } fclose(file); for(i=1;i<=n;i++){ e[i]=b[i]; for(j=1;j<=n;j++){ c[i][j]=a[i][j]; } } gausspivot(a,b,x,n); ludcmp(c,n,indx,&d); lubksb(c,n,indx,e); for(i=1;i<=n;i++){ printf("%f %f\n",e[i],x[i]); } free_ivector(indx,1,n); free_vector(e,1,n); free_matrix(c,1,n,1,n); free_matrix(a,1,n,1,n); free_vector(b,1,n); free_vector(x,1,n); return 0; } void gaussel(float **a, float *b, float *x, int n){ int i,j,k; for(i=1;i<=n;i++){ a[i][n+1]=b[i]; } for(i=1;i<=n;i++){ for(j=i+1;j<=n+1;j++){ a[i][j]=a[i][j]/a[i][i]; } for(j=i+1;j<=n;j++){ for(k=i+1;k<=n+1;k++){ a[j][k]=a[j][k]-a[j][i]*a[i][k]; } } } x[n]=a[n][n+1]; for(i=n-1;i>=1;i--){ x[i]=a[i][n+1]; for(j=i+1;j<=n;j++){ x[i]=x[i]-a[i][j]*x[j]; } } } void gausspivot(float **a, float *b, float *x, int n){ int i,j,k,ind; float dummy, tiny=1.0e-20; for(i=1;i<=n;i++){ a[i][n+1]=b[i]; } for(i=1;i<=n;i++){ ind=i; for(j=i+1;j<=n;j++){ if(fabs(a[j][i])>fabs(a[ind][i])) ind=j; } if(ind != i){ for(k=i;k<=n+1;k++){ dummy = a[i][k]; a[i][k] = a[ind][k]; a[ind][k] = dummy; } } for(j=i+1;j<=n+1;j++){ a[i][j]=a[i][j]/a[i][i]; } for(j=i+1;j<=n;j++){ for(k=i+1;k<=n+1;k++){ a[j][k]=a[j][k]-a[j][i]*a[i][k]; } } } x[n]=a[n][n+1]; for(i=n-1;i>=1;i--){ x[i]=a[i][n+1]; for(j=i+1;j<=n;j++){ x[i]=x[i]-a[i][j]*x[j]; } } }