[Help-gsl] bug in gsl_linalg_LU_invert?

[kdsfinger]

hi, there
I am using a symmetric matrix sigma and the inverse matrix of that
sigma_inv, the dimension is 16x16. The multiplication of sigma and
sigma_inv should be 1 on diagonal and zero elsewhere. I am using
double and am expecting the epsilon of the result I get should be
small (< 1e-10) but actually some of them are very large (~0.001). The
magnitude of the value in the sigma_inv is about 1. Here is the part
of the code:
--------------------------------------------
//dimension = 16; covariance matrix sigma has been loaded and enforced
to be symmetric.
gsl_matrix* sigma_inv = gsl_matrix(dimension, dimension);
gsl_matrix_set_zero(sigma_inv);

for (uint i = 1; i < dimension; i++)
    for (uint j = 0; j < i; j++)
        gsl_matrix_set(sigma, i, j, gsl_matrix_get(sigma, j, i));
                                
gsl_matrix* work = gsl_matrix(dimension, dimension);
gsl_matrix_memcpy(work, sigma);
gsl_permutation *p = gsl_permutation_alloc(dimension);
gsl_permutation_init(p);
int s = 0;
gsl_linalg_LU_decomp( work, p, &s );
gsl_linalg_LU_invert( work, p, sigma_inv );

gsl_matrix* tmp = gsl_matrix_alloc(dimension, dimension);
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1, sigma, sigma_inv, 0, tmp);
--------------------------------------------

I've checked that the det(sigma) > 0.

Now I check the value in tmp and the difference of tmp(x,y) and
tmp(y,x) and I list a few of them here: 1st column is x, 2nd column is
y, 3rd is tmp(x,y), 4th is tmp(y,x), 5th is tmp(x,y)-tmp(y,x)

0 1 -0.0429039 0.000897743 -0.0438017 // is the difference too large?
0 2 0.00112287 0.00126048 -0.000137616
0 3 -1.67638e-16 0.00287359 -0.00287359 //this is very asymmetric
0 4 0.0639474 -0.00404521 0.0679926
12 15 1.79978e-17 -2.34188e-17 4.14165e-17 //1e-17 looks reasonable
13 14 5.3481e-16 6.27122e-17 4.72098e-16
13 15 1.93639e-16 -3.32037e-17 2.26842e-16
14 15 2.7616e-16 -5.70182e-16 8.46342e-16


The diagonal value is:
0.994551 1.00339 0.999861 1 1.03005 1 1 1 1 1 1 1 1 1 1 1

My question is: is the result in the normal error range? Also, since
the sigma is symmetric, I expect that sigma_inv should be also
symmetric (within small epsilon). I did the same thing to sigma_inv as
to tmp. The list is long and I only show a few here:

0 3 1.73063 1.72869 0.00193355
0 5 -3.41303 -3.47368 0.060645
4 11 2.84698 3.03546 -0.18848
4 12 0.758666 0.860237 -0.101571
5 14 -2.27969 -2.27969 -4.44089e-16
5 15 -2.1416 -2.1416 -4.44089e-16
6 7 -0.0664263 -0.0664263 -8.32667e-17
6 8 0.196091 0.196091 -1.94289e-16

Thanks for your comments.

zl2k