Re: [Help-gsl] gsl_linalg_LU_invert

[Philip Ogunbona]

Warren,
I found my mistake, I had already transposed the matrix in-place using
gsl_matrix_transpose() at a point in the source code before using the
gsl_linalg_LU_invert.
I'm sorry about being a bother and thanks for the help.

Regards,
Philip

Warren Weckesser wrote:
> Philip,
> 
> How do you set up the matrix to be inverted? Are you simply reversing
> the interpretation of row and column? Could you show your code?
> 
> This code shows an example of computing an inverse with the LU
> library function:
> --------------------
> #include <stdio.h>
> #include <gsl/gsl_linalg.h>
>     
> int
> main (void)
>     {
>     double a_data[] = { 1.0, 0.6, 0.0,
>                         0.0, 1.5, 1.0,
>                         0.0, 1.0, 1.0 };
>     /*
>      * Inverse is
>      *    1  -1.2   1.2
>      *    0   2.0  -2.0
>      *    0  -2.0   3.0
>      */
>     double inva[9];
>     int s, i, j;
>     
>     gsl_matrix_view m
>          = gsl_matrix_view_array(a_data, 3, 3);
>     gsl_matrix_view inv
>          = gsl_matrix_view_array(inva,3,3);
>     gsl_permutation * p = gsl_permutation_alloc (3);
> 
>     printf("The matrix is\n");
>     for (i = 0; i < 3; ++i)
>         for (j = 0; j < 3; ++j)
>             printf(j==2?"%6.3f\n":"%6.3f ", gsl_matrix_get(&m.matrix,i,j));
> 
>     gsl_linalg_LU_decomp (&m.matrix, p, &s);    
>     gsl_linalg_LU_invert (&m.matrix, p, &inv.matrix);
> 
>     printf("The inverse is\n");
>     for (i = 0; i < 3; ++i)
>         for (j = 0; j < 3; ++j)
>             printf(j==2?"%6.3f\n":"%6.3f ",gsl_matrix_get(&inv.matrix,i,j));
> 
>     gsl_permutation_free (p);
>     return 0;
>     }
> --------------------
> 
> Compile and run:
> $ gcc gslinv.c -lgsl -lgslcblas -o gslinv
> $ ./gslinv
> The matrix is
>  1.000  0.600  0.000
>  0.000  1.500  1.000
>  0.000  1.000  1.000
> The inverse is
>  1.000 -1.200  1.200
> -0.000  2.000 -2.000
> -0.000 -2.000  3.000
> $
> 
> 
> 
> On Fri, Oct 31, 2008 at 2:22 PM, Philip Ogunbona
> <address@hidden <mailto:address@hidden>> wrote:
> 
>     While testing the linear algebra functions provided in gsl I came across
>     the following.
> 
>     With appropriate setup the call to,
> 
>     gsl_linalg_LU_invert(squareMatrix, perm, squareInverse);
> 
>     gave the inverse of the transpose of squareMatrix. I verified this using
>     octave. Note that squareMatrix is a simple 3x3 matrix of Gaussian
>     variates.
>     Is this the expected behaviour?
> 
>     Regards,
> 
>     Philip
> 
> 
>     _______________________________________________
>     Help-gsl mailing list
>     address@hidden <mailto:address@hidden>
>     http://lists.gnu.org/mailman/listinfo/help-gsl
> 
>