/*        ATK IV - Numerical Programming (2006)        */
/*    Excercise 8.2 -  Matrix Inverse and Determinant  */

#include<stdio.h>
#include<math.h>
#define NRANSI
#include<nr.h>
#include<nrutil.h>

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];
  }
}