Re: [Help-gsl] gsl_linalg_LU_invert

[Warren Weckesser]

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>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
> http://lists.gnu.org/mailman/listinfo/help-gsl
>