[Help-gsl] Please run these programm. I dont know where the mistake is.

[buc-family]

Please run these programm. I dont know where the mistake is.

Regards

Andrzej Buczynski

Moosburg, Germany

//   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.  } ;
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 (aa_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 (aa_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 (aa_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

*/
};