[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] Please run these programms.
From: |
buc-family |
Subject: |
[Bug-gsl] Please run these programms. |
Date: |
Sun, 17 Dec 2006 22:40:46 +0100 |
User-agent: |
Thunderbird 1.5.0.8 (X11/20060911) |
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;
};
- [Bug-gsl] Please run these programms.,
buc-family <=