[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
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
>
>
>