#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#include<gsl/gsl_linalg.h>
#include <gsl/gsl_math.h>

int main(int argc, char *argv[]){
  int i, j, ndim, signum;
  float  *A;
  FILE *fp;
  gsl_matrix * B;
  gsl_permutation * p;
  double det; 

  if (argc < 2){
    fprintf(stderr,"./ludecomp <.mat file>\n");
    return(-1);
  }

  fp = fopen(argv[1],"r");
  fscanf(fp,"%d\n",&ndim);
  A = (float *)malloc(ndim*ndim*sizeof(float));
  for(i=0; i < ndim; i++){
    for(j=0; j < ndim; j++)
      fscanf(fp," %f",&A[j+ndim*i]);
    fscanf(fp,"\n");
  }
  fclose(fp);
  
  B =  gsl_matrix_alloc (ndim, ndim);
  p =  gsl_permutation_alloc(ndim);
  
  for(i=0; i < ndim; i++){
    for(j=0; j < ndim; j++){
      printf(" %2.6e ",A[j+ndim*i]);
      gsl_matrix_set(B,i,j,A[j+ndim*i]);
    }
    printf("\n");
  }
  
  gsl_linalg_LU_decomp (B,p,&signum);
 // det = gsl_linalg_LU_det (B,&signum); 
  
//  printf("det = %2.6e\n",det); 
  det = 1.0; 
  for(i=0; i < ndim; i++){
    for(j=0; j < ndim; j++){
  //     printf(" %2.6e ",gsl_matrix_get(B,i,j));
      if (i==j)
       det*=gsl_matrix_get(B,i,j); 

    }
    printf("\n");
  }
  
 if (signum < 0.0)
     det *=-1.0;  

  printf("det = %2.6e\n",det); 
  gsl_permutation_free(p);
  gsl_matrix_free(B);
   
   
  return(0);
}