#include <stdio.h>
/* No need for these headers.
 */
/* #include <gsl/gsl_math.h>*/ 
/* #include <stdlib.h> */
#include <gsl/gsl_matrix.h>
/* No need to include this header; double is the default type for
 *  matrices (via internal macro hackery).
 */
/*#include <gsl/gsl_matrix_double.h> */

/* Needed for BLAS operations. */
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h> 
#include <iostream>

using namespace std;

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);

    /*save a copy for later. m is about to modify double[] data*/
    gsl_matrix* m2 = gsl_matrix_alloc(n,n);
    gsl_matrix_memcpy(m2,&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 matrix and inverse */
    
    /* Bug lies below. m above has already modified double[] data. m2
     * now doesn't contain [0.1 0.3 ; 0.1 0.25], but rather the LU
     * factorisation [0.1 0.3; 1 -0.05] (using GNU Octave notation for
     * matrices).
     */

    /* gsl_matrix_view m2 = gsl_matrix_view_array(data, n, n); */
    gsl_matrix* tmp = gsl_matrix_alloc(n,n);

    /* Use BLAS function instead. */
    /* gsl_linalg_matmult(m2, inverse, tmp); */
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, m2, inverse, 0, 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);
    gsl_matrix_free(m2);
}