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

int main(int argn, char*argv[])
{
    /* construct a simple matrix */
    int n=2;
    double data[4] = {0.1,0.3,0.1,0.25}; 
    gsl_matrix_view m = gsl_matrix_view_array(data, n, n);

    /* Make a copy of the original matrix. */
    gsl_matrix* morig = gsl_matrix_alloc(n,n);
    gsl_matrix_memcpy(morig,&m.matrix);

    /* calculate the inverse */
    gsl_matrix*inverse = gsl_matrix_alloc(n, n);
    gsl_permutation*perm = gsl_permutation_alloc(n);
    int s=0;
    gsl_linalg_LU_decomp(&m.matrix, perm, &s);
    gsl_linalg_LU_invert(&m.matrix, perm, inverse);

    /* multiply original matrix and inverse */
    gsl_matrix* tmp = gsl_matrix_alloc(n,n);
    gsl_linalg_matmult(morig, inverse, tmp);

    /* the following should now output the product of matrix
       and inverted matrix- a unit matrix, if the inverse is
       correct */
    int x,y;
    for(y=0;y<n;y++) {
	for(x=0;x<n;x++) {
	    printf("%5.3f\t", gsl_matrix_get(tmp, y, x));
	}
	printf("\n");
    }

    /* free up memory */
    gsl_permutation_free(perm);
    gsl_matrix_free(inverse);
    gsl_matrix_free(tmp);
}