bug with dgesdd

[Bruno]

Hi,

in the octave bug list there is a matrix example where dgesdd performs badly
(see http://savannah.gnu.org/bugs/?func=deta ... 87#options).
I don't know if lapack team is aware of this problem (at least I didn't see it
on this forum).

I have transformed the matrix in a ascii format and written a small f90 program
which reads it and use dgesdd or dgesvd but I am unable to attach
the 2 files (what is possible to attach to a message for this forum ?) So I
have put the 2 files in my homepage :

http://www.iecn.u-nancy.fr/~pincon/

On the given matrix :

1/ dgesvd works quite well, I got :

|| U*diag(sigma)*VT - A ||_1 / || A ||_1 near 10^(-14)

2/ but with dgesdd, first the last singular value is negative !!! and
the previous relative error is about 48.4 !

I have tested using 2 differents blas lib and got the same results.

hth
Bruno

[admin]

Thanks Bruno!
We actually knew that there was a serious bug in dgesdd but we couldn't find a matrix to reproduce it.
For the current problems found in LAPACK, you can always refer to http://www.netlib.org/lapack/Errata/
The dgesdd problem is bug0025 as you can see.
The good news is that I was able to reproduce it and so I confirm the problem.
Your code and your matrix are really great and straightforward, so we will start debugging the dgesdd routine.
We keep you posted.
Merci beaucoup
Julie et Julien Langou

[Bruno]

The matrix has been posted on the octave bug tracker
by Jaroslav Hajek, I have only translated it in ascii and
wrote the small f90 code.

Thanks for the ref bug,

Bruno

[Julien Langou]

Hello,

Thanks for the bug report!

So I have taken a look at all this.

Good news: I think I have found where the bug is.
Bad news: I think the routine should have failed with INFO = 1.
So does not really help. But that helps somewhat!

So line 353 from 374 of DBDSC reads:

Code: Select all

            IF( ICOMPQ.EQ.2 ) THEN
               CALL DLASD0( NSIZE, SQRE, D( START ), E( START ),
     $                      U( START, START ), LDU, VT( START, START ),
     $                      LDVT, SMLSIZ, IWORK, WORK( WSTART ), INFO )
            ELSE
               CALL DLASDA( ICOMPQ, SMLSIZ, NSIZE, SQRE, D( START ),
     $                      E( START ), Q( START+( IU+QSTART-2 )*N ), N,
     $                      Q( START+( IVT+QSTART-2 )*N ),
     $                      IQ( START+K*N ), Q( START+( DIFL+QSTART-2 )*
     $                      N ), Q( START+( DIFR+QSTART-2 )*N ),
     $                      Q( START+( Z+QSTART-2 )*N ),
     $                      Q( START+( POLES+QSTART-2 )*N ),
     $                      IQ( START+GIVPTR*N ), IQ( START+GIVCOL*N ),
     $                      N, IQ( START+PERM*N ),
     $                      Q( START+( GIVNUM+QSTART-2 )*N ),
     $                      Q( START+( IC+QSTART-2 )*N ),
     $                      Q( START+( IS+QSTART-2 )*N ),
     $                      WORK( WSTART ), IWORK, INFO )
               IF( INFO.NE.0 ) THEN
                  RETURN
               END IF
            END IF


and I strongly believe the check on INFO should be outside the ELSE statement.

In your case, DBDSDC is calling DLASD0 twice (well with vecLib not with reference BLAS !!!), the first call returns an INFO = 1 so everything should stop. The second call returns INFO = 0 and overrides the first INFO. This is wrong.

So I'll send this to the lapackers and hopefully we will patch this soon.

Julien

[Bruno]

Thanks Julien. I have another question : when dgesdd "normally" outputs
with info=0, can we considered that the singular values and vectors have
near about the same accuray than with dgesvd ? (is there some paper
about the stability of dgessd ?)

Cordiales salutations
Bruno

[Julien Langou]

Hello, update: Last week, I forwarded your question around because I do not know the answer for sure. I did not receive any answer yet. Bottom line there should be not much difference in term of the quality of the answer. The first step for both is a reduction to bidiagonal form. The second step whether you use GESDD or GESVD has about the same accuracy or more accurate than the first. So at the end, you'll get the same accuracy more or less. Cheers, Julien.

[Bruno]

Thanks for the information Juien (and sorry for the delay I suppose
I forgot the dgesdd pb...).

Finally this means that for a singular elements computations one can
try first dgesdd and then dgesvd in case the first fails.

Regards
Bruno

[vukicevic]

Hi. I'm trying to find square root of the folowing matrix:

A =
(9,-16)
(-16,3)

I use dgesdd function, and I've got:

U =
(-0.769509, 0.638636)
(0.638636, 0.769509)

S =
(22.2788, 0)
(0, 10.2788)

VT =
(-0.769509, 0.638636)
(-0.638636, -0.769509)

Until now everything seems ok, because

U * S * VT = A =
(9,-16)
(-16,3)

Now it should be possible to get square root of A simply by taking square roots of eigenvalues in matrix S:

S**(1/2) =
(4.72004, 0)
(0, 3.20606)

and square root of A should be:

A**(1/2) = U * S**(1/2) * VT

But when I had tried to check the result, I've got:

A**(1/2) * A**(1/2) =
(17.3846, -5.89723)
(-5.89723, 15.1731)

which is obviously wrong.

I'm a litle confused with that result because the same procedure works fine with e.g.

B =
(2.2,0.4)
(0.4,2.8)

I suppose that I can try to solve the problem with dsyevr function, but it is interesting to know what is wrong.

Best regards, Vlatko

[Julien Langou]

You will not be able to compute the square root of a matrix using an SVD decomposition!