//   check of inversion of real unity matrix UNITY
#include <stdio.h>
#include <math.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_complex.h>
#include <gsl/gsl_complex_math.h>
#include <gsl/gsl_mode.h>
#include <gsl/gsl_permutation.h>
#include<gsl/gsl_permute_complex_double.h>
#include <gsl/gsl_types.h>
#include <gsl/gsl_errno.h>
#include <gsl/gsl_check_range.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_vector_complex.h>
#include <gsl/gsl_block_complex_double.h>
/*
Brian Gough, (GSL Maintainer)  wrote:
You can get help debugging your own program on the 
                address@hidden list
Please use the bug-gsl mailing list only for reporting actual 
bugs in GSL -thanks. 
Note that you do need to call
          gsl_linalg_LU_decomp 
to decompose the matrix before calling any other functions. 
*/
int main (void)
{  
  int N = 4 , i, j, signum ;
  double det, temp ;
  gsl_permutation * p ;
  double aa_data[] = {   -5.,  1., 3.,   2.,
    6.,  0., -2.,  0.,
    -2., 3.,  1.,  -2.,
    4., 5.,   4.,   3.  } ;

  double bb_data[] = {   -5.,  1., 3.,   2.,
    6.,  0., -2.,  0.,
    -2., 3.,  1.,  -2.,
    4., 5.,   4.,   3.  } ;

  double cc_data[] = {   -5.,  1., 3.,   2.,
    6.,  0., -2.,  0.,
    -2., 3.,  1.,  -2.,
    4., 5.,   4.,   3.  } ;

  double dd_data[] = {   -5.,  1., 3.,   2.,
    6.,  0., -2.,  0.,
    -2., 3.,  1.,  -2.,
    4., 5.,   4.,   3.  } ;

  gsl_matrix_view AA = gsl_matrix_view_array (aa_data, N, N);
  /*printf("\n   AA    \n");
    for(i=0; i<N ;i++)  
    {   printf("\n ");
    for(j=0; j<N ;j++)
    {  temp = (double) (gsl_matrix_get( &AA.matrix, i, j )  );
    printf(" %g  " , temp ) ; }; 
    };    printf("\n");    */
  // determinant  of AA
  p = gsl_permutation_alloc (N); 
  gsl_linalg_LU_decomp (&AA.matrix, p, &signum );
  det = gsl_linalg_LU_det ( &AA.matrix, signum )  ;
  printf("\n     signum = %d \n", signum ) ; 
  printf("\n    det ( A )= %g  \n", det );   // det aa = -90  it is OK

  // inversion of matrix AA = inverseAA   
  gsl_matrix  *inverseAA =  gsl_matrix_calloc( N, N );
  AA = gsl_matrix_view_array (bb_data, N, N) ;
  p = gsl_permutation_alloc (N);
  gsl_linalg_LU_decomp (&AA.matrix, p, &signum );
  gsl_linalg_LU_invert (&AA.matrix ,  p,  inverseAA ) ;
  //free memory
  gsl_permutation_free (p);

  printf("\n    inverseAA    \n");
    for(i=0; i<N ;i++)  
    {   printf("\n ");
    for(j=0; j<N ;j++)
    {  temp = (double) (gsl_matrix_get( inverseAA,i,j)  );
    printf(" %g  " , temp ) ; }; 
    };     printf("\n");  
  
    /*   Compute C = AA * inverseAA   */

  AA = gsl_matrix_view_array (cc_data, N, N) ;
  gsl_matrix  *C =  gsl_matrix_calloc( N, N );
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
      1.0,  &AA.matrix, inverseAA, 
      0.0, C );
  printf("\n   C = AA * inverseAA   should be UNIT matrix \n");
  for(i=0; i<N ;i++)  
  {   printf("\n ");
    for(j=0; j<N ;j++)
    {  temp = (double) (gsl_matrix_get( C,i,j)  );
      printf(" %g " , temp ) ; }; 
  };  printf("\n");


  AA = gsl_matrix_view_array (dd_data, N, N) ;
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
      1.0,  inverseAA, &AA.matrix, 
      0.0, C );

  printf("\n   C =  inverseAA * AA  should be UNIT matrix  \n");
  for(i=0; i<N ;i++)  
  {   printf("\n ");
    for(j=0; j<N ;j++)
    {  temp = (double) (gsl_matrix_get( C,i,j)  );
      printf(" %g " , temp ) ; }; 
  };  printf("\n");


  // free memory
  gsl_matrix_free (inverseAA) ; 
  gsl_matrix_free ( C ) ; 
  return 0;
  /*

     Matrix A 
     1 0 0 0 1 
     0 1 0 0 2 
     0 0 1 0 3 
     0 0 0 1 4 
     1 2 3 4 5 

     B = inverse( A )

     Matrix C = A * B:
     1 0 0 0 0 
     0 1 0 0 0 
     0 0 1 0 0 
     0 0 0 1 0 
     0 0 0 0 1 
     If the matrix C = A * B is unity then check is OK

*/
};