bug-gsl
[Top][All Lists]

## [Bug-gsl] LU_inverse bug

 From: Georges Szegezdi Subject: [Bug-gsl] LU_inverse bug Date: Mon, 11 Oct 2004 16:05:38 +0200 User-agent: Mozilla Thunderbird 0.8 (Windows/20040913)

Hi, I'm using the GSL library to do the pseudo-inverse calcul of a matrix (and then do some other stuff on it like transpose etc..).
```I have a problem with some type of matrices, I reproduce the example above.
```
When inverting a matrix, the result sometimes comes out to NaN or Inf. I tried both of the example above in Matlab and Octave and it worked, giving a warning tho. I think LU_invert should catch these errors(or impossible calculations) and reply with a very high number instead of a NaN, so we can continue the calculation with the inverted matrix. Or return an error code or something, I've looked around but haven't found anything about such return code(apologize me if that exists).
```Here is the code that I use:

```
gsl_matrix *Z = gsl_matrix_calloc (N,N); gsl_blas_dgemm(CblasTrans,CblasNoTrans,1.0,anotherMatrix,anotherMatrix,0.0,Z);
```gsl_matrix *inverse = gsl_matrix_calloc (N,N);
gsl_permutation *p = gsl_permutation_calloc(N);
int s;
gsl_linalg_LU_decomp(Z, p, &s);
gsl_linalg_LU_invert(Z, p, inverse);

The values and result of the matrices are:
Z:
```
1 0 0 0 0 0 0 0.988097 -0.821691 1.07345e-05 0.000741653
```0       0       0       0       0       0       0       0       0       0       0
```
0 0 3 0 0 0 0 2.97242 -1.80443 0.00179485 0.000916459 0 0 0 2 0 0 0 1.97459 -1.69657 0.00118993 0.00113353 0 0 0 0 1 0 0 0.981981 -0.803057 0.000930131 0.000391623
```0       0       0       0       0       0       0       0       0       0       0
```
0 0 0 0 0 0 2 1.91503 -1.38932 1.00122 1.00186 0.988097 0 2.97242 1.97459 0.981981 0 1.91503 8.67102 -6.39055 0.932407 0.99256 -0.821691 0 -1.80443 -1.69657 -0.803057 0 -1.38932 -6.39055 4.81479 -0.748694 -0.648082 1.07345e-05 0 0.00179485 0.00118993 0.000930131 0 1.00122 0.932407 -0.748694 1.00199 0.00115457 0.000741653 0 0.000916459 0.00113353 0.000391623 0 1.00186 0.99256 -0.648082 0.00115457 1.00186
```
inverse:
NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN
NaN
NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN
NaN
NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN
NaN
NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN
NaN
NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN     NaN
NaN
NaN     NaN     NaN     NaN     NaN     Inf     NaN     NaN     NaN     NaN
NaN
```
1.87121e+08 0 1.77284e+08 1.91665e+08 1.89854e+08 0 8.18867e+09 -1.40127e+08 5.18902e+07 -8.00467e+09 -8.00764e+09 -3.21066e+06 -0 -3.04114e+06 -3.28855e+06 -3.25729e+06 -0 -1.40127e+08 2.40206e+06 -893389 1.3697e+08 1.37021e+08 1.21369e+06 -0 1.14531e+06 1.2432e+06 1.23044e+06 -0 5.18902e+07 -893389 356305 -5.06987e+07 -5.072e+07 -1.82888e+08 0 -1.73278e+08 -1.87329e+08 -1.8556e+08 0 -8.00467e+09 1.3697e+08 -5.06987e+07 7.82482e+09 7.82773e+09 -1.82958e+08 0 -1.73344e+08 -1.87401e+08 -1.85631e+08 0 -8.00764e+09 1.37021e+08 -5.072e+07 7.82773e+09 7.83063e+09
```
```
If I add tiny little noise (about 0.0001) to the 0 and 1's, then Z and inverse are the following, and the result can be computed or used later on.
```
Z:
```
1.00196 0.000371977 0.00153351 0.00237999 0.00160393 0.000778553 0.00173078 0.145228 0
```.181509 0.00139875      0.001183
```
0.000371977 3.28104e-06 0.00235466 0.00119582 7.42628e-05 2.57241e-06 0.000684671 -1.93966e
```-05     0.000254357     8.50162e-05     0.000603617
```
0.00153351 0.00235466 3.00054 0.00249093 0.00166448 0.00132093 0.00275591 0.0597847 -
```0.0701054       0.00276492      0.00303812
```
0.00237999 0.00119582 0.00249093 2.00327 0.00152111 0.00140153 0.00159401 -0.404662 0
```.72661  0.00201412      0.00115836
```
0.00160393 7.42628e-05 0.00166448 0.00152111 1.00194 0.000226316 0.000708809 0.085846 0
```.101533 0.000619092     0.000642688
```
0.000778553 2.57241e-06 0.00132093 0.00140153 0.000226316 3.46488e-06 0.00125373 0.0001779
```1       0.000336786     0.000780943     0.000476505
```
0.00173078 0.000684671 0.00275591 0.00159401 0.000708809 0.00125373 2.00203 0.448053 -
```0.537788        1.00138 1.00199
```
0.145228 -1.93966e-05 0.0597847 -0.404662 0.085846 0.00017791 0.448053 0.220735-
```0.215039        0.289155        0.159107
```
0.181509 0.000254357 -0.0701054 0.72661 0.101533 0.000336786 -0.537788 -0.215039 0
```.517096 -0.129895       -0.407587
```
0.00139875 8.50162e-05 0.00276492 0.00201412 0.000619092 0.000780943 1.00138 0.289155 -
```0.129895        1.00097 0.00103852
```
0.001183 0.000603617 0.00303812 0.00115836 0.000642688 0.000476505 1.00199 0.159107 -
```0.407587        0.00103852      1.00168

inverse:
```
-2.11121e+14 1.37696e+17 -1.35357e+14 4.36086e+14 -1.15437e+14 -3.6991e+16 -5.1383e+14 1.76733e+
```15      -3.70722e+14    -2.75925e+13    1.72982e+13
```
1.37981e+17 -8.99931e+19 8.84663e+16 -2.85002e+17 7.54453e+16 2.41745e+19 3.37537e+17 -1.15504e
```+18     2.42288e+17     1.63131e+16     -1.30239e+16
```
-1.3788e+14 8.99291e+16 -8.84179e+13 2.84732e+14 -7.53876e+13 -2.41454e+16 -3.50761e+14 1.15403e+
```15      -2.42091e+14    -2.79554e+12    2.65045e+13
```
4.26692e+14 -2.78286e+17 2.73498e+14 -8.81619e+14 2.33318e+14 7.48099e+16 9.81945e+14 -3.57263e
```+15     7.49341e+14     1.12455e+14     2.16622e+13
```
-1.15066e+14 7.50469e+16 -7.37697e+13 2.37687e+14 -6.29163e+13 -2.01628e+16 -2.77811e+14 9.63263e+
```14      -2.02055e+14    -1.72821e+13    7.18693e+12
```
-3.52257e+16 2.29732e+19 -2.25715e+16 7.28093e+16 -1.92627e+16 -6.18101e+18 -7.51066e+16 2.95014e+
```17      -6.18707e+16    -1.5257e+16     -7.75471e+15
```
-2.59225e+15 1.69241e+18 -1.67716e+15 5.29793e+15 -1.41516e+15 -4.43568e+17 -1.88197e+16 2.15435e+
```16      -4.53382e+15    1.22034e+16     1.27402e+16
```
1.74131e+15 -1.13568e+18 1.11622e+15 -3.59751e+15 9.52145e+14 3.05232e+17 4.08137e+15 -1.45788e
```+16     3.05791e+15     3.84631e+14     1.41962e+13
```
-3.67718e+14 2.39828e+17 -2.35734e+14 7.5963e+14 -2.01065e+14 -6.44441e+16 -8.76872e+14 3.07846e+
```15      -6.45728e+14    -6.6191e+13     1.20177e+13
```
2.05606e+15 -1.34271e+18 1.33343e+15 -4.19025e+15 1.12197e+15 3.49594e+17 1.75463e+16 -1.70546e
```+16     3.59225e+15     -1.23051e+16    -1.27279e+16
```
2.09856e+15 -1.37043e+18 1.36067e+15 -4.27812e+15 1.14522e+15 3.57054e+17 1.76354e+16 -1.74106e
```+16     3.66691e+15     -1.22851e+16    -1.2717e+16

```
Here is a last example, I also had NaN on that matrix, I don't know why tho. this one doesn't use pure 0 or 1's.
```
Z:
297929  3805.55 -3825.03        3825.03 -3825.03        -3825.03        1241.64
3805.55 50.4217 -50.7088        50.7088 -50.7088        -50.7088        14.0719
-3825.03        -50.7088        51      -51     51      51      -13.7362
3825.03 50.7088 -51     51      -51     -51     13.7362
-3825.03        -50.7088        51      -51     51      51      -13.7362
-3825.03        -50.7088        51      -51     51      51      -13.7362
1241.64 14.0719 -13.7362        13.7362 -13.7362        -13.7362        102.012

inverse:
NaN     NaN     NaN     NaN     NaN     NaN     NaN
NaN     NaN     NaN     NaN     NaN     NaN     NaN
NaN     NaN     NaN     NaN     NaN     NaN     NaN
NaN     NaN     NaN     NaN     NaN     NaN     NaN
NaN     NaN     -Inf    NaN     Inf     NaN     NaN
```
5.81559e+12 -5.24055e+16 2.06566e+30 -0 -0 -2.06566e+30 2.82898e+14
```0       0       -2.81475e+14    0       0       2.81475e+14     0

I don't know if it is a real bug or a math issue with matrices, please let me
know.
I have looked into the mailing forum but haven't found anything about that
topic.
Thanks.
G

```