gsl matrix invert
gsl_matrix *
invert_a_matrix(gsl_matrix *matrix)
{
gsl_permutation *p = gsl_permutation_alloc(size);
int s;
// Compute the LU decomposition of this matrix
gsl_linalg_LU_decomp(matrix, p, &s);
// Compute the inverse of the LU decomposition
gsl_matrix *inv = gsl_matrix_alloc(size, size);
gsl_linalg_LU_invert(matrix, p, inv);
gsl_permutation_free(p);
return inv;
}