[Bug-gsl] Please run these programms.

[buc-family]

Hi,
I don't know, it's bug or it's my mistake.
I'm sending now three programms.  The output is buggy.

Thanks for taking the time.
 
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>
int main (void)
{  
   int N = 4 ;
   int i, j ; // current variable
   double det ;
   double temp ;
    int signum ;
   gsl_permutation * p ;

double aa_data[] = {   -5.,  1., 3.,   2.,
                                   6.,  0., -2.,  0.,
                                  -2., 3.,  1.,  -2.,
                                   4., 5.,   4.,   3.  } ;
// det aa = -90 

gsl_matrix_view AA = gsl_matrix_view_array (aa_data, 4, 4);
   det = gsl_linalg_LU_det ( &AA.matrix, signum )  ;
   printf("\n    det ( A )= %g  \n", det );
   printf("\n     signum = %d \n", signum ) ;

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");  printf("\n    ");

gsl_matrix  *inverseAA =  gsl_matrix_calloc( N, N );

// inversion of matrix AA = inverseAA
    p = gsl_permutation_alloc (N); 
   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");  printf("\n    ");

 gsl_matrix  *C =  gsl_matrix_calloc( N, N );
/* Compute C = A *inverseA */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                1.0,  &AA.matrix, inverseAA, 
                0.0, C );

printf("\n   C = AA * inverseAA    \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");  printf("\n    ");
// free memory

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

Matrix B = inv(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 
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

*/
};





//   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>

/*   STEP 1
 matrix B :
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 

STEP 2
matrix  inverseB  from routine gsl_linalg_LU_invert ( ... )

STEP 3
Matrix C = inverseB * 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 = inverseA * B is unity then check is OK
*/

int main (void)
{  
   int N = 5 ;
   int i, j ; // current variable
   double det ;
   double temp ;
    int signum ;
   gsl_permutation * p ;

//STEP 1
double B_data[] = { 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 } ; 
gsl_matrix_view B = gsl_matrix_view_array (B_data, N, N);
printf("\n  control print of matrix  B    \n");
    for(i=0; i<N ;i++)  
          {   printf("\n ");
              for(j=0; j<N ;j++)
                 {  temp = (double) (gsl_matrix_get( &B.matrix, i, j )  );
                   printf(" %g  " , temp ) ; }; 
                 };
printf("\n");  printf("\n    ");
// determinant of B
   det = gsl_linalg_LU_det ( &B.matrix, signum )  ;
   printf("\n    det ( B )= %g  \n", det );
   printf("\n     signum = %d \n", signum ) ;

//  STEP 2

gsl_matrix  *inverseB =  gsl_matrix_calloc( N, N );

// inversion of matrix B = inverseB
    p = gsl_permutation_alloc (N); 
   gsl_linalg_LU_invert (&B.matrix ,  p,  inverseB ) ;
   //free memory
   gsl_permutation_free (p);

printf("\n   control print of  inverseB    \n");
    for(i=0; i<N ;i++)  
          {   printf("\n ");
              for(j=0; j<N ;j++)
                 {  temp = (double) (gsl_matrix_get( inverseB,i,j)  );
                   printf(" %g  " , temp ) ; }; 
                 };

// determinant of inverseB
   det = gsl_linalg_LU_det ( inverseB, signum )  ;
   printf("\n    det ( inverseB )= %g  \n", det );
   printf("\n     signum = %d \n", signum ) ;

printf("\n");  printf("\n    ");

// STEP 3

 gsl_matrix  *C =  gsl_matrix_calloc( N, N );
/* Compute C = B *inverseB */
  gsl_blas_dgemm (CblasNoTrans, CblasNoTrans,
                1.0,  &B.matrix, inverseB, 
                0.0, C );

printf("\n control print of C=AA*inverseAA,     C  should be unity   \n");
    
    for(i=0; i<N ;i++)  
          {   printf("\n ");
              for(j=0; j<N ;j++)
                 {  temp = (double) (gsl_matrix_get( C,i,j)  );
                   printf(" %f  " , temp ) ; }; 
                 };
printf("\n");  printf("\n    ");
// free memory

gsl_matrix_free (inverseB) ; 
gsl_matrix_free ( C ) ; 
return 0;
};





#include <stdio.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>
#include <gsl/gsl_rng.h>
/*     produce uniform random numbers 
        in the range [0.0, 1.0)   */
int main (void)
{
   int N =3  ; // matrix N x N
   int i, j ; // current variable
   gsl_complex temp , temp2 ; // temporary variable
   gsl_matrix_complex  *A =  gsl_matrix_complex_alloc(N, N );
   gsl_matrix_complex  *inverseA =  gsl_matrix_complex_alloc(N, N );
   gsl_matrix_complex  *inverseinverseA =  gsl_matrix_complex_alloc(N, N );

// setting of  random matrix A

   const gsl_rng_type * T;
   gsl_rng * r;
   gsl_rng_env_setup();
   T = gsl_rng_default;
   r = gsl_rng_alloc (T);
   double u1, u2 ;  // temporary variable
   for(i = 0; i <N ; i++) 
     for(j = 0; j <N ; j++) 
           { u1= gsl_rng_uniform (r) ; 
         u2= gsl_rng_uniform (r) ; 
         gsl_matrix_complex_set(A,i,j,gsl_complex_rect(u1,u2));
           }; 
   
// determinant 
    
gsl_complex det ;
int signum;
gsl_linalg_complex_LU_det ( A,  signum ) ;

printf("\n  det A = %g +i * %g \n", GSL_REAL(det) , GSL_IMAG(det) ) ;


// inversion of matrix A = inverseA
  
   {
   gsl_permutation * p = gsl_permutation_alloc (N); 
   gsl_linalg_complex_LU_invert ( A ,  p,  inverseA ) ;
   //free memory
   gsl_permutation_free (p);
   }
  
// inversion of inversion of A = inverseinverseA
 
    {
     gsl_permutation * p = gsl_permutation_alloc (N); 
     gsl_linalg_complex_LU_invert ( inverseA ,  p,  inverseinverseA ) ;
    //free memory
     gsl_permutation_free (p);
    } 
  
printf("\n\n   A=A-(inverse_inverseA)  should be  ZERO   \n");
    for(i=0; i<N ;i++)  
          { printf("\n");
        for(j=0; j<N ;j++)
                  { temp  =  gsl_matrix_complex_get ( inverseinverseA, i, j);
            temp2 = gsl_matrix_complex_get ( A, i, j);
printf("A(%d,%d)=%f+i*%f \n", i, j, GSL_REAL(gsl_complex_sub(temp2, temp)),
                                     GSL_IMAG(gsl_complex_sub(temp2, temp) ) ); 
                  }; 
          };
printf("\n");
// free memory

gsl_matrix_complex_free (A) ;
gsl_matrix_complex_free (inverseA) ;
gsl_matrix_complex_free (inverseinverseA) ;
gsl_rng_free (r);
return 0;
};