/* ATK IV - Numerical Programming (2006) */ /* Excercise 8.2 - Matrix Inverse and Determinant */ #include #include #define NRANSI #include #include void luinv(float **a, float **ai, int *indx, int n); void ludet(float **a, int n, float *d); int main(void){ float **a,**ai; int *indx; static int n; float d; FILE *file; int i,j; file=fopen("matrix_01.dat","rt"); fscanf(file,"%d\n",&n); indx = ivector(1,n); a=matrix(1,n,1,n); ai=matrix(1,n,1,n); for(i=1;i<=n;i++){ for(j=1;j<=n;j++){ fscanf(file,"%f",&a[i][j]); } } fclose(file); ludcmp(a,n,indx,&d); luinv(a,ai,indx,n); ludet(a,n,&d); for(i=1;i<=n;i++){ for(j=1;j<=n;j++){ printf("%f ",ai[i][j]); } printf("\n"); } printf("%f\n",d); free_ivector(indx,1,n); free_matrix(a,1,n,1,n); return 0; } void luinv(float **a, float **ai, int *indx, int n){ int i,j; float *b; b=vector(1,n); for(i=1;i<=n;i++){ for(j=1;j<=n;j++){ ai[i][j]=0.0; } ai[i][i]=1.0; } for(j=1;j<=n;j++){ for(i=1;i<=n;i++){ b[i]=ai[i][j]; } lubksb(a,n,indx,b); for(i=1;i<=n;i++){ ai[i][j]=b[i]; } } free_vector(b,1,n); } void ludet(float **a, int n, float *d){ int i; for(i=1;i<=n;i++){ *d=*d*a[i][i]; } }