Hi,
Am Montag, den 18.12.2006, 22:28 +0100 schrieb [EMAIL PROTECTED]:
> Please run these programm. I dont know where the mistake is.
Run it under DDD and look at the entries of aa_data[] .
> 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);
AA is a view; changes to the elements of AA will change the entries of
the underlying aa_data[]!
So, your aa_data[] changes during your program. The attached version of
your program uses four different arrays to circumvent this (which
obviously is not the way to do it in a real world application).
Thomas
// 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
[email protected] 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
*/
};
_______________________________________________
Help-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/help-gsl
