#include <stdio.h>
#include "/usr/include/gsl/gsl_math.h"
#include "/usr/include/gsl/gsl_linalg.h"

int main() {
    int n, i, j, k, s, pi;
    double val, sum, Lik, Ukj;
    FILE *fp;

    n = 4;
    gsl_matrix *A = gsl_matrix_calloc(n,n);
    gsl_matrix *Ainv = gsl_matrix_calloc(n,n);
    gsl_vector *b = gsl_vector_calloc(n);
    gsl_vector *x = gsl_vector_calloc(n);
    gsl_permutation *p = gsl_permutation_calloc(n);

    if ((fp = fopen("A.txt", "rt")) == NULL) {
        printf ("\n*** Could not open file A.txt ***\n");
        return 1;
    }

    printf("\n*** Input: A ***\n");
    for (i=0; i < n; i++) {
        for(j=0; j < n; j++) {
            fscanf(fp, "%lf", &val);
            gsl_matrix_set(A,i,j,val);
            printf("%lf ", gsl_matrix_get(A,i,j));
        }
	printf("\n");
    }

    fclose(fp);

    gsl_linalg_LU_decomp(A, p, &s);

    for (i=0; i < n; i++) gsl_vector_set(b,i,0);
    for (i=0; i < n; i++) {
	gsl_vector_set(b,i,1);
	if (i > 0) gsl_vector_set(b,i-1,0);
	gsl_linalg_LU_solve(A, p, b, x);
	for (j=0; j < n; j++) gsl_matrix_set(Ainv,j,i,gsl_vector_get(x,j));
    }

    printf("\n*** Output: Ainv ***\n");
    for (i=0; i < n; i++) {
	for (j=0; j < n; j++)
	    printf("%lf ", gsl_matrix_get(Ainv,i,j));
	printf("\n");
    }
    printf("\n");

    return 0;
}