help-gsl
[Top][All Lists]

## Re: [Help-gsl] permutation matrix with lu decomposition problem

 From: Patrick Alken Subject: Re: [Help-gsl] permutation matrix with lu decomposition problem Date: Sat, 1 Jun 2019 15:00:08 +0000 User-agent: Mozilla/5.0 (X11; Linux x86_64; rv:60.0) Gecko/20100101 Thunderbird/60.7.0

```GSL defines the permutation p acting on a vector as
(https://www.gnu.org/software/gsl/doc/html/permutation.html#permutations):

vi' = P v = v_{pi}

So for p = (2,0,1) acting on [ v1 v2 v3]^T, we would obtain [ v3 v1 v2 ]

The matrix representation of this is:

[ 0 0 1 ]
[ 1 0 0 ]
[ 0 1 0 ]

so it seems things are working correctly. Of course the two matrices you
listed are just transposes (inverses) of each other, so it really just
depends on the convention of the permutation definition. But I think in
this case things are correct and consistent.

Patrick

On 5/31/19 11:39 PM, Elias Posoukidis wrote:
> Hi everyone,
>
> i have the following problem with the permutation matrix after a lu
> -decomposition. The following code
>
> //----------------------------------------------
>
>        double a_data[] = { 25, 5, 1, 64, 8, 1, 144, 12, 1 };
>
>         gsl_matrix_view m = gsl_matrix_view_array(a_data, 3, 3);
>
>         int s;
>
>         gsl_permutation *p = gsl_permutation_alloc(3);
>
>         gsl_linalg_LU_decomp(&m.matrix, p, &s);
>
>         for (size_t i = 0; i < 3; i++) {
>             for (size_t j = 0; j < 3; j++) {
>                 printf("%3.2f,", gsl_matrix_get(&m.matrix, i, j));
>             }
>             printf("\n");
>         }
>
>         gsl_permutation_fprintf(stdout, p, " %u");
>         printf("\n");
>         gsl_permutation_free(p);
>
> //----------------------------------------------
>
> gives for the lu-matrix
>
> 144.00,12.00,1.00,
> 0.17,2.92,0.83,
> 0.44,0.91,-0.20,
>
> which is the correct result, but for the permutation vector the
> following  {2 0 1}, which means that
>
> the permutation matrix is
>
> 0,1,0
>
> 0,0,1
>
> 1,0,0
>
> but the correct permutation matrix should be
>
> 0,0,1
>
> 1,0,0
>
> 0,1,0
>
> What am i doing wrong?
>
> kind regards,
>
> eliasp
>
>
>

```